File: /nfs/home/0/users/jenkins/mfix.git/model/solve_scalar_eq.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 SUBROUTINE SOLVE_Scalar_EQ(IER)
21
22
23
24
25
26
27
28
29 USE param
30 USE param1
31 USE toleranc
32 USE run
33 USE physprop
34 USE geometry
35 USE fldvar
36 USE output
37 USE indices
38 USE drag
39 USE residual
40 USE ur_facs
41 USE pgcor
42 USE pscor
43 USE leqsol
44 USE bc
45 USE energy
46 USE rxns
47 USE scalars
48 Use ambm
49 Use tmp_array, S_p => Array1, S_c => Array2, EPs => Array3, VxGama => Array4
50 USE compar
51 USE mflux
52 USE functions
53 IMPLICIT NONE
54
55
56
57
58
59
60
61
62 INTEGER IER
63
64
65
66
67
68
69
70
71
72 INTEGER m
73
74
75 INTEGER n
76
77 DOUBLE PRECISION apo
78
79
80
81
82 DOUBLE PRECISION res1, mres1, num_res, den_res
83 INTEGER ires1
84
85
86 INTEGER IJK
87
88
89 INTEGER LEQM, LEQI
90
91 character(LEN=8) :: Vname
92
93
94
95 call lock_ambm
96 call lock_tmp_array
97
98 RESID(RESID_sc,0) = ZERO
99 NUM_RESID(RESID_sc,0) = ZERO
100 DEN_RESID(RESID_sc,0) = ZERO
101 MAX_RESID(RESID_sc,0) = ZERO
102 IJK_RESID(RESID_sc,0) = 0
103
104
105
106 DO N = 1, NScalar
107
108 M = Phase4Scalar(N)
109
110 CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
111
112 IF(M == 0) THEN
113
114 DO IJK = IJKSTART3, IJKEND3
115
116 IF (FLUID_AT(IJK)) THEN
117 APO = ROP_GO(IJK)*VOL(IJK)*ODT
118 S_P(IJK) = APO + ( ZMAX(SUM_R_G(IJK)) &
119 + Scalar_p(IJK, N) )*VOL(IJK)
120 S_C(IJK) = APO*ScalarO(IJK,N) &
121 + Scalar(IJK,N)*ZMAX((-SUM_R_G(IJK)))*VOL(IJK) &
122 + Scalar_c(IJK, N)*VOL(IJK)
123 ELSE
124
125 (IJK) = ZERO
126 S_C(IJK) = ZERO
127
128 ENDIF
129 END DO
130
131
132 IF(.NOT.ADDED_MASS) THEN
133 CALL CONV_DIF_PHI (Scalar(1,N), DIF_Scalar(1,N), DISCRETIZE(9), &
134 U_G, V_G, W_G, Flux_gE, Flux_gN, Flux_gT, M, A_M, B_M, IER)
135 ELSE
136 CALL CONV_DIF_PHI (Scalar(1,N), DIF_Scalar(1,N), DISCRETIZE(9), &
137 U_G, V_G, W_G, Flux_gSE, Flux_gSN, Flux_gST, M, A_M, B_M, IER)
138 ENDIF
139
140
141 CALL BC_PHI (Scalar(1,N), BC_Scalar(1,N), BC_ScalarW(1,N), BC_HW_Scalar(1,N), &
142 BC_C_Scalar(1,N), M, A_M, B_M, IER)
143
144
145 CALL SOURCE_PHI (S_P, S_C, EP_G, Scalar(1,N), M, A_M, B_M, IER)
146
147 CALL CALC_RESID_S (Scalar(1,N), A_M, B_M, M, num_res, den_res, res1, &
148 mres1, ires1, ZERO, IER)
149 RESID(RESID_sc,0) = RESID(RESID_sc,0)+res1
150 NUM_RESID(RESID_sc,0) = NUM_RESID(RESID_sc,0)+num_res
151 DEN_RESID(RESID_sc,0) = DEN_RESID(RESID_sc,0)+den_res
152 if(mres1 .gt. MAX_RESID(RESID_sc,0))then
153 MAX_RESID(RESID_sc,0) = mres1
154 IJK_RESID(RESID_sc,0) = ires1
155 endif
156
157 CALL UNDER_RELAX_S (Scalar(1,N), A_M, B_M, M, UR_FAC(9), IER)
158
159
160
161
162
163
164
165
166
167 CALL ADJUST_LEQ (res1, LEQ_IT(9), LEQ_METHOD(9), &
168 LEQI, LEQM, IER)
169
170 write(Vname, '(A,I2)')'Scalar',N
171 CALL SOLVE_LIN_EQ (Vname, 9, Scalar(1,N), A_M, B_M, M, LEQI, LEQM, &
172 LEQ_SWEEP(9), LEQ_TOL(9), LEQ_PC(9), IER)
173
174
175 ELSE
176
177 DO IJK = IJKSTART3, IJKEND3
178
179 IF (FLUID_AT(IJK)) THEN
180 APO = ROP_sO(IJK, M)*VOL(IJK)*ODT
181 S_P(IJK) = APO + ( ZMAX(SUM_R_s(IJK, M)) &
182 + Scalar_p(IJK, N) )*VOL(IJK)
183 S_C(IJK) = APO*ScalarO(IJK,N) &
184 + Scalar(IJK,N)*ZMAX((-SUM_R_s(IJK, M)))*VOL(IJK)&
185 + Scalar_c(IJK, N)*VOL(IJK)
186 EPs(IJK) = EP_s(IJK, M)
187 ELSE
188
189 (IJK) = ZERO
190 S_C(IJK) = ZERO
191 EPS(IJK) = ZERO
192
193 ENDIF
194 END DO
195
196
197 IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
198 CALL CONV_DIF_PHI (Scalar(1,N), DIF_Scalar(1,N), DISCRETIZE(9), &
199 U_s(1,m), V_s(1,m), W_s(1,m), Flux_sE(1,M), Flux_sN(1,M), Flux_sT(1,M), M, &
200 A_M, B_M, IER)
201 ELSE
202 CALL CONV_DIF_PHI (Scalar(1,N), DIF_Scalar(1,N), DISCRETIZE(9), &
203 U_s(1,m), V_s(1,m), W_s(1,m), Flux_sSE, Flux_sSN, Flux_sST, M, &
204 A_M, B_M, IER)
205 ENDIF
206
207
208 CALL BC_PHI (Scalar(1,N), BC_Scalar(1,N), BC_ScalarW(1,N), BC_HW_Scalar(1,N), &
209 BC_C_Scalar(1,N), M, A_M, B_M, IER)
210
211
212 CALL SOURCE_PHI (S_P, S_C, EPs, Scalar(1,N), M, A_M, B_M, IER)
213
214 CALL CALC_RESID_S (Scalar(1,N), A_M, B_M, M, num_res, den_res, res1, &
215 mres1, ires1, ZERO, IER)
216 RESID(RESID_sc,0) = RESID(RESID_sc,0)+res1
217 NUM_RESID(RESID_sc,0) = NUM_RESID(RESID_sc,0)+num_res
218 DEN_RESID(RESID_sc,0) = DEN_RESID(RESID_sc,0)+den_res
219 if(mres1 .gt. MAX_RESID(RESID_sc,0))then
220 MAX_RESID(RESID_sc,0) = mres1
221 IJK_RESID(RESID_sc,0) = ires1
222 endif
223
224 CALL UNDER_RELAX_S (Scalar(1,N), A_M, B_M, M, UR_FAC(9), IER)
225
226
227
228
229
230
231
232
233 CALL ADJUST_LEQ (res1, LEQ_IT(9), LEQ_METHOD(9), &
234 LEQI, LEQM, IER)
235
236 write(Vname, '(A,I2)')'Scalar',N
237 CALL SOLVE_LIN_EQ (Vname, 9, Scalar(1,N), A_M, B_M, M, LEQI, LEQM, &
238 LEQ_SWEEP(9), LEQ_TOL(9), LEQ_PC(9), IER)
239
240
241 END IF
242 END DO
243
244 call unlock_ambm
245 call unlock_tmp_array
246
247 RETURN
248 END SUBROUTINE SOLVE_Scalar_EQ
249
250
251
252
253
254