solve_lin_eq.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOLVE_LIN_EQ                                            C
!  Purpose: Interface for linear equation solver                       C
!                                                                      C
!  Author: M. Syamlal                                 Date: 31-MAY-96  C
!  Reviewer:                                          Date:            C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE solve_lin_eq(VNAME, Vno, VAR, A_M, B_M, M, ITMAX,&
                              method, sweep, tol1, pc, ier)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE geometry
      USE compar
      USE residual
      USE toleranc
      USE leqsol
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! variable name
      CHARACTER(LEN=*), INTENT(IN) :: Vname
! variable number
!     Note: not really used beyond this subroutine. here it is
!     used for potentially adjusting the tolerances but it is
!     currently disabled code.
!     1 = pressure correction equation
!     2 = solids correction equation or gas/solids continuity
!     3 = gas/solids u-momentum
!     4 = gas/solids v-momentum
!     5 = gas/solids w-momentum
!     6 = temperature
!     7 = species
!     8 = granular temperature
!     9 = scalar, E_Turb_G, k_Turb_G
      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, 0:dimension_m)
! vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
! phase index
      INTEGER, INTENT(IN) :: M
! maximum number of iterations (generally leq_it)
      INTEGER, INTENT(IN) :: ITMAX
! linear equation solver method (generally leq_method)
!     1 = sor
!     2 = bicgstab (default)
!     3 = gmres
!     5 = cg
      INTEGER, INTENT(IN) :: METHOD
! sweep direction of leq solver (leq_sweep)
!     e.g., options = 'isis', 'rsrs' (default), 'asas'
      CHARACTER(LEN=4), INTENT(IN) :: SWEEP
! convergence tolerance for leq solver (leq_tol)
      DOUBLE PRECISION, INTENT(IN) :: TOL1
! preconditioner (leq_pc)
!     options = 'line' (default), 'diag', 'none'
      CHARACTER(LEN=4), INTENT(IN) :: PC
! error index
      INTEGER, INTENT(INOUT) :: IER
!-----------------------------------------------
! Local parameters
!-----------------------------------------------
! Adjust LEQ tolerance flag
      LOGICAL, PARAMETER :: adjust_leq_tol = .false.
      LOGICAL, PARAMETER :: leq_tol_scheme1 = .false.
! currently only used for gmres routine
      INTEGER, PARAMETER :: MAX_IT = 1
!-----------------------------------------------
! Local variables
!-----------------------------------------------
!
      DOUBLE PRECISION :: max_resid_local, tol_resid_max
! convergence tolerance for leq solver
      DOUBLE PRECISION :: TOL
! transpose of septadiaganol matrix A_M
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: A_MT
! indices
      INTEGER :: II, IJK
! for constructing local character strings
      CHARACTER(LEN=80) :: LINE0, LINE1
!-----------------------------------------------


! Adjusting the tolerances
! ---------------------------------------------------------------->>>
      IF(adjust_leq_tol) THEN
         max_resid_local = maxval(resid(:,m),1)
         tol_resid_max   = max(tol_resid, tol_resid_t, tol_resid_th, tol_resid_x)
         IF(leq_tol_scheme1.AND.resid(vno,m).LT.1.0d-1) THEN
            if(vno.le.5) then
               tol = max(tol1,tol1*resid(vno,m)/tol_resid)
            elseif (vno.eq.6) then
               tol = max(tol1,tol1*resid(vno,m)/tol_resid_t)
            elseif (vno.eq.7) then
               tol = max(tol1,tol1*resid(vno,m)/tol_resid_x)
            elseif (vno.eq.8) then
               tol = max(tol1,tol1*resid(vno,m)/tol_resid_th)
            endif
            Write(*,*) 'Adjusting LEQ_Tolerance', vname, tol, resid(vno,m)
         ELSEIF(max_resid_local.LT.1.0d-1) THEN
            tol = max(tol1,tol1*max_resid_local/tol_resid_max)
            Write(*,*) 'Adjusting LEQ_Tolerance', vname, tol, max_resid_local
         ENDIF
      ELSE
        tol = tol1
      ENDIF
! ----------------------------------------------------------------<<<


! Solve the linear system of equations
! ---------------------------------------------------------------->>>
      SELECT CASE (method)
      CASE (1)
! SOR: Successive Sver Relaxation method from Templates
        CALL leq_sor (vname, vno, var, a_m(:,:,m), b_m(:,m), &
                      itmax, ier)

      CASE (2)
! BICGSTAB: BIConjugate Gradients STabilized method
         IF(do_transpose) THEN  ! mfix.dat keyword default=false
            allocate( a_mt(-3:3, ijkstart3:ijkend3 ))
!!$omp parallel do private(ijk,ii)
            DO ijk=ijkstart3,ijkend3
               do ii=-3,3
                  a_mt(ii,ijk) = a_m(ijk,ii,m)
               enddo
            ENDDO
            call leq_bicgst(vname, vno, var, a_mt(:,:), b_m(:,m), &
                            sweep, tol, pc, itmax, ier)
            deallocate( a_mt )
         ELSE
            call leq_bicgs(vname, vno, var, a_m(:,:,m), b_m(:,m),&
                           sweep, tol, pc, itmax, ier)
         ENDIF


      CASE (3)
! GMRES: A Generalized Minimal RESidual Algorithm
         call leq_gmres(vname, vno, var, a_m(:,:,m), b_m(:,m),&
                        sweep, tol, itmax, max_it, ier)

      CASE (4)
! Mix:
         ier = 0
         call leq_bicgs(vname,vno, var, a_m(:,:,m), b_m(:,m), sweep,&
                       tol, pc, itmax, ier)
         IF (ier .eq. -2) THEN
            ier = 0
            print*,'calling leq_gmres', vname
            call leq_gmres(vname, vno, var, a_m(:,:,m), b_m(:,m),&
                           sweep, tol, itmax, max_it, ier)
         ENDIF


      CASE (5)
! CG: Conjugate Gradients
         call leq_cg(vname, vno, var, a_m(:,:,m), b_m(:,m), sweep,&
                     tol, pc, itmax, ier)

!     CASE (6) - Disabled
! LSOR: Line Successive Over Relaxation method
!       CALL LEQ_LSOR(VNAME, VAR, A_M(:,:,M), B_M(:,M), ITMAX, IER)


      CASE DEFAULT
         line0(1:14) = 'SOLVE_LIN_EQ: '
         line0(15:80)= vname
         WRITE(line1,'(A, I2, A)') &
             'Error: LEQ_METHOD = ', method, ' is invalid'
         CALL write_error(line0, line1, 1)
         CALL mfix_exit(mype)
      END SELECT
! ----------------------------------------------------------------<<<

      RETURN
      END SUBROUTINE solve_lin_eq