File: /nfs/home/0/users/jenkins/mfix.git/model/check_mass_balance.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 SUBROUTINE CHECK_Mass_balance (init)
19
20
21
22
23
24
25 USE param
26 USE param1
27 USE toleranc
28 USE fldvar
29 USE rxns
30 USE geometry
31 USE run
32 USE bc
33 USE constant
34 USE physprop
35 USE indices
36 USE funits
37 USE compar
38 USE mpi_utility
39 USE output
40 USE check
41 USE mflux
42 USE functions
43 IMPLICIT NONE
44
45
46
47
48
49
50
51
52
53
54
55
56 INTEGER init
57
58 INTEGER IER
59
60
61 INTEGER M
62
63
64 INTEGER N
65
66
67 INTEGER L
68
69
70 INTEGER I1, I2, J1, J2, K1, K2
71
72
73 DOUBLE PRECISION Accumulation_old, Accumulation_delta, flux, flux_in_tot, flux_out_tot, error_percent
74 DOUBLE PRECISION fin, fout
75
76
77 DOUBLE PRECISION Accumulation, Accumulation_sp
78
79
80 if(report_mass_balance_dt == UNDEFINED) return
81
82 if(init == 0) then
83
84
85 = time
86 report_time = time + report_mass_balance_dt
87
88
89 DO L = 1, DIMENSION_BC
90 flux_out_g(L) = ZERO
91 flux_in_g(L) = ZERO
92 DO N = 1, NMAX(0)
93 flux_out_X_g(L, N) = ZERO
94 flux_in_X_g(L, N) = ZERO
95 END DO
96
97 DO M = 1, MMAX
98 flux_out_s(L, M) = ZERO
99 flux_in_s(L, M) = ZERO
100 DO N = 1, NMAX(M)
101 flux_out_X_s(L, M, N) = ZERO
102 flux_in_X_s(L, M, N) = ZERO
103 END DO
104 END DO
105
106 END DO
107
108 Integral_SUM_R_g = ZERO
109 DO N = 1, NMAX(0)
110 Integral_R_g(N) = ZERO
111 END DO
112
113 DO M = 1, MMAX
114 Integral_SUM_R_s(M) = ZERO
115 DO N = 1, NMAX(M)
116 Integral_R_s(M, N) = ZERO
117 END DO
118 END DO
119
120
121 = Accumulation(ROP_g)
122 DO N = 1, NMAX(0)
123 Accumulation_X_g(N) = Accumulation_sp(ROP_g, X_g(1, N))
124 END DO
125
126 DO M = 1, MMAX
127 Accumulation_s(M) = Accumulation(ROP_s(1,M))
128 DO N = 1, NMAX(M)
129 Accumulation_X_s(M, N) = Accumulation_sp(ROP_s(1, M), X_s(1, M, N))
130 END DO
131 END DO
132
133 return
134
135
136 else
137
138
139 DO L = 1, DIMENSION_BC
140 IF (BC_DEFINED(L)) THEN
141 I1 = BC_I_W(L)
142 I2 = BC_I_E(L)
143 J1 = BC_J_S(L)
144 J2 = BC_J_N(L)
145 K1 = BC_K_B(L)
146 K2 = BC_K_T(L)
147
148
149 IF(.NOT.Added_Mass) THEN
150 call Calc_mass_fluxHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_gE, Flux_gN, Flux_gT, fin, fout, IER)
151 ELSE
152 call Calc_mass_fluxHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_gSE, Flux_gSN, Flux_gST, fin, fout, IER)
153 ENDIF
154 flux_out_g(L) = flux_out_g(L) + fout * dt_prev
155 flux_in_g(L) = flux_in_g(L) + fin * dt_prev
156
157 DO M = 1, MMAX
158
159 IF(.NOT.Added_Mass .OR. M /= M_AM) THEN
160 call Calc_mass_fluxHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_sE(1,m), Flux_sN(1,m), Flux_sT(1,m), fin, fout, IER)
161 ELSE
162 call Calc_mass_fluxHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_sSE, Flux_sSN, Flux_sST, fin, fout, IER)
163 ENDIF
164 flux_out_s(L, M) = flux_out_s(L, M) + fout * dt_prev
165 flux_in_s(L, M) = flux_in_s(L, M) + fin * dt_prev
166 END DO
167
168 ENDIF
169 END DO
170
171 Integral_SUM_R_g = Integral_SUM_R_g + Accumulation(SUM_R_g) * dt_prev
172 IF(SPECIES_EQ(0))THEN
173 DO N = 1, NMAX(0)
174 DO L = 1, DIMENSION_BC
175 IF (BC_DEFINED(L)) THEN
176 I1 = BC_I_W(L)
177 I2 = BC_I_E(L)
178 J1 = BC_J_S(L)
179 J2 = BC_J_N(L)
180 K1 = BC_K_B(L)
181 K2 = BC_K_T(L)
182
183 IF(.NOT.Added_Mass) THEN
184 call Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_gE, &
185 Flux_gN, Flux_gT, X_g(1, N), DISCRETIZE(7), fin, fout, IER)
186 ELSE
187 call Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), Flux_gSE, &
188 Flux_gSN, Flux_gST, X_g(1, N), DISCRETIZE(7), fin, fout, IER)
189 ENDIF
190 flux_out_X_g(L, N) = flux_out_X_g(L, N) + fout * dt_prev
191 flux_in_X_g(L, N) = flux_in_X_g(L, N) + fin * dt_prev
192 ENDIF
193 END DO
194
195 Integral_R_g(N) = Integral_R_g(N) + (Accumulation(R_gp(1,N)) - &
196 Accumulation_sp(ROX_gc(1,N), X_g(1,N)) )* dt_prev
197 END DO
198 ENDIF
199
200 DO M = 1, MMAX
201 IF(SPECIES_EQ(M))THEN
202 Integral_SUM_R_s(M) = Integral_SUM_R_s(M) + Accumulation(SUM_R_s(1,M)) * dt_prev
203 DO N = 1, NMAX(M)
204 DO L = 1, DIMENSION_BC
205 IF (BC_DEFINED(L)) THEN
206 I1 = BC_I_W(L)
207 I2 = BC_I_E(L)
208 J1 = BC_J_S(L)
209 J2 = BC_J_N(L)
210 K1 = BC_K_B(L)
211 K2 = BC_K_T(L)
212
213
214 IF(.NOT.Added_Mass .OR. M /= M_AM) THEN
215 call Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), &
216 Flux_sE(1,M), Flux_sN(1,M), Flux_sT(1,M), X_s(1, M, N), DISCRETIZE(7), fin, fout, &
217 IER)
218 ELSE
219 call Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, BC_PLANE(L), &
220 Flux_sSE, Flux_sSN, Flux_sST, X_s(1, M, N), DISCRETIZE(7), fin, fout, IER)
221 ENDIF
222 flux_out_X_s(L, M, N) = flux_out_X_s(L, M, N) + fout * dt_prev
223 flux_in_X_s(L, M, N) = flux_in_X_s(L, M, N) + fin * dt_prev
224 ENDIF
225 END DO
226
227 Integral_R_s(M, N) = Integral_R_s(M, N) + (Accumulation(R_sp(1,M,N)) - &
228 Accumulation_sp(ROX_sc(1,M,N), X_s(1,M,N)) )* dt_prev
229 END DO
230 ENDIF
231 END DO
232
233 endif
234
235
236 if (time >= report_time)then
237
238 CALL START_LOG
239 IF(DMP_LOG)WRITE(UNIT_LOG, '(/A,G12.5,A,G12.5)') 'Mass balance for interval ', start_time, ' to ', time
240
241 Accumulation_old = Accumulation_g
242 Accumulation_g = Accumulation(ROP_g)
243 Accumulation_delta = Accumulation_g - Accumulation_old - Integral_SUM_R_g
244 IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Total Fluid Accumulation (g)'
245 IF(DMP_LOG)WRITE(UNIT_LOG, '(4(A,G12.5))') ' Old = ', Accumulation_old, ', New = ', &
246 Accumulation_g, ', Production = ', Integral_SUM_R_g, ', net accu(New - Old - Production) = ', Accumulation_delta
247
248 IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Integral of boundary flux (g)'
249 IF(DMP_LOG)WRITE(Unit_log, '(A, T8, A, T21, A, T34, A)')' BC#', 'in', 'out', '(in - out)'
250 flux = zero
251 flux_in_tot = zero
252 flux_out_tot = zero
253 DO L = 1, DIMENSION_BC
254 if(flux_out_g(L) /= ZERO .OR. flux_in_g(L) /= ZERO) then
255 IF(DMP_LOG)WRITE(Unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')L, flux_in_g(L), flux_out_g(L), (flux_in_g(L)-flux_out_g(L))
256 endif
257 flux = flux + flux_in_g(L) - flux_out_g(L)
258 flux_in_tot = flux_in_tot + flux_in_g(L)
259 flux_out_tot = flux_out_tot + flux_out_g(L)
260 END DO
261 if((flux - Accumulation_delta) /= zero) then
262 error_percent = undefined
263 if(flux_in_tot /= zero) error_percent = (flux - Accumulation_delta)*100./flux_in_tot
264 else
265 error_percent = zero
266 endif
267 IF(DMP_LOG)WRITE(Unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
268 IF(DMP_LOG)WRITE(Unit_log, '(A, G12.5, A, G12.5)')'Error (net influx - net accu) = ', &
269 (flux - Accumulation_delta), ' %Error_g = ', error_percent
270
271 DO M = 1, MMAX
272 Accumulation_old = Accumulation_s(M)
273 Accumulation_s(M) = Accumulation(ROP_s(1, M))
274 Accumulation_delta = Accumulation_s(M) - Accumulation_old - Integral_SUM_R_s(M)
275 IF(DMP_LOG)WRITE(UNIT_LOG, '(/A, I1, A)') 'Total Solids-', M, ' Accumulation (g)'
276 IF(DMP_LOG)WRITE(UNIT_LOG, '(4(A,G12.5))') ' Old = ', Accumulation_old, ', New = ', &
277 Accumulation_s(M), ', Production = ', Integral_SUM_R_s(M), ', net accu(New - Old - Production) = ', Accumulation_delta
278
279 IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Integral of boundary flux (g)'
280 IF(DMP_LOG)WRITE(Unit_log, '(A, T8, A, T21, A, T34, A)')' BC#', 'in', 'out', '(in - out)'
281 flux = zero
282 flux_in_tot = zero
283 flux_out_tot = zero
284 DO L = 1, DIMENSION_BC
285 if(flux_out_s(L,M) /= ZERO .OR. flux_in_s(L,M) /= ZERO) then
286 IF(DMP_LOG) THEN
287 WRITE(Unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')L, &
288 flux_in_s(L,M), flux_out_s(L,M), (flux_in_s(L,M)-flux_out_s(L,M))
289 ENDIF
290 endif
291 flux = flux + flux_in_s(L,M) - flux_out_s(L,M)
292 flux_in_tot = flux_in_tot + flux_in_s(L,M)
293 flux_out_tot = flux_out_tot + flux_out_s(L,M)
294 END DO
295 if((flux - Accumulation_delta) /= zero) then
296 if(flux_in_tot /= zero) then
297 error_percent = (flux - Accumulation_delta)*100./flux_in_tot
298 else
299 error_percent = (flux - Accumulation_delta)*100./Accumulation_old
300 endif
301 else
302 error_percent = zero
303 endif
304 IF(DMP_LOG)WRITE(Unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
305 IF(DMP_LOG)WRITE(Unit_log, '(A, G12.5, A, I1, A, G12.5)')'Error (net influx - net accu) = ', &
306 (flux - Accumulation_delta), ' %Error_',M,' = ', error_percent
307 END DO
308
309
310 IF(SPECIES_EQ(0) )THEN
311 DO N = 1, NMAX(0)
312
313 IF(DMP_LOG)WRITE(UNIT_LOG, '(/A,I2)') 'Gas species - ', N
314
315 Accumulation_old = Accumulation_X_g(N)
316 Accumulation_X_g(N) = Accumulation_sp(ROP_g, X_g(1, N))
317 Accumulation_delta = Accumulation_X_g(N) - Accumulation_old - Integral_R_g(N)
318 IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Species Accumulation (g)'
319 IF(DMP_LOG)WRITE(UNIT_LOG, '(4(A,G12.5))') ' Old = ', Accumulation_old, ', New = ', &
320 Accumulation_X_g(N), ', Production = ', Integral_R_g(N), ', net accu(New - Old - Production) = ', Accumulation_delta
321
322 IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Integral of boundary flux (g)'
323 IF(DMP_LOG)WRITE(Unit_log, '(A, T8, A, T21, A, T34, A)')' BC#', 'in', 'out', '(in - out)'
324 flux = zero
325 flux_in_tot = zero
326 flux_out_tot = zero
327 DO L = 1, DIMENSION_BC
328 if(flux_out_X_g(L, N) /= ZERO .OR. flux_in_X_g(L, N) /= ZERO) then
329 IF(DMP_LOG)WRITE(Unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')L, flux_in_X_g(L, N), flux_out_X_g(L, N), &
330 (flux_in_X_g(L, N)-flux_out_X_g(L, N))
331 endif
332 flux = flux + flux_in_X_g(L, N) - flux_out_X_g(L, N)
333 flux_in_tot = flux_in_tot + flux_in_X_g(L, N)
334 flux_out_tot = flux_out_tot + flux_out_X_g(L, N)
335 END DO
336 if((flux - Accumulation_delta) /= zero) then
337 error_percent = undefined
338 if(flux_in_tot /= zero) then
339 error_percent = (flux - Accumulation_delta)*100./flux_in_tot
340 elseif (Integral_R_g(N) /= zero) then
341 error_percent = (flux - Accumulation_delta)*100./Integral_R_g(N)
342 endif
343 else
344 error_percent = zero
345 endif
346 IF(DMP_LOG)WRITE(Unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
347 IF(DMP_LOG)WRITE(Unit_log, '(A, G12.5, A, I2, A, G12.5)')'Error (net influx - net accu) = ', &
348 (flux - Accumulation_delta), ' %Error_g(',N,') = ', error_percent
349
350 END DO
351 ENDIF
352
353 DO M = 1, MMAX
354 IF(SPECIES_EQ(M) )THEN
355 DO N = 1, NMAX(M)
356
357 IF(DMP_LOG)WRITE(UNIT_LOG, '(/A,I1, A, I2)') 'Solids-', M, ' species - ', N
358
359 Accumulation_old = Accumulation_X_s(M,N)
360 Accumulation_X_s(M,N) = Accumulation_sp(ROP_s(1,M), X_s(1, M, N))
361 Accumulation_delta = Accumulation_X_s(M,N) - Accumulation_old - Integral_R_s(M,N)
362 IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Species Accumulation (g)'
363 IF(DMP_LOG)WRITE(UNIT_LOG, '(4(A,G12.5))') ' Old = ', Accumulation_old, ', New = ', &
364 Accumulation_X_s(M,N), ', Production = ', Integral_R_s(M,N), &
365 ', net accu(New - Old - Production) = ', Accumulation_delta
366
367 IF(DMP_LOG)WRITE(UNIT_LOG, '(A)') 'Integral of boundary flux (g)'
368 IF(DMP_LOG)WRITE(Unit_log, '(A, T8, A, T21, A, T34, A)')' BC#', 'in', 'out', '(in - out)'
369 flux = zero
370 flux_in_tot = zero
371 flux_out_tot = zero
372 DO L = 1, DIMENSION_BC
373 if(flux_out_X_s(L, M, N) /= ZERO .OR. flux_in_X_s(L, M, N) /= ZERO) then
374 IF(DMP_LOG)WRITE(Unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')L, flux_in_X_s(L, M, N), flux_out_X_s(L, M, N), &
375 (flux_in_X_s(L, M, N)-flux_out_X_s(L, M, N))
376 endif
377 flux = flux + flux_in_X_s(L, M, N) - flux_out_X_s(L, M, N)
378 flux_in_tot = flux_in_tot + flux_in_X_s(L, M, N)
379 flux_out_tot = flux_out_tot + flux_out_X_s(L, M, N)
380 END DO
381 if((flux - Accumulation_delta) /= zero) then
382 error_percent = undefined
383 if(flux_in_tot /= zero) then
384 error_percent = (flux - Accumulation_delta)*100./flux_in_tot
385 elseif (Integral_R_s(M,N) /= zero)then
386 error_percent = (flux - Accumulation_delta)*100./Integral_R_s(M,N)
387 endif
388 else
389 error_percent = zero
390 endif
391 IF(DMP_LOG)WRITE(Unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
392 IF(DMP_LOG)WRITE(Unit_log, '(A, G12.5, A, I1, A, I2, A, G12.5)')'Error (net influx - net accu) = ', &
393 (flux - Accumulation_delta), ' %Error_',M,'(',N,') = ', error_percent
394
395 END DO
396 ENDIF
397 END DO
398 IF(DMP_LOG)WRITE(UNIT_LOG, '(/)')
399 CALL END_LOG
400
401 start_time = time
402 report_time = time + report_mass_balance_dt
403
404
405 DO L = 1, DIMENSION_BC
406 flux_out_g(L) = ZERO
407 flux_in_g(L) = ZERO
408 DO N = 1, NMAX(0)
409 flux_out_X_g(L, N) = ZERO
410 flux_in_X_g(L, N) = ZERO
411 END DO
412
413 DO M = 1, MMAX
414 flux_out_s(L, M) = ZERO
415 flux_in_s(L, M) = ZERO
416 DO N = 1, NMAX(M)
417 flux_out_X_s(L, M, N) = ZERO
418 flux_in_X_s(L, M, N) = ZERO
419 END DO
420 END DO
421
422 END DO
423
424 Integral_SUM_R_g = ZERO
425 DO N = 1, NMAX(0)
426 Integral_R_g(N) = ZERO
427 END DO
428
429 DO M = 1, MMAX
430 Integral_SUM_R_s(M) = ZERO
431 DO N = 1, NMAX(M)
432 Integral_R_s(M, N) = ZERO
433 END DO
434 END DO
435
436
437 endif
438
439
440 RETURN
441 END SUBROUTINE CHECK_Mass_balance
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462 SUBROUTINE Calc_mass_flux(I1, I2, J1, J2, K1, K2, Plane, U, V, W, ROP, flux_in, flux_out, IER)
463 USE param
464 USE param1
465 IMPLICIT NONE
466
467
468
469 INTEGER :: I1, I2, J1, J2, K1, K2
470
471
472 Character :: Plane
473
474
475 DOUBLE PRECISION :: U(DIMENSION_3), V(DIMENSION_3), W(DIMENSION_3)
476
477
478 DOUBLE PRECISION :: ROP(DIMENSION_3)
479
480
481 DOUBLE PRECISION :: Xn(DIMENSION_3)
482
483
484
485
486 DOUBLE PRECISION :: flux_in, flux_out
487
488 INTEGER :: IER
489
490 Xn = one
491 call Calc_mass_flux_sp(I1, I2, J1, J2, K1, K2, PLANE, U, V, W, ROP, Xn, flux_in, flux_out, IER)
492 return
493 end SUBROUTINE Calc_mass_flux
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512 SUBROUTINE Calc_mass_flux_sp(I1, I2, J1, J2, K1, K2, Plane, U, V, W, ROP, Xn, flux_in, flux_out, IER)
513 USE param
514 USE param1
515 USE geometry
516 USE physprop
517 USE indices
518 USE compar
519 USE mpi_utility
520 USE functions
521 IMPLICIT NONE
522
523
524
525 INTEGER :: I1, I2, J1, J2, K1, K2
526
527
528 Character :: Plane
529
530
531 DOUBLE PRECISION :: U(DIMENSION_3), V(DIMENSION_3), W(DIMENSION_3)
532
533
534 DOUBLE PRECISION :: ROP(DIMENSION_3)
535
536
537 DOUBLE PRECISION :: Xn(DIMENSION_3)
538
539
540
541
542 DOUBLE PRECISION :: flux_in, flux_out
543
544 INTEGER :: IER
545
546
547 INTEGER I, J, K, IJK
548
549 flux_in = zero
550 flux_out = zero
551 IER = 0
552
553
554 DO K = K1, K2
555 DO J = J1, J2
556 DO I = I1, I2
557 IF(.NOT.IS_ON_myPE_owns(I, J, K)) cycle
558 IF(DEAD_CELL_AT(I,J,K)) cycle
559
560 SELECT CASE (TRIM(PLANE))
561 CASE ('W')
562 IJK = IM_OF(FUNIJK(I,J,K))
563 IF(U(IJK) > ZERO)THEN
564 flux_out = flux_out + AYZ(IJK) * U(IJK) * ROP(IJK) * Xn(IJK)
565 ELSE
566 flux_in = flux_in - AYZ(IJK) * U(IJK) * ROP(IP_OF(IJK)) * Xn(IP_OF(IJK))
567 ENDIF
568
569 CASE ('E')
570 IJK = FUNIJK(I,J,K)
571 IF(U(IJK) > ZERO)THEN
572 flux_in = flux_in + AYZ(IJK) * U(IJK) * ROP(IJK) * Xn(IJK)
573 ELSE
574 flux_out = flux_out - AYZ(IJK) * U(IJK) * ROP(IP_OF(IJK)) * Xn(IP_OF(IJK))
575 ENDIF
576
577 CASE ('S')
578 IJK = JM_OF(FUNIJK(I,J,K))
579 IF(V(IJK) > ZERO)THEN
580 flux_out = flux_out + AXZ(IJK) * V(IJK) * ROP(IJK) * Xn(IJK)
581 ELSE
582 flux_in = flux_in - AXZ(IJK) * V(IJK) * ROP(JP_OF(IJK)) * Xn(JP_OF(IJK))
583 ENDIF
584
585 CASE ('N')
586 IJK = FUNIJK(I,J,K)
587 IF(V(IJK) > ZERO)THEN
588 flux_in = flux_in + AXZ(IJK) * V(IJK) * ROP(IJK) * Xn(IJK)
589 ELSE
590 flux_out = flux_out - AXZ(IJK) * V(IJK) * ROP(JP_OF(IJK)) * Xn(JP_OF(IJK))
591 ENDIF
592
593 CASE ('B')
594 IJK = KM_OF(FUNIJK(I,J,K))
595 IF(W(IJK) > ZERO)THEN
596 flux_out = flux_out + AXY(IJK) * W(IJK) * ROP(IJK) * Xn(IJK)
597 ELSE
598 flux_in = flux_in - AXY(IJK) * W(IJK) * ROP(KP_OF(IJK)) * Xn(KP_OF(IJK))
599 ENDIF
600
601 CASE ('T')
602 IJK = FUNIJK(I,J,K)
603 IF(W(IJK) > ZERO)THEN
604 flux_in = flux_in + AXY(IJK) * W(IJK) * ROP(IJK) * Xn(IJK)
605 ELSE
606 flux_out = flux_out - AXY(IJK) * W(IJK) * ROP(KP_OF(IJK)) * Xn(KP_OF(IJK))
607 ENDIF
608
609 CASE DEFAULT
610
611
612
613
614
615 END SELECT
616 ENDDO
617 ENDDO
618 ENDDO
619
620 call global_all_sum(flux_in)
621 call global_all_sum(flux_out)
622
623 return
624 end SUBROUTINE Calc_mass_flux_sp
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642 SUBROUTINE Calc_mass_fluxHR(I1, I2, J1, J2, K1, K2, Plane, Flux_E, Flux_N, Flux_T, flux_in, flux_out, IER)
643 USE param
644 USE param1
645 IMPLICIT NONE
646
647
648
649 INTEGER :: I1, I2, J1, J2, K1, K2
650
651
652 Character :: Plane
653
654
655 DOUBLE PRECISION :: Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
656
657
658 DOUBLE PRECISION :: Xn(DIMENSION_3)
659
660
661
662
663 DOUBLE PRECISION :: flux_in, flux_out
664
665 INTEGER :: IER
666
667 Xn = one
668 call Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, PLANE, Flux_E, Flux_N, Flux_T, Xn, 0, flux_in, flux_out, IER)
669 return
670 end SUBROUTINE Calc_mass_fluxHR
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689 SUBROUTINE Calc_mass_flux_spHR(I1, I2, J1, J2, K1, K2, Plane, Flux_E, Flux_N, Flux_T, Xn, disc, flux_in, flux_out, IER)
690 USE param
691 USE param1
692 USE geometry
693 USE physprop
694 USE indices
695 USE compar
696 Use xsi_array
697 USE mpi_utility
698 USE xsi
699 USE functions
700 IMPLICIT NONE
701
702
703
704 INTEGER :: I1, I2, J1, J2, K1, K2
705
706
707 Character :: Plane
708
709
710 DOUBLE PRECISION :: Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
711
712
713 DOUBLE PRECISION :: Xn(DIMENSION_3)
714
715
716
717 INTEGER Disc
718
719
720
721
722 DOUBLE PRECISION :: flux_in, flux_out
723
724
725 DOUBLE PRECISION :: Flux, PHI_HO
726
727 INTEGER :: IER
728
729
730 INTEGER I, J, K, IJK
731
732 call lock_xsi_array
733
734
735
736
737 CALL CALC_XSI (DISC, Xn, Flux_E, Flux_N, Flux_T, XSI_E, XSI_N, XSI_T, 0)
738
739 flux_in = zero
740 flux_out = zero
741 IER = 0
742
743 DO K = K1, K2
744 DO J = J1, J2
745 DO I = I1, I2
746 IF(.NOT.IS_ON_myPE_owns(I, J, K)) cycle
747 IF(DEAD_CELL_AT(I,J,K)) cycle
748
749 SELECT CASE (TRIM(PLANE))
750 CASE ('W')
751 IJK = IM_OF(FUNIJK(I,J,K))
752 PHI_HO = XSI_E(IJK)*Xn(IP_OF(IJK))+(1.0-XSI_E(IJK))*Xn(IJK)
753 Flux = FLUX_E(IJK)
754
755 CASE ('E')
756 IJK = FUNIJK(I,J,K)
757 PHI_HO = XSI_E(IJK)*Xn(IP_OF(IJK))+(1.0-XSI_E(IJK))*Xn(IJK)
758 Flux = -FLUX_E(IJK)
759
760 CASE ('S')
761 IJK = JM_OF(FUNIJK(I,J,K))
762 PHI_HO = XSI_N(IJK)*Xn(JP_OF(IJK))+(1.0-XSI_N(IJK))*Xn(IJK)
763 Flux = FLUX_N(IJK)
764
765 CASE ('N')
766 IJK = FUNIJK(I,J,K)
767 PHI_HO = XSI_N(IJK)*Xn(JP_OF(IJK))+(1.0-XSI_N(IJK))*Xn(IJK)
768 Flux = -FLUX_N(IJK)
769
770 CASE ('B')
771 IJK = KM_OF(FUNIJK(I,J,K))
772 PHI_HO = XSI_T(IJK)*Xn(KP_OF(IJK))+(1.0-XSI_T(IJK))*Xn(IJK)
773 Flux = FLUX_T(IJK)
774
775 CASE ('T')
776 IJK = FUNIJK(I,J,K)
777 PHI_HO = XSI_T(IJK)*Xn(KP_OF(IJK))+(1.0-XSI_T(IJK))*Xn(IJK)
778 Flux = -FLUX_T(IJK)
779
780 CASE DEFAULT
781 PHI_HO = ZERO
782 Flux = ZERO
783
784 END SELECT
785
786 IF(FLUX > ZERO)THEN
787 flux_out = flux_out + FLUX * PHI_HO
788 ELSE
789 flux_in = flux_in - FLUX * PHI_HO
790 ENDIF
791
792 ENDDO
793 ENDDO
794 ENDDO
795
796 call global_all_sum(flux_in)
797 call global_all_sum(flux_out)
798 call unlock_xsi_array
799
800 return
801 end SUBROUTINE Calc_mass_flux_spHR
802
803
804
805
806
807
808
809
810
811
812
813 DOUBLE PRECISION FUNCTION Accumulation(ro)
814
815
816
817
818
819
820
821
822 USE param
823 USE param1
824 IMPLICIT NONE
825
826
827
828
829
830
831
832 DOUBLE PRECISION :: RO(DIMENSION_3)
833
834 DOUBLE PRECISION :: Xn(DIMENSION_3)
835
836
837 DOUBLE PRECISION Accumulation_sp
838
839
840
841 = ONE
842 Accumulation = Accumulation_sp(ro, xn)
843
844 RETURN
845 END FUNCTION Accumulation
846
847
848
849
850
851
852
853
854
855
856
857 DOUBLE PRECISION FUNCTION Accumulation_sp(ro, Xn)
858
859
860
861
862
863
864
865
866 USE param
867 USE param1
868 USE parallel
869 USE physprop
870 USE geometry
871 USE indices
872 USE compar
873 USE mpi_utility
874 USE functions
875 IMPLICIT NONE
876
877
878
879
880
881
882
883 INTEGER IJK
884
885
886 DOUBLE PRECISION :: RO(DIMENSION_3)
887
888 DOUBLE PRECISION :: Xn(DIMENSION_3)
889
890 DOUBLE PRECISION SUM
891
892
893
894 = ZERO
895
896
897
898 DO IJK = ijkstart3, ijkend3
899 IF(.NOT.IS_ON_myPE_wobnd(I_OF(IJK), J_OF(IJK), K_OF(IJK))) cycle
900 IF (FLUID_AT(IJK)) SUM = SUM + RO(IJK) * Xn(IJK) * VOL(IJK)
901 END DO
902
903 call global_all_sum(sum)
904 Accumulation_sp = sum
905
906 RETURN
907 END FUNCTION Accumulation_sp
908
909
910
911 DOUBLE PRECISION FUNCTION check_conservation (Phi, A_m, B_m, IJK)
912 USE param
913 USE param1
914 USE geometry
915 USE indices
916 USE matrix
917 USE compar
918 USE functions
919 IMPLICIT NONE
920
921
922 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
923
924
925 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
926
927 DOUBLE PRECISION Phi(DIMENSION_3)
928
929
930 INTEGER IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP
931
932 IMJK = IM_OF(IJK)
933 IJMK = JM_OF(IJK)
934 IJKM = KM_OF(IJK)
935 IPJK = IP_OF(IJK)
936 IJPK = JP_OF(IJK)
937 IJKP = KP_OF(IJK)
938
939 check_conservation = Phi(IJK) * A_m(IJK, 0, 0) &
940 + Phi(IPJK) * A_m(IJK, E, 0) &
941 + Phi(IMJK) * A_m(IJK, W, 0) &
942 + Phi(IJPK) * A_m(IJK, N, 0) &
943 + Phi(IJMK) * A_m(IJK, S, 0) &
944 + Phi(IJKP) * A_m(IJK, T, 0) &
945 + Phi(IJKM) * A_m(IJK, B, 0) &
946 - B_m(IJK,0)
947
948 END FUNCTION check_conservation
949
950
951
952