File: /nfs/home/0/users/jenkins/mfix.git/model/solve_energy_eq.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 SUBROUTINE SOLVE_ENERGY_EQ(IER)
23
24
25
26
27 USE param
28 USE param1
29 USE toleranc
30 USE run
31 USE physprop
32 USE geometry
33 USE fldvar
34 USE output
35 USE indices
36 USE drag
37 USE residual
38 USE ur_facs
39 USE pgcor
40 USE pscor
41 USE leqsol
42 USE bc
43 USE energy
44 USE rxns
45 Use ambm
46 Use tmp_array, S_p => ARRAY1, S_C => ARRAY2, EPs => ARRAY3
47 Use tmp_array1, VxGama => ARRAYm1
48 USE compar
49 USE discretelement
50 USE des_thermo
51 USE mflux
52 USE mpi_utility
53 USE sendrecv
54 USE ps
55 USE mms
56 USE functions
57
58 IMPLICIT NONE
59
60
61
62
63 INTEGER, INTENT(INOUT) :: IER
64
65
66
67
68 INTEGER :: M
69 INTEGER :: TMP_SMAX
70
71 DOUBLE PRECISION :: CpxFlux_E(DIMENSION_3), &
72 CpxFlux_N(DIMENSION_3), &
73 CpxFlux_T(DIMENSION_3)
74
75 DOUBLE PRECISION :: apo
76
77 INTEGER :: IJK, I, J, K
78
79 INTEGER :: LEQM, LEQI
80
81
82
83
84
85 INTEGER :: Err_l(0:numPEs-1)
86 INTEGER :: Err_g(0:numPEs-1)
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108 call lock_ambm
109 call lock_tmp_array
110 call lock_tmp_array1
111
112
113 = 0
114
115 TMP_SMAX = SMAX
116 IF(DISCRETE_ELEMENT) THEN
117 TMP_SMAX = 0
118 ENDIF
119
120
121 DO M = 0, TMP_SMAX
122 CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
123 ENDDO
124
125 DO IJK = IJKSTART3, IJKEND3
126 I = I_OF(IJK)
127 J = J_OF(IJK)
128 K = K_OF(IJK)
129
130 IF (IS_ON_myPE_plus2layers(IP1(I),J,K)) THEN
131 IF(.NOT.ADDED_MASS) THEN
132 CpxFlux_E(IJK) = HALF * (C_pg(IJK) + C_pg(IP_OF(IJK))) * Flux_gE(IJK)
133 ELSE
134 CpxFlux_E(IJK) = HALF * (C_pg(IJK) + C_pg(IP_OF(IJK))) * Flux_gSE(IJK)
135 ENDIF
136 ENDIF
137
138 IF (IS_ON_myPE_plus2layers(I,JP1(J),K)) THEN
139 IF(.NOT.ADDED_MASS) THEN
140 CpxFlux_N(IJK) = HALF * (C_pg(IJK) + C_pg(JP_OF(IJK))) * Flux_gN(IJK)
141 ELSE
142 CpxFlux_N(IJK) = HALF * (C_pg(IJK) + C_pg(JP_OF(IJK))) * Flux_gSN(IJK)
143 ENDIF
144 ENDIF
145
146 IF (IS_ON_myPE_plus2layers(I,J,KP1(K))) THEN
147 IF(.NOT.ADDED_MASS) THEN
148 CpxFlux_T(IJK) = HALF * (C_pg(IJK) + C_pg(KP_OF(IJK))) * Flux_gT(IJK)
149 ELSE
150 CpxFlux_T(IJK) = HALF * (C_pg(IJK) + C_pg(KP_OF(IJK))) * Flux_gST(IJK)
151 ENDIF
152 ENDIF
153
154 IF (FLUID_AT(IJK)) THEN
155 APO = ROP_GO(IJK)*C_PG(IJK)*VOL(IJK)*ODT
156 S_P(IJK) = APO + S_RPG(IJK)*VOL(IJK)
157 S_C(IJK) = APO*T_GO(IJK)-HOR_G(IJK)*VOL(IJK)+S_RCG(IJK)*VOL(IJK)
158 IF(USE_MMS) S_C(IJK) = S_C(IJK) + MMS_T_G_SRC(IJK)*VOL(IJK)
159 ELSE
160 S_P(IJK) = ZERO
161 S_C(IJK) = ZERO
162 ENDIF
163 ENDDO
164
165
166 IF(DES_CONTINUUM_COUPLED) CALL DES_Hgm(S_C, S_P)
167
168
169 CALL CONV_DIF_PHI (T_g, K_G, DISCRETIZE(6), U_G, V_G, W_G, &
170 CpxFlux_E, CpxFlux_N, CpxFlux_T, 0, A_M, B_M, IER)
171
172
173 CALL BC_PHI (T_g, BC_T_G, BC_TW_G, BC_HW_T_G, BC_C_T_G, 0, A_M, B_M, IER)
174
175
176 CALL SOURCE_PHI (S_P, S_C, EP_G, T_G, 0, A_M, B_M, IER)
177
178
179 IF(POINT_SOURCE) CALL POINT_SOURCE_PHI (T_g, PS_T_g, &
180 PS_CpxMFLOW_g, 0, A_M, B_M, IER)
181
182 DO M = 1, TMP_SMAX
183 DO IJK = IJKSTART3, IJKEND3
184 I = I_OF(IJK)
185 J = J_OF(IJK)
186 K = K_OF(IJK)
187
188 IF (IS_ON_myPE_plus2layers(IP1(I),J,K)) THEN
189 IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
190 CpxFlux_E(IJK) = HALF * (C_ps(IJK,M) + C_ps(IP_OF(IJK),M)) * Flux_sE(IJK,M)
191 ELSE
192 (IJK) = HALF * (C_ps(IJK,M) + C_ps(IP_OF(IJK),M)) * Flux_sSE(IJK)
193 ENDIF
194 ENDIF
195
196 IF (IS_ON_myPE_plus2layers(I,JP1(J),K)) THEN
197 IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
198 CpxFlux_N(IJK) = HALF * (C_ps(IJK,M) + C_ps(JP_OF(IJK),M)) * Flux_sN(IJK,M)
199 ELSE
200 CpxFlux_N(IJK) = HALF * (C_ps(IJK,M) + C_ps(JP_OF(IJK),M)) * Flux_sSN(IJK)
201 ENDIF
202 ENDIF
203
204 IF (IS_ON_myPE_plus2layers(I,J,KP1(K))) THEN
205 IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
206 CpxFlux_T(IJK) = HALF * (C_ps(IJK,M) + C_ps(KP_OF(IJK),M)) * Flux_sT(IJK,M)
207 ELSE
208 CpxFlux_T(IJK) = HALF * (C_ps(IJK,M) + C_ps(KP_OF(IJK),M)) * Flux_sST(IJK)
209 ENDIF
210 ENDIF
211
212 IF (FLUID_AT(IJK)) THEN
213 APO = ROP_SO(IJK,M)*C_PS(IJK,M)*VOL(IJK)*ODT
214 S_P(IJK) = APO + S_RPS(IJK,M)*VOL(IJK)
215 S_C(IJK) = APO*T_SO(IJK,M) - HOR_S(IJK,M)*VOL(IJK) + &
216 S_RCS(IJK,M)*VOL(IJK)
217 VXGAMA(IJK,M) = GAMA_GS(IJK,M)*VOL(IJK)
218 EPS(IJK) = EP_S(IJK,M)
219 IF(USE_MMS) S_C(IJK) = S_C(IJK) + MMS_T_S_SRC(IJK)*VOL(IJK)
220 ELSE
221 S_P(IJK) = ZERO
222 S_C(IJK) = ZERO
223 VXGAMA(IJK,M) = ZERO
224 EPS(IJK) = ZERO
225 IF(USE_MMS) EPS(IJK) = EP_S(IJK,M)
226 ENDIF
227 ENDDO
228
229
230 CALL CONV_DIF_PHI (T_s(1,M), K_S(1,M), DISCRETIZE(6), &
231 U_S(1,M), V_S(1,M), W_S(1,M), CpxFlux_E, CpxFlux_N, &
232 CpxFlux_T, M, A_M, B_M, IER)
233
234
235 CALL BC_PHI (T_s(1,M), BC_T_S(1,M), BC_TW_S(1,M), &
236 BC_HW_T_S(1,M), BC_C_T_S(1,M), M, A_M, B_M, IER)
237
238
239 CALL SOURCE_PHI (S_P, S_C, EPS, T_S(1,M), M, A_M, B_M, IER)
240
241
242 IF(POINT_SOURCE) CALL POINT_SOURCE_PHI (T_s(:,M), PS_T_s(:,M),&
243 PS_CpxMFLOW_s(:,M), M, A_M, B_M, IER)
244
245 ENDDO
246
247
248 IF (TMP_SMAX > 0 .AND. .NOT.USE_MMS) &
249 CALL PARTIAL_ELIM_S (T_G, T_S, VXGAMA, A_M, B_M, IER)
250
251 CALL CALC_RESID_S (T_G, A_M, B_M, 0, NUM_RESID(RESID_T,0),&
252 DEN_RESID(RESID_T,0), RESID(RESID_T,0), MAX_RESID(RESID_T,&
253 0), IJK_RESID(RESID_T,0), ZERO, IER)
254
255 CALL UNDER_RELAX_S (T_G, A_M, B_M, 0, UR_FAC(6), IER)
256
257
258
259
260
261
262
263
264 DO M = 1, TMP_SMAX
265 CALL CALC_RESID_S (T_S(1,M), A_M, B_M, M, NUM_RESID(RESID_T,M), &
266 DEN_RESID(RESID_T,M), RESID(RESID_T,M), MAX_RESID(&
267 RESID_T,M), IJK_RESID(RESID_T,M), ZERO, IER)
268
269 CALL UNDER_RELAX_S (T_S(1,M), A_M, B_M, M, UR_FAC(6), IER)
270 ENDDO
271
272
273 CALL ADJUST_LEQ(RESID(RESID_T,0), LEQ_IT(6), LEQ_METHOD(6), &
274 LEQI, LEQM, IER)
275
276
277 CALL SOLVE_LIN_EQ ('T_g', 6, T_G, A_M, B_M, 0, LEQI, LEQM, &
278 LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6), IER)
279
280 IF(ier == -2) Err_l(myPE) = 120
281
282
283 DO IJK = IJKSTART3, IJKEND3
284 IF(.NOT.WALL_AT(IJK))&
285 T_g(IJK) = MIN(TMAX, MAX(TMIN, T_g(IJK)))
286 ENDDO
287
288
289
290 DO M = 1, TMP_SMAX
291 CALL ADJUST_LEQ (RESID(RESID_T,M), LEQ_IT(6), LEQ_METHOD(6), &
292 LEQI, LEQM, IER)
293
294 CALL SOLVE_LIN_EQ ('T_s', 6, T_S(1,M), A_M, B_M, M, LEQI, &
295 LEQM, LEQ_SWEEP(6), LEQ_TOL(6), LEQ_PC(6), IER)
296
297
298 IF(ier == -2) Err_l(myPE) = 121
299
300
301 DO IJK = IJKSTART3, IJKEND3
302 IF(.NOT.WALL_AT(IJK))&
303 T_s(IJK, M) = MIN(TMAX, MAX(TMIN, T_s(IJK, M)))
304 ENDDO
305 ENDDO
306
307 call unlock_ambm
308 call unlock_tmp_array
309 call unlock_tmp_array1
310
311
312
313
314
315 CALL global_all_sum(Err_l, Err_g)
316 IER = maxval(Err_g)
317
318
319 RETURN
320 END SUBROUTINE SOLVE_ENERGY_EQ
321