leq_sor.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: LEQ_SOR                                                 C
!  Purpose: Solve system of linear system using SOR method             C
!           Successive over-relaxation                                 C
!                                                                      C
!  Author: M. Syamlal                                 Date: 19-AUG-96  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE leq_sor(VNAME, VNO, VAR, A_M, B_M, ITMAX, IER)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE geometry
      USE indices
      USE compar
      USE sendrecv
      USE leqsol
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! variable name
      CHARACTER(LEN=*), INTENT(IN) :: Vname
! variable number (not really used here; see calling subroutine)
      INTEGER, INTENT(IN) :: VNO
! variable
!     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
!     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
!     e_Turb_G
      DOUBLE PRECISION, INTENT(INOUT) :: Var(dimension_3)
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: &
                        A_m(dimension_3, -3:3)
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: &
                        B_m(dimension_3)
! maximum number of iterations (generally leq_it)
      INTEGER, INTENT(IN) :: ITMAX
! error indicator
      INTEGER, INTENT(INOUT) :: IER
!-----------------------------------------------
! Local parameters
!-----------------------------------------------
! OVERRELAXATION FACTOR
      DOUBLE PRECISION, PARAMETER :: OMEGA = 1.0  !1.2
      integer :: iidebug
      parameter( iidebug = 0 )
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! Variable
      DOUBLE PRECISION :: Var_tmp(dimension_3)
! Indices
      INTEGER :: IJK
      INTEGER :: ITER

      DOUBLE PRECISION oAm
!-----------------------------------------------

!!$omp parallel do private(IJK,OAM)
      DO ijk = ijkstart3, ijkend3
         IF(.NOT.is_on_mype_owns(i_of(ijk),j_of(ijk), k_of(ijk))) cycle

         oam = one/a_m(ijk,0)
         a_m(ijk,0) = one
         a_m(ijk,-2) = a_m(ijk,-2)*oam
         a_m(ijk,-1) = a_m(ijk,-1)*oam
         a_m(ijk,1) = a_m(ijk,1)*oam
         a_m(ijk,2) = a_m(ijk,2)*oam
         a_m(ijk,-3) = a_m(ijk,-3)*oam
         a_m(ijk,3) = a_m(ijk,3)*oam
         b_m(ijk) = b_m(ijk)*oam
      ENDDO

      DO iter = 1, itmax
         IF (do_k) THEN

!!$omp parallel do private(IJK)
            DO ijk = ijkstart3, ijkend3
              IF(.NOT.is_on_mype_owns(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
              var_tmp(ijk) = var(ijk) + omega*(b_m(ijk)-&
                 a_m(ijk,-1)*var(im_of(ijk))-a_m(ijk,1)*var(ip_of(ijk))-&
                 a_m(ijk,-2)*var(jm_of(ijk))-a_m(ijk,2)*var(jp_of(ijk))-&
                 a_m(ijk,-3)*var(km_of(ijk))-a_m(ijk,3)*var(kp_of(ijk))-&
                 var(ijk))
            ENDDO
         ELSE

!!$omp parallel do private(IJK)
           DO ijk = ijkstart3, ijkend3
             IF(.NOT.is_on_mype_owns(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
             var_tmp(ijk) = var(ijk) + omega*(b_m(ijk)-&
                  a_m(ijk,-2)*var(jm_of(ijk))-a_m(ijk,2)*var(jp_of(ijk))-&
                  a_m(ijk,-1)*var(im_of(ijk))-a_m(ijk,1)*var(ip_of(ijk))-&
                  var(ijk))
           ENDDO
         ENDIF

      call send_recv(var,2)
      ENDDO

!!$omp parallel do private(IJK)
      DO ijk = ijkstart3, ijkend3
        var(ijk) = var_tmp(ijk)
      ENDDO

      iter_tot(vno) = iter


      RETURN
      END SUBROUTINE leq_sor