File: /nfs/home/0/users/jenkins/mfix.git/model/solve_granular_energy.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 SUBROUTINE SOLVE_GRANULAR_ENERGY(IER)
21
22
23
24
25 USE param
26 USE param1
27 USE toleranc
28 USE run
29 USE physprop
30 USE visc_s
31 USE geometry
32 USE fldvar
33 USE constant
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 compar
48 USE mflux
49 USE mpi_utility
50 USE mms
51 USE functions
52 IMPLICIT NONE
53
54
55
56
57 INTEGER :: IER
58
59
60
61
62 INTEGER :: M, L
63
64 DOUBLE PRECISION CpxFlux_E(DIMENSION_3), CpxFlux_N(DIMENSION_3), CpxFlux_T(DIMENSION_3)
65
66 DOUBLE PRECISION :: apo
67
68
69 DOUBLE PRECISION :: sourcelhs, sourcerhs
70
71 INTEGER :: IJK
72
73 INTEGER :: LEQM, LEQI
74
75 DOUBLE PRECISION :: M_PM,D_PM
76
77 DOUBLE PRECISION :: TOT_EPS(DIMENSION_3)
78
79 DOUBLE PRECISION :: TOT_SUM_RS(DIMENSION_3)
80
81 DOUBLE PRECISION :: TOT_NO(DIMENSION_3)
82
83 DOUBLE PRECISION :: smallTheta
84
85
86 DOUBLE PRECISION :: VxTC_ss(DIMENSION_3,DIMENSION_LM)
87
88
89
90
91
92 INTEGER :: Err_l(0:numPEs-1)
93 INTEGER :: Err_g(0:numPEs-1)
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111 call lock_ambm
112 call lock_tmp_array
113
114
115
116 = 0
117
118 smallTheta = (to_SI)**4 * ZERO_EP_S
119
120 DO M = 0, MMAX
121 CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
122 ENDDO
123
124 SELECT CASE (KT_TYPE_ENUM)
125 CASE(LUN_1984, AHMADI_1995, SIMONIN_1996, GD_1999, GTSH_2012)
126
127 DO M = 1, SMAX
128
129 DO IJK = ijkstart3, ijkend3
130
131 IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
132 CpxFlux_E(IJK) = 1.5D0 * Flux_sE(IJK,M)
133 CpxFlux_N(IJK) = 1.5D0 * Flux_sN(IJK,M)
134 CpxFlux_T(IJK) = 1.5D0 * Flux_sT(IJK,M)
135 ELSE
136 CpxFlux_E(IJK) = 1.5D0 * Flux_sSE(IJK)
137 CpxFlux_N(IJK) = 1.5D0 * Flux_sSN(IJK)
138 CpxFlux_T(IJK) = 1.5D0 * Flux_sST(IJK)
139 ENDIF
140
141 IF (FLUID_AT(IJK)) THEN
142
143 IF (KT_TYPE_ENUM == GD_1999 .OR. &
144 KT_TYPE_ENUM == GTSH_2012) THEN
145 CALL SOURCE_GRANULAR_ENERGY_GD(SOURCELHS, &
146 SOURCERHS, IJK, M, IER)
147 ELSE
148 CALL SOURCE_GRANULAR_ENERGY(SOURCELHS, &
149 SOURCERHS, IJK, M, IER)
150 ENDIF
151 APO = 1.5D0*ROP_SO(IJK,M)*VOL(IJK)*ODT
152 S_P(IJK) = APO + SOURCELHS + 1.5D0* &
153 ZMAX(SUM_R_S(IJK,M)) * VOL(IJK)
154 S_C(IJK) = APO*THETA_MO(IJK,M) + SOURCERHS + 1.5d0* &
155 THETA_M(IJK,M)*ZMAX((-SUM_R_S(IJK,M))) * VOL(IJK)
156 EPS(IJK) = EP_S(IJK,M)
157
158 IF(USE_MMS) S_C(IJK) = S_C(IJK) + &
159 MMS_THETA_M_SRC(IJK)*VOL(IJK)
160 ELSE
161 EPS(IJK) = ZERO
162 S_P(IJK) = ZERO
163 S_C(IJK) = ZERO
164 IF(USE_MMS) EPS(IJK) = EP_S(IJK,M)
165 ENDIF
166 ENDDO
167
168
169 CALL CONV_DIF_PHI (THETA_M(1,M), KTH_S(1,M), DISCRETIZE(8), &
170 U_S(1,M), V_S(1,M), W_S(1,M), &
171 CpxFlux_E, CpxFlux_N, CpxFlux_T, M, A_M, B_M, IER)
172
173
174 CALL BC_PHI (THETA_M(1,M), BC_THETA_M(1,M), BC_THETAW_M(1,M), &
175 BC_HW_THETA_M(1,M), BC_C_THETA_M(1,M), M, A_M, B_M, IER)
176
177
178 CALL BC_THETA (M, A_M, B_M, IER)
179
180
181 CALL SOURCE_PHI (S_P, S_C, EPS, THETA_M(1,M), M, A_M, B_M)
182
183
184
185
186 IF (SCHAEFFER) THEN
187 DO IJK = ijkstart3, ijkend3
188 IF (FLUID_AT(IJK) .AND. &
189 EP_g(IJK) .LT. EP_g_blend_start(ijk)) THEN
190
191 A_M(IJK,1,M) = ZERO
192 A_M(IJK,-1,M) = ZERO
193 A_M(IJK,2,M) = ZERO
194 A_M(IJK,-2,M) = ZERO
195 A_M(IJK,3,M) = ZERO
196 A_M(IJK,-3,M) = ZERO
197 A_M(IJK,0,M) = -ONE
198 B_M(IJK,M) = -smallTheta
199 ENDIF
200 ENDDO
201 ENDIF
202
203 CALL CALC_RESID_S (THETA_M(1,M), A_M, B_M, M, &
204 NUM_RESID(RESID_TH,M), &
205 DEN_RESID(RESID_TH,M), RESID(RESID_TH,M), &
206 MAX_RESID(RESID_TH,M), IJK_RESID(RESID_TH,M), ZERO, IER)
207
208 CALL UNDER_RELAX_S (THETA_M(1,M), A_M, B_M, M, UR_FAC(8), IER)
209
210
211
212
213
214
215
216
217
218 CALL ADJUST_LEQ (RESID(RESID_TH,M), LEQ_IT(8), LEQ_METHOD(8),&
219 LEQI, LEQM, IER)
220
221 CALL SOLVE_LIN_EQ ('Theta_m', 8, THETA_M(1,M), A_M, B_M, M,&
222 LEQI, LEQM, LEQ_SWEEP(8), LEQ_TOL(8), LEQ_PC(8), IER)
223
224
225
226 IF(ier == -2) Err_l(myPE) = 140
227
228
229 CALL ADJUST_THETA (M, IER)
230
231 IF (IER /= 0) Err_l(myPE) = 141
232
233 ENDDO
234
235
236
237 CASE(IA_2005)
238
239 DO M = 1, SMAX
240 DO IJK = ijkstart3, ijkend3
241
242
243 IF(WALL_AT(IJK)) cycle
244
245 D_PM = D_P(IJK,M)
246 M_PM = (PI/6.d0)*(D_PM**3)*RO_S(IJK,M)
247
248
249
250 IF(.NOT.ADDED_MASS .OR. M /= M_AM) THEN
251 CpxFlux_E(IJK) = (1.5D0/M_PM) * Flux_sE(IJK,M)
252 CpxFlux_N(IJK) = (1.5D0/M_PM) * Flux_sN(IJK,M)
253 CpxFlux_T(IJK) = (1.5D0/M_PM) * Flux_sT(IJK,M)
254 ELSE
255 (IJK) = (1.5D0/M_PM) * Flux_sSE(IJK)
256 CpxFlux_N(IJK) = (1.5D0/M_PM) * Flux_sSN(IJK)
257 CpxFlux_T(IJK) = (1.5D0/M_PM) * Flux_sST(IJK)
258 ENDIF
259
260 IF (FLUID_AT(IJK)) THEN
261
262 CALL SOURCE_GRANULAR_ENERGY_IA(SOURCELHS, &
263 SOURCERHS, IJK, M, IER)
264 APO = (1.5D0/M_PM)*ROP_SO(IJK,M)*VOL(IJK)*ODT
265 S_P(IJK) = APO + SOURCELHS + (1.5d0/M_PM)*&
266 ZMAX(SUM_R_S(IJK,M)) * VOL(IJK)
267 S_C(IJK) = APO*THETA_MO(IJK,M) + SOURCERHS + &
268 (1.50d0/M_PM)*THETA_M(IJK,M)*&
269 ZMAX((-SUM_R_S(IJK,M))) * VOL(IJK)
270 EPS(IJK) = EP_S(IJK,M)
271 ELSE
272 EPS(IJK) = ZERO
273 S_P(IJK) = ZERO
274 S_C(IJK) = ZERO
275 ENDIF
276 ENDDO
277
278
279 CALL CONV_DIF_PHI (THETA_M(1,M), KTH_S(1,M), DISCRETIZE(8),&
280 U_S(1,M), V_S(1,M), W_S(1,M), &
281 CpxFlux_E, CpxFlux_N, CpxFlux_T, M, A_M, B_M, IER)
282
283
284 CALL BC_PHI (THETA_M(1,M), BC_THETA_M(1,M), BC_THETAW_M(1,M),&
285 BC_HW_THETA_M(1,M),BC_C_THETA_M(1,M), M, A_M, B_M, IER)
286
287
288 CALL BC_THETA (M, A_M, B_M, IER)
289
290
291 CALL SOURCE_PHI (S_P, S_C, EPS, THETA_M(1,M), M, A_M, B_M)
292 ENDDO
293
294
295
296 IF (SMAX > 1) THEN
297 CALL CALC_VTC_SS (VXTC_SS, IER)
298 CALL PARTIAL_ELIM_IA (THETA_M, VXTC_SS, A_M, B_M, IER)
299 ENDIF
300
301
302
303
304 IF (SCHAEFFER) THEN
305 DO M = 1, SMAX
306 DO IJK = ijkstart3, ijkend3
307 IF (FLUID_AT(IJK) .AND. &
308 EP_g(IJK) .LT. EP_g_blend_start(ijk)) THEN
309
310 D_PM = D_P(IJK,M)
311 M_PM = (PI/6.d0)*(D_PM**3)*RO_S(IJK,M)
312 A_M(IJK,1,M) = ZERO
313 A_M(IJK,-1,M) = ZERO
314 A_M(IJK,2,M) = ZERO
315 A_M(IJK,-2,M) = ZERO
316 A_M(IJK,3,M) = ZERO
317 A_M(IJK,-3,M) = ZERO
318 A_M(IJK,0,M) = -ONE
319
320
321
322
323 (IJK,M) = -smallTheta*M_PM
324 ENDIF
325 ENDDO
326 ENDDO
327 ENDIF
328
329 DO M = 1, SMAX
330 CALL CALC_RESID_S (THETA_M(1,M), A_M, B_M, M, &
331 NUM_RESID(RESID_TH,M), DEN_RESID(RESID_TH,M), RESID(RESID_TH,M),&
332 MAX_RESID(RESID_TH,M), IJK_RESID(RESID_TH,M), ZERO, IER)
333
334 CALL UNDER_RELAX_S (THETA_M(1,M), A_M, B_M, M, UR_FAC(8), IER)
335
336 CALL ADJUST_LEQ (RESID(RESID_TH,M), LEQ_IT(8), &
337 LEQ_METHOD(8), LEQI, LEQM, IER)
338
339 CALL SOLVE_LIN_EQ ('Theta_m', 8, THETA_M(1,M), A_M, B_M, M,&
340 LEQI, LEQM, LEQ_SWEEP(8), LEQ_TOL(8), LEQ_PC(8), IER)
341
342
343
344 IF(ier == -2) Err_l(myPE) = 140
345
346
347 CALL ADJUST_THETA (M, IER)
348
349 IF (IER /= 0) Err_l(myPE) = 141
350
351 ENDDO
352
353
354
355
356 CASE (GHD_2007)
357
358
359 = MMAX
360
361
362 (:) = ZERO
363 TOT_SUM_RS(:) = ZERO
364 TOT_EPS(:) = ZERO
365
366 CALL CALC_NFLUX (IER)
367
368 DO IJK = ijkstart3, ijkend3
369
370
371 (IJK) = 1.5D0 * Flux_nE(IJK)
372 CpxFlux_N(IJK) = 1.5D0 * Flux_nN(IJK)
373 CpxFlux_T(IJK) = 1.5D0 * Flux_nT(IJK)
374
375 IF (FLUID_AT(IJK)) THEN
376 DO L = 1,SMAX
377 M_PM = (PI/6.d0)*(D_P(IJK,L)**3)*RO_S(IJK,L)
378 TOT_SUM_RS(IJK) = TOT_SUM_RS(IJK) + SUM_R_S(IJK,L)/M_PM
379 TOT_NO(IJK) = TOT_NO(IJK) + ROP_SO(IJK,L)/M_PM
380 TOT_EPS(IJK) = TOT_EPS(IJK) + EP_S(IJK,L)
381 ENDDO
382
383
384 CALL SOURCE_GHD_GRANULAR_ENERGY (SOURCELHS, &
385 SOURCERHS, IJK, IER)
386 APO = 1.5D0 * TOT_NO(IJK)*VOL(IJK)*ODT
387 S_P(IJK) = APO + SOURCELHS + 1.5d0 *&
388 ZMAX(TOT_SUM_RS(IJK)) * VOL(IJK)
389 S_C(IJK) = APO*THETA_MO(IJK,M) + SOURCERHS + 1.50d0 *&
390 THETA_M(IJK,M)*ZMAX(-TOT_SUM_RS(IJK)) * VOL(IJK)
391 EPS(IJK) = TOT_EPS(IJK)
392 ELSE
393 EPS(IJK) = ZERO
394 S_P(IJK) = ZERO
395 S_C(IJK) = ZERO
396 ENDIF
397 ENDDO
398
399
400 CALL CONV_DIF_PHI (THETA_M(1,M), KTH_S(1,M), DISCRETIZE(8),&
401 U_S(1,M), V_S(1,M), W_S(1,M), &
402 CpxFlux_E, CpxFlux_N, CpxFlux_T, M, A_M, B_M, IER)
403
404
405 CALL BC_PHI (THETA_M(1,M), BC_THETA_M(1,M), BC_THETAW_M(1,M), &
406 BC_HW_THETA_M(1,M), BC_C_THETA_M(1,M), M, A_M, B_M, IER)
407
408
409 CALL BC_THETA (M, A_M, B_M, IER)
410
411
412 CALL SOURCE_PHI (S_P, S_C, EPS, THETA_M(1,M), M, A_M, B_M)
413
414 CALL CALC_RESID_S (THETA_M(1,M), A_M, B_M, M, &
415 NUM_RESID(RESID_TH,M), DEN_RESID(RESID_TH,M), RESID(RESID_TH,M),&
416 MAX_RESID(RESID_TH,M), IJK_RESID(RESID_TH,M), ZERO, IER)
417
418 CALL UNDER_RELAX_S (THETA_M(1,M), A_M, B_M, M, UR_FAC(8), IER)
419
420 CALL ADJUST_LEQ (RESID(RESID_TH,M), LEQ_IT(8), LEQ_METHOD(8),&
421 LEQI, LEQM, IER)
422
423 CALL SOLVE_LIN_EQ ('Theta_m', 8, THETA_M(1,M), A_M, B_M, M, &
424 LEQI, LEQM, LEQ_SWEEP(8), LEQ_TOL(8), LEQ_PC(8), IER)
425
426
427
428 IF(ier == -2) Err_l(myPE) = 140
429
430
431 CALL ADJUST_THETA (M, IER)
432
433
434 IF (IER /= 0) Err_l(myPE) = 141
435
436
437
438
439 CASE DEFAULT
440
441 WRITE (*, '(A)') 'ADJUST_THETA'
442 WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
443 call mfix_exit(myPE)
444 END SELECT
445
446 call unlock_ambm
447 call unlock_tmp_array
448
449
450 CALL global_all_sum(Err_l, Err_g)
451 IER = maxval(Err_g)
452
453
454 RETURN
455 END SUBROUTINE SOLVE_GRANULAR_ENERGY
456
457
458