File: /nfs/home/0/users/jenkins/mfix.git/model/solve_lin_eq.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOLVE_LIN_EQ                                            C
4     !  Purpose: Interface for linear equation solver                       C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 31-MAY-96  C
7     !  Reviewer:                                          Date:            C
8     !                                                                      C
9     !                                                                      C
10     !  Literature/Document References:                                     C
11     !  Variables referenced:                                               C
12     !  Variables modified:                                                 C
13     !  Local variables:                                                    C
14     !                                                                      C
15     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
16     
17           SUBROUTINE SOLVE_LIN_EQ(VNAME, Vno, VAR, A_M, B_M, M, ITMAX,&
18                                   METHOD, SWEEP, TOL1, PC, IER)
19     
20     !-----------------------------------------------
21     ! Modules
22     !-----------------------------------------------
23           USE param
24           USE param1
25           USE geometry
26           USE compar
27           USE residual
28           USE toleranc
29           USE leqsol
30           IMPLICIT NONE
31     !-----------------------------------------------
32     ! Dummy arguments
33     !-----------------------------------------------
34     ! variable name
35           CHARACTER(LEN=*), INTENT(IN) :: Vname
36     ! variable number
37     !     Note: not really used beyond this subroutine. here it is
38     !     used for potentially adjusting the tolerances but it is
39     !     currently disabled code.
40     !     1 = pressure correction equation
41     !     2 = solids correction equation or gas/solids continuity
42     !     3 = gas/solids u-momentum
43     !     4 = gas/solids v-momentum
44     !     5 = gas/solids w-momentum
45     !     6 = temperature
46     !     7 = species
47     !     8 = granular temperature
48     !     9 = scalar, E_Turb_G, k_Turb_G
49           INTEGER, INTENT(IN) :: Vno
50     ! variable
51     !     e.g., pp_g, epp, rop_g, rop_s, u_g, u_s, v_g, v_s, w_g,
52     !     w_s, T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
53     !     e_Turb_G
54           DOUBLE PRECISION, INTENT(INOUT) :: Var(DIMENSION_3)
55     ! septadiagonal matrix A_m
56           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
57     ! vector b_m
58           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
59     ! phase index
60           INTEGER, INTENT(IN) :: M
61     ! maximum number of iterations (generally leq_it)
62           INTEGER, INTENT(IN) :: ITMAX
63     ! linear equation solver method (generally leq_method)
64     !     1 = sor
65     !     2 = bicgstab (default)
66     !     3 = gmres
67     !     5 = cg
68           INTEGER, INTENT(IN) :: METHOD
69     ! sweep direction of leq solver (leq_sweep)
70     !     e.g., options = 'isis', 'rsrs' (default), 'asas'
71           CHARACTER(LEN=4), INTENT(IN) :: SWEEP
72     ! convergence tolerance for leq solver (leq_tol)
73           DOUBLE PRECISION, INTENT(IN) :: TOL1
74     ! preconditioner (leq_pc)
75     !     options = 'line' (default), 'diag', 'none'
76           CHARACTER(LEN=4), INTENT(IN) :: PC
77     ! error index
78           INTEGER, INTENT(INOUT) :: IER
79     !-----------------------------------------------
80     ! Local parameters
81     !-----------------------------------------------
82     ! Adjust LEQ tolerance flag
83           LOGICAL, PARAMETER :: adjust_leq_tol = .FALSE.
84           LOGICAL, PARAMETER :: leq_tol_scheme1 = .FALSE.
85     ! currently only used for gmres routine
86           INTEGER, PARAMETER :: MAX_IT = 1
87     !-----------------------------------------------
88     ! Local variables
89     !-----------------------------------------------
90     !
91           DOUBLE PRECISION :: max_resid_local, tol_resid_max
92     ! convergence tolerance for leq solver
93           DOUBLE PRECISION :: TOL
94     ! transpose of septadiaganol matrix A_M
95           DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: A_MT
96     ! indices
97           INTEGER :: II, IJK
98     ! for constructing local character strings
99           CHARACTER(LEN=80) :: LINE0, LINE1
100     !-----------------------------------------------
101     
102     
103     ! Adjusting the tolerances
104     ! ---------------------------------------------------------------->>>
105           IF(adjust_leq_tol) THEN
106              max_resid_local = maxval(resid(:,M),1)
107              tol_resid_max   = max(TOL_RESID, TOL_RESID_T, TOL_RESID_TH, TOL_RESID_X)
108              IF(leq_tol_scheme1.AND.resid(Vno,M).LT.1.0D-1) THEN
109                 if(Vno.le.5) then
110                    TOL = MAX(TOL1,TOL1*RESID(Vno,M)/TOL_RESID)
111                 elseif (Vno.eq.6) then
112                    TOL = MAX(TOL1,TOL1*RESID(Vno,M)/TOL_RESID_T)
113                 elseif (Vno.eq.7) then
114                    TOL = MAX(TOL1,TOL1*RESID(Vno,M)/TOL_RESID_X)
115                 elseif (Vno.eq.8) then
116                    TOL = MAX(TOL1,TOL1*RESID(Vno,M)/TOL_RESID_Th)
117                 endif
118                 Write(*,*) 'Adjusting LEQ_Tolerance', Vname, tol, resid(Vno,M)
119              ELSEIF(max_resid_local.LT.1.0D-1) THEN
120                 TOL = MAX(TOL1,TOL1*max_resid_local/TOL_RESID_max)
121                 Write(*,*) 'Adjusting LEQ_Tolerance', Vname, tol, max_resid_local
122              ENDIF
123           ELSE
124             TOL = TOL1
125           ENDIF
126     ! ----------------------------------------------------------------<<<
127     
128     
129     ! Solve the linear system of equations
130     ! ---------------------------------------------------------------->>>
131           SELECT CASE (METHOD)
132           CASE (1)
133     ! SOR: Successive Sver Relaxation method from Templates
134             CALL LEQ_SOR (VNAME, VNO, VAR, A_M(:,:,M), B_M(:,M), &
135                           ITMAX, IER)
136     
137           CASE (2)
138     ! BICGSTAB: BIConjugate Gradients STabilized method
139              IF(do_transpose) THEN  ! mfix.dat keyword default=false
140                 allocate( A_mt(-3:3, ijkstart3:ijkend3 ))
141     !!$omp parallel do private(ijk,ii)
142                 DO ijk=ijkstart3,ijkend3
143                    do ii=-3,3
144                       A_mt(ii,ijk) = A_m(ijk,ii,M)
145                    enddo
146                 ENDDO
147                 call leq_bicgst(VNAME, VNO, VAR, A_Mt(:,:), B_M(:,M), &
148                                 SWEEP, TOL, PC, ITMAX, IER)
149                 deallocate( A_mt )
150              ELSE
151                 call leq_bicgs(VNAME, VNO, VAR, A_M(:,:,M), B_M(:,M),&
152                                SWEEP, TOL, PC, ITMAX, IER)
153              ENDIF
154     
155     
156           CASE (3)
157     ! GMRES: A Generalized Minimal RESidual Algorithm
158              call leq_gmres(VNAME, VNO, VAR, A_M(:,:,M), B_M(:,M),&
159                             SWEEP, TOL, ITMAX, MAX_IT, IER)
160     
161           CASE (4)
162     ! Mix:
163              IER = 0
164              call leq_bicgs(VNAME,VNO, VAR, A_M(:,:,M), B_M(:,M), SWEEP,&
165                            TOL, PC, ITMAX, IER)
166              IF (IER .eq. -2) THEN
167                 IER = 0
168                 print*,'calling leq_gmres', Vname
169                 call leq_gmres(VNAME, VNO, VAR, A_M(:,:,M), B_M(:,M),&
170                                SWEEP, TOL, ITMAX, MAX_IT, IER)
171              ENDIF
172     
173     
174           CASE (5)
175     ! CG: Conjugate Gradients
176              call leq_cg(VNAME, VNO, VAR, A_M(:,:,M), B_M(:,M), SWEEP,&
177                          TOL, PC, ITMAX, IER)
178     
179     !     CASE (6) - Disabled
180     ! LSOR: Line Successive Over Relaxation method
181     !       CALL LEQ_LSOR(VNAME, VAR, A_M(:,:,M), B_M(:,M), ITMAX, IER)
182     
183     
184           CASE DEFAULT
185              LINE0(1:14) = 'SOLVE_LIN_EQ: '
186              LINE0(15:80)= VName
187              WRITE(LINE1,'(A, I2, A)') &
188                  'Error: LEQ_METHOD = ', METHOD, ' is invalid'
189              CALL WRITE_ERROR(LINE0, LINE1, 1)
190              CALL mfix_exit(myPE)
191           END SELECT
192     ! ----------------------------------------------------------------<<<
193     
194           RETURN
195           END SUBROUTINE SOLVE_LIN_EQ
196