MFIX  2016-1
check_mod.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Module name: CHECK(init) C
4 ! Purpose: Check global species and elemental balances C
5 ! C
6 ! Author: M. Syamlal Date: 27-Nov-02 C
7 ! Reviewer: Date: C
8 ! C
9 ! Literature/Document References: C
10 ! C
11 ! Variables referenced: C
12 ! Variables modified: C
13 ! C
14 ! Local variables: ABORT, SUM C
15 ! C
16 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
17 !
18  MODULE check
19 
20  Use param, only: dimension_bc, dim_m, dim_n_g, dim_n_s
21  ! Use param1
22 ! variables for check_mass_balance
23  DOUBLE PRECISION start_time, report_time
24 
29 
33  DOUBLE PRECISION flux_out_x_s(dimension_bc,dim_m, dim_n_s), flux_in_x_s(dimension_bc,dim_m, dim_n_s)
34 
35  CONTAINS
36 
37  SUBROUTINE check_mass_balance (init)
38 !
39 ! Include param.inc file to specify parameter values
40 !
41 !-----------------------------------------------
42 ! M o d u l e s
43 !-----------------------------------------------
44  USE bc
45  USE compar
46  USE constant
47  USE fldvar
48  USE functions
49  USE funits
50  USE geometry
51  USE indices
52  USE machine, only: start_log, end_log
53  USE mflux
54  USE mpi_utility
55  USE output
56  USE param
57  USE param1
58  USE physprop
60  USE rxns
61  USE toleranc
62  IMPLICIT NONE
63 !-----------------------------------------------
64 ! G l o b a l P a r a m e t e r s
65 !-----------------------------------------------
66 !-----------------------------------------------
67 ! L o c a l P a r a m e t e r s
68 !-----------------------------------------------
69 !-----------------------------------------------
70 ! L o c a l V a r i a b l e s
71 !-----------------------------------------------
72 !
73 ! Flag to check whether this is an initialization call
74 ! 0 -> initialization call; 1 -> integration call
75  INTEGER init
76 !
77  INTEGER IER
78 !
79 ! Solids phase
80  INTEGER M
81 !
82 ! Species index
83  INTEGER NN
84 !
85 ! Do-loop counter
86  INTEGER L
87 !
88 ! Starting and ending indices (Make sure these represent a plane and NOT a volume)
89  INTEGER I1, I2, J1, J2, K1, K2
90 !
91 !
92  DOUBLE PRECISION Accumulation_old, Accumulation_delta, flux, flux_in_tot, flux_out_tot, error_percent
93  DOUBLE PRECISION fin, fout
94 
95 !-----------------------------------------------
96 
97  if(report_mass_balance_dt == undefined) return
98 
99  if(init == 0) then
100 ! allocate arrays
101 ! Initilaize this routine
102  start_time = time
104 
105  !initialize flux and reaction rate arrays
106  DO l = 1, dimension_bc
107  flux_out_g(l) = zero
108  flux_in_g(l) = zero
109  DO nn = 1, nmax(0)
110  flux_out_x_g(l, nn) = zero
111  flux_in_x_g(l, nn) = zero
112  END DO
113 
114  DO m = 1, mmax
115  flux_out_s(l, m) = zero
116  flux_in_s(l, m) = zero
117  DO nn = 1, nmax(m)
118  flux_out_x_s(l, m, nn) = zero
119  flux_in_x_s(l, m, nn) = zero
120  END DO
121  END DO
122 
123  END DO
124 
126  DO nn = 1, nmax(0)
127  integral_r_g(nn) = zero
128  END DO
129 
130  DO m = 1, mmax
132  DO nn = 1, nmax(m)
133  integral_r_s(m, nn) = zero
134  END DO
135  END DO
136 
137  !Accumulation
139  DO nn = 1, nmax(0)
141  END DO
142 
143  DO m = 1, mmax
145  DO nn = 1, nmax(m)
146  accumulation_x_s(m, nn) = accumulation_sp(rop_s(1, m), x_s(1, m, nn))
147  END DO
148  END DO
149 
150  return
151 
152 
153  else
154 ! Flow in and out of the boundaries
155 ! Use dt_prev for these integrations because adjust_dt may have changed the dt
156  DO l = 1, dimension_bc
157  IF (bc_defined(l)) THEN
158  i1 = bc_i_w(l)
159  i2 = bc_i_e(l)
160  j1 = bc_j_s(l)
161  j2 = bc_j_n(l)
162  k1 = bc_k_b(l)
163  k2 = bc_k_t(l)
164 
165 ! call Calc_mass_flux(I1, I2, J1, J2, K1, K2, BC_PLANE(L), U_g, V_g, W_g, ROP_g, fin, fout, IER)
166  IF(.NOT.added_mass) THEN
167  call calc_mass_fluxhr(i1, i2, j1, j2, k1, k2, bc_plane(l), flux_ge, flux_gn, flux_gt, fin, fout, ier)
168  ELSE
169  call calc_mass_fluxhr(i1, i2, j1, j2, k1, k2, bc_plane(l), flux_gse, flux_gsn, flux_gst, fin, fout, ier)
170  ENDIF
171  flux_out_g(l) = flux_out_g(l) + fout * dt_prev
172  flux_in_g(l) = flux_in_g(l) + fin * dt_prev
173 
174  DO m = 1, mmax
175 ! call Calc_mass_flux(I1, I2, J1, J2, K1, K2, BC_PLANE(L), U_s(1,m), V_s(1,m), W_s(1,m), ROP_s(1,m), fin, fout, IER)
176  IF(.NOT.added_mass .OR. m /= m_am) THEN
177  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)
178  ELSE
179  call calc_mass_fluxhr(i1, i2, j1, j2, k1, k2, bc_plane(l), flux_sse, flux_ssn, flux_sst, fin, fout, ier)
180  ENDIF
181  flux_out_s(l, m) = flux_out_s(l, m) + fout * dt_prev
182  flux_in_s(l, m) = flux_in_s(l, m) + fin * dt_prev
183  END DO
184 
185  ENDIF
186  END DO
187 
189  IF(species_eq(0))THEN
190  DO nn = 1, nmax(0)
191  DO l = 1, dimension_bc
192  IF (bc_defined(l)) THEN
193  i1 = bc_i_w(l)
194  i2 = bc_i_e(l)
195  j1 = bc_j_s(l)
196  j2 = bc_j_n(l)
197  k1 = bc_k_b(l)
198  k2 = bc_k_t(l)
199 ! call Calc_mass_flux_sp(I1, I2, J1, J2, K1, K2, BC_PLANE(L), U_g, V_g, W_g, ROP_g, X_g(1, N), fin, fout, IER)
200  IF(.NOT.added_mass) THEN
201  call calc_mass_flux_sphr(i1, i2, j1, j2, k1, k2, bc_plane(l), flux_ge, &
202  flux_gn, flux_gt, x_g(1, nn), discretize(7), fin, fout, ier)
203  ELSE
204  call calc_mass_flux_sphr(i1, i2, j1, j2, k1, k2, bc_plane(l), flux_gse, &
205  flux_gsn, flux_gst, x_g(1, nn), discretize(7), fin, fout, ier)
206  ENDIF
207  flux_out_x_g(l, nn) = flux_out_x_g(l, nn) + fout * dt_prev
208  flux_in_x_g(l, nn) = flux_in_x_g(l, nn) + fin * dt_prev
209  ENDIF
210  END DO
211 
212  integral_r_g(nn) = integral_r_g(nn) + (accumulation(r_gp(1,nn)) - &
213  accumulation_sp(rox_gc(1,nn), x_g(1,nn)) )* dt_prev
214  END DO
215  ENDIF
216 
217  DO m = 1, mmax
218  IF(species_eq(m))THEN
220  DO nn = 1, nmax(m)
221  DO l = 1, dimension_bc
222  IF (bc_defined(l)) THEN
223  i1 = bc_i_w(l)
224  i2 = bc_i_e(l)
225  j1 = bc_j_s(l)
226  j2 = bc_j_n(l)
227  k1 = bc_k_b(l)
228  k2 = bc_k_t(l)
229 ! call Calc_mass_flux_sp(I1, I2, J1, J2, K1, K2, BC_PLANE(L), &
230 ! U_s(1,M), V_s(1,M), W_s(1,M), ROP_s(1,M), X_s(1, M, N), fin, fout, &
231  IF(.NOT.added_mass .OR. m /= m_am) THEN
232  call calc_mass_flux_sphr(i1, i2, j1, j2, k1, k2, bc_plane(l), &
233  flux_se(1,m), flux_sn(1,m), flux_st(1,m), x_s(1, m, nn), discretize(7), fin, fout, &
234  ier)
235  ELSE
236  call calc_mass_flux_sphr(i1, i2, j1, j2, k1, k2, bc_plane(l), &
237  flux_sse, flux_ssn, flux_sst, x_s(1, m, nn), discretize(7), fin, fout, ier)
238  ENDIF
239  flux_out_x_s(l, m, nn) = flux_out_x_s(l, m, nn) + fout * dt_prev
240  flux_in_x_s(l, m, nn) = flux_in_x_s(l, m, nn) + fin * dt_prev
241  ENDIF
242  END DO
243 
244  integral_r_s(m, nn) = integral_r_s(m, nn) + (accumulation(r_sp(1,m,nn)) - &
245  accumulation_sp(rox_sc(1,m,nn), x_s(1,m,nn)) )* dt_prev
246  END DO
247  ENDIF
248  END DO
249 
250  endif
251 
252 
253  if (time >= report_time)then
254 
255  CALL start_log
256  IF(dmp_log)WRITE(unit_log, '(/A,G12.5,A,G12.5)') 'Mass balance for interval ', start_time, ' to ', time
257 
258  accumulation_old = accumulation_g
260  accumulation_delta = accumulation_g - accumulation_old - integral_sum_r_g
261  IF(dmp_log)WRITE(unit_log, '(A)') 'Total Fluid Accumulation (g)'
262  IF(dmp_log)WRITE(unit_log, '(4(A,G12.5))') ' Old = ', accumulation_old, ', New = ', &
263  accumulation_g, ', Production = ', integral_sum_r_g, ', net accu(New - Old - Production) = ', accumulation_delta
264 
265  IF(dmp_log)WRITE(unit_log, '(A)') 'Integral of boundary flux (g)'
266  IF(dmp_log)WRITE(unit_log, '(A, T8, A, T21, A, T34, A)')' BC#', 'in', 'out', '(in - out)'
267  flux = zero
268  flux_in_tot = zero
269  flux_out_tot = zero
270  DO l = 1, dimension_bc
271  if(flux_out_g(l) /= zero .OR. flux_in_g(l) /= zero) then
272  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))
273  endif
274  flux = flux + flux_in_g(l) - flux_out_g(l)
275  flux_in_tot = flux_in_tot + flux_in_g(l)
276  flux_out_tot = flux_out_tot + flux_out_g(l)
277  END DO
278  if((flux - accumulation_delta) /= zero) then
279  error_percent = undefined
280  if(flux_in_tot /= zero) error_percent = (flux - accumulation_delta)*100./flux_in_tot
281  else
282  error_percent = zero
283  endif
284  IF(dmp_log)WRITE(unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
285  IF(dmp_log)WRITE(unit_log, '(A, G12.5, A, G12.5)')'Error (net influx - net accu) = ', &
286  (flux - accumulation_delta), ' %Error_g = ', error_percent
287 
288  DO m = 1, mmax
289  accumulation_old = accumulation_s(m)
290  accumulation_s(m) = accumulation(rop_s(1, m))
291  accumulation_delta = accumulation_s(m) - accumulation_old - integral_sum_r_s(m)
292  IF(dmp_log)WRITE(unit_log, '(/A, I1, A)') 'Total Solids-', m, ' Accumulation (g)'
293  IF(dmp_log)WRITE(unit_log, '(4(A,G12.5))') ' Old = ', accumulation_old, ', New = ', &
294  accumulation_s(m), ', Production = ', integral_sum_r_s(m), ', net accu(New - Old - Production) = ', accumulation_delta
295 
296  IF(dmp_log)WRITE(unit_log, '(A)') 'Integral of boundary flux (g)'
297  IF(dmp_log)WRITE(unit_log, '(A, T8, A, T21, A, T34, A)')' BC#', 'in', 'out', '(in - out)'
298  flux = zero
299  flux_in_tot = zero
300  flux_out_tot = zero
301  DO l = 1, dimension_bc
302  if(flux_out_s(l,m) /= zero .OR. flux_in_s(l,m) /= zero) then
303  IF(dmp_log) THEN
304  WRITE(unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')l, &
305  flux_in_s(l,m), flux_out_s(l,m), (flux_in_s(l,m)-flux_out_s(l,m))
306  ENDIF
307  endif
308  flux = flux + flux_in_s(l,m) - flux_out_s(l,m)
309  flux_in_tot = flux_in_tot + flux_in_s(l,m)
310  flux_out_tot = flux_out_tot + flux_out_s(l,m)
311  END DO
312  if((flux - accumulation_delta) /= zero) then
313  if(flux_in_tot /= zero) then
314  error_percent = (flux - accumulation_delta)*100./flux_in_tot
315  else
316  error_percent = (flux - accumulation_delta)*100./accumulation_old
317  endif
318  else
319  error_percent = zero
320  endif
321  IF(dmp_log)WRITE(unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
322  IF(dmp_log)WRITE(unit_log, '(A, G12.5, A, I1, A, G12.5)')'Error (net influx - net accu) = ', &
323  (flux - accumulation_delta), ' %Error_',m,' = ', error_percent
324  END DO
325 
326 
327  IF(species_eq(0) )THEN
328  DO nn = 1, nmax(0)
329 
330  IF(dmp_log)WRITE(unit_log, '(/A,I2)') 'Gas species - ', nn
331 
332  accumulation_old = accumulation_x_g(nn)
334  accumulation_delta = accumulation_x_g(nn) - accumulation_old - integral_r_g(nn)
335  IF(dmp_log)WRITE(unit_log, '(A)') 'Species Accumulation (g)'
336  IF(dmp_log)WRITE(unit_log, '(4(A,G12.5))') ' Old = ', accumulation_old, ', New = ', &
337  accumulation_x_g(nn), ', Production = ', integral_r_g(nn), ', net accu(New - Old - Production) = ', accumulation_delta
338 
339  IF(dmp_log)WRITE(unit_log, '(A)') 'Integral of boundary flux (g)'
340  IF(dmp_log)WRITE(unit_log, '(A, T8, A, T21, A, T34, A)')' BC#', 'in', 'out', '(in - out)'
341  flux = zero
342  flux_in_tot = zero
343  flux_out_tot = zero
344  DO l = 1, dimension_bc
345  if(flux_out_x_g(l, nn) /= zero .OR. flux_in_x_g(l, nn) /= zero) then
346  IF(dmp_log)WRITE(unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')l, flux_in_x_g(l, nn), flux_out_x_g(l, nn), &
347  (flux_in_x_g(l, nn)-flux_out_x_g(l, nn))
348  endif
349  flux = flux + flux_in_x_g(l, nn) - flux_out_x_g(l, nn)
350  flux_in_tot = flux_in_tot + flux_in_x_g(l, nn)
351  flux_out_tot = flux_out_tot + flux_out_x_g(l, nn)
352  END DO
353  if((flux - accumulation_delta) /= zero) then
354  error_percent = undefined
355  if(flux_in_tot /= zero) then
356  error_percent = (flux - accumulation_delta)*100./flux_in_tot
357  elseif (integral_r_g(nn) /= zero) then
358  error_percent = (flux - accumulation_delta)*100./integral_r_g(nn)
359  endif
360  else
361  error_percent = zero
362  endif
363  IF(dmp_log)WRITE(unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
364  IF(dmp_log)WRITE(unit_log, '(A, G12.5, A, I2, A, G12.5)')'Error (net influx - net accu) = ', &
365  (flux - accumulation_delta), ' %Error_g(',nn,') = ', error_percent
366 
367  END DO
368  ENDIF
369 
370  DO m = 1, mmax
371  IF(species_eq(m) )THEN
372  DO nn = 1, nmax(m)
373 
374  IF(dmp_log)WRITE(unit_log, '(/A,I1, A, I2)') 'Solids-', m, ' species - ', nn
375 
376  accumulation_old = accumulation_x_s(m,nn)
377  accumulation_x_s(m,nn) = accumulation_sp(rop_s(1,m), x_s(1, m, nn))
378  accumulation_delta = accumulation_x_s(m,nn) - accumulation_old - integral_r_s(m,nn)
379  IF(dmp_log)WRITE(unit_log, '(A)') 'Species Accumulation (g)'
380  IF(dmp_log)WRITE(unit_log, '(4(A,G12.5))') ' Old = ', accumulation_old, ', New = ', &
381  accumulation_x_s(m,nn), ', Production = ', integral_r_s(m,nn), &
382  ', net accu(New - Old - Production) = ', accumulation_delta
383 
384  IF(dmp_log)WRITE(unit_log, '(A)') 'Integral of boundary flux (g)'
385  IF(dmp_log)WRITE(unit_log, '(A, T8, A, T21, A, T34, A)')' BC#', 'in', 'out', '(in - out)'
386  flux = zero
387  flux_in_tot = zero
388  flux_out_tot = zero
389  DO l = 1, dimension_bc
390  if(flux_out_x_s(l, m, nn) /= zero .OR. flux_in_x_s(l, m, nn) /= zero) then
391  IF(dmp_log)WRITE(unit_log, '(2X, I5, 1X, 3(G12.5, 1x))')l, flux_in_x_s(l, m, nn), flux_out_x_s(l, m, nn), &
392  (flux_in_x_s(l, m, nn)-flux_out_x_s(l, m, nn))
393  endif
394  flux = flux + flux_in_x_s(l, m, nn) - flux_out_x_s(l, m, nn)
395  flux_in_tot = flux_in_tot + flux_in_x_s(l, m, nn)
396  flux_out_tot = flux_out_tot + flux_out_x_s(l, m, nn)
397  END DO
398  if((flux - accumulation_delta) /= zero) then
399  error_percent = undefined
400  if(flux_in_tot /= zero) then
401  error_percent = (flux - accumulation_delta)*100./flux_in_tot
402  elseif (integral_r_s(m,nn) /= zero)then
403  error_percent = (flux - accumulation_delta)*100./integral_r_s(m,nn)
404  endif
405  else
406  error_percent = zero
407  endif
408  IF(dmp_log)WRITE(unit_log, '(2X, A, 1X, 3(G12.5, 1x))')'Total', flux_in_tot, flux_out_tot, flux
409  IF(dmp_log)WRITE(unit_log, '(A, G12.5, A, I1, A, I2, A, G12.5)')'Error (net influx - net accu) = ', &
410  (flux - accumulation_delta), ' %Error_',m,'(',nn,') = ', error_percent
411 
412  END DO
413  ENDIF
414  END DO
415  IF(dmp_log)WRITE(unit_log, '(/)')
416  CALL end_log
417 
418  start_time = time
420 
421  !initialize flux and reaction rate arrays
422  DO l = 1, dimension_bc
423  flux_out_g(l) = zero
424  flux_in_g(l) = zero
425  DO nn = 1, nmax(0)
426  flux_out_x_g(l, nn) = zero
427  flux_in_x_g(l, nn) = zero
428  END DO
429 
430  DO m = 1, mmax
431  flux_out_s(l, m) = zero
432  flux_in_s(l, m) = zero
433  DO nn = 1, nmax(m)
434  flux_out_x_s(l, m, nn) = zero
435  flux_in_x_s(l, m, nn) = zero
436  END DO
437  END DO
438 
439  END DO
440 
442  DO nn = 1, nmax(0)
443  integral_r_g(nn) = zero
444  END DO
445 
446  DO m = 1, mmax
448  DO nn = 1, nmax(m)
449  integral_r_s(m, nn) = zero
450  END DO
451  END DO
452 
453 
454  endif
455 
456 
457  RETURN
458  END SUBROUTINE check_mass_balance
459 
460 
461 
462 
463 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
464 ! C
465 ! Module name: Calc_mass_flux C
466 ! Purpose: Calculate the mass flux (g/s) across a plane C
467 ! C
468 ! Author: M. Syamlal Date: 9-Dec-02 C
469 ! Reviewer: Date: C
470 ! C
471 ! Literature/Document References: C
472 ! C
473 ! Variables referenced: C
474 ! Variables modified: C
475 ! C
476 ! C
477 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
478 !
479  SUBROUTINE calc_mass_flux(I1, I2, J1, J2, K1, K2, Plane, U, V, W, ROP, flux_in, flux_out, IER)
480  USE param, only: dimension_3
481  USE param1, only: one
482  IMPLICIT NONE
483 
484 !
485 ! Starting and ending indices (Make sure these represent a plane and NOT a volume)
486  INTEGER :: I1, I2, J1, J2, K1, K2
487 
488 ! For each (scalar) cell the plane (W, east, south, N, bottom, T) to be used for flux calc
489  Character :: Plane
490 
491 ! Components of Velocity
492  DOUBLE PRECISION :: U(dimension_3), V(dimension_3), W(dimension_3)
493 
494 ! Macroscopic density
495  DOUBLE PRECISION :: ROP(dimension_3)
496 
497 ! Species Mass fraction
498  DOUBLE PRECISION :: Xn(dimension_3)
499 
500 ! flux across the plane (composed of many cell faces)
501 ! flux_in: flux from the scalar cell into the rest of the domain
502 ! flux_out: flux into the scalar cell from the rest of the domain
503  DOUBLE PRECISION :: flux_in, flux_out
504 
505  INTEGER :: IER
506 
507  xn = one
508  call calc_mass_flux_sp(i1, i2, j1, j2, k1, k2, plane, u, v, w, rop, xn, flux_in, flux_out, ier)
509  return
510  end SUBROUTINE calc_mass_flux
511 
512 
513 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
514 ! C
515 ! Module name: Calc_mass_flux_sp C
516 ! Purpose: Calculate the species mass flux (g/s) across a plane C
517 ! C
518 ! Author: M. Syamlal Date: 9-Dec-02 C
519 ! Reviewer: Date: C
520 ! C
521 ! Literature/Document References: C
522 ! C
523 ! Variables referenced: C
524 ! Variables modified: C
525 ! C
526 ! C
527 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
528 !
529  SUBROUTINE calc_mass_flux_sp(I1, I2, J1, J2, K1, K2, Plane, U, V, W, ROP, Xn, flux_in, flux_out, IER)
530  USE compar
531  USE functions
532  USE geometry
533  USE indices
534  USE mpi_utility
535  USE param, only: dimension_3
536  USE param1, only: zero
537  USE physprop
538  IMPLICIT NONE
539 
540 !
541 ! Starting and ending indices (Make sure these represent a plane and NOT a volume)
542  INTEGER :: I1, I2, J1, J2, K1, K2
543 
544 ! For each (scalar) cell the plane (W, east, south, N, bottom, T) to be used for flux calc
545  Character :: Plane
546 
547 ! Components of Velocity
548  DOUBLE PRECISION :: U(dimension_3), V(dimension_3), W(dimension_3)
549 
550 ! Macroscopic density
551  DOUBLE PRECISION :: ROP(dimension_3)
552 
553 ! Species Mass fraction
554  DOUBLE PRECISION :: Xn(dimension_3)
555 
556 ! flux across the plane (composed of many cell faces)
557 ! flux_in: flux from the scalar cell into the rest of the domain
558 ! flux_out: flux into the scalar cell from the rest of the domain
559  DOUBLE PRECISION :: flux_in, flux_out
560 
561  INTEGER :: IER
562 !
563 ! Indices
564  INTEGER I, J, K, IJK
565 
566  flux_in = zero
567  flux_out = zero
568  ier = 0
569 
570 
571  DO k = k1, k2
572  DO j = j1, j2
573  DO i = i1, i2
574  IF(.NOT.is_on_mype_owns(i, j, k)) cycle
575  IF(dead_cell_at(i,j,k)) cycle
576 
577  SELECT CASE (trim(plane))
578  CASE ('W')
579  ijk = im_of(funijk(i,j,k))
580  IF(u(ijk) > zero)THEN
581  flux_out = flux_out + ayz(ijk) * u(ijk) * rop(ijk) * xn(ijk)
582  ELSE
583  flux_in = flux_in - ayz(ijk) * u(ijk) * rop(ip_of(ijk)) * xn(ip_of(ijk))
584  ENDIF
585 
586  CASE ('E')
587  ijk = funijk(i,j,k)
588  IF(u(ijk) > zero)THEN
589  flux_in = flux_in + ayz(ijk) * u(ijk) * rop(ijk) * xn(ijk)
590  ELSE
591  flux_out = flux_out - ayz(ijk) * u(ijk) * rop(ip_of(ijk)) * xn(ip_of(ijk))
592  ENDIF
593 
594  CASE ('S')
595  ijk = jm_of(funijk(i,j,k))
596  IF(v(ijk) > zero)THEN
597  flux_out = flux_out + axz(ijk) * v(ijk) * rop(ijk) * xn(ijk)
598  ELSE
599  flux_in = flux_in - axz(ijk) * v(ijk) * rop(jp_of(ijk)) * xn(jp_of(ijk))
600  ENDIF
601 
602  CASE ('N')
603  ijk = funijk(i,j,k)
604  IF(v(ijk) > zero)THEN
605  flux_in = flux_in + axz(ijk) * v(ijk) * rop(ijk) * xn(ijk)
606  ELSE
607  flux_out = flux_out - axz(ijk) * v(ijk) * rop(jp_of(ijk)) * xn(jp_of(ijk))
608  ENDIF
609 
610  CASE ('B')
611  ijk = km_of(funijk(i,j,k))
612  IF(w(ijk) > zero)THEN
613  flux_out = flux_out + axy(ijk) * w(ijk) * rop(ijk) * xn(ijk)
614  ELSE
615  flux_in = flux_in - axy(ijk) * w(ijk) * rop(kp_of(ijk)) * xn(kp_of(ijk))
616  ENDIF
617 
618  CASE ('T')
619  ijk = funijk(i,j,k)
620  IF(w(ijk) > zero)THEN
621  flux_in = flux_in + axy(ijk) * w(ijk) * rop(ijk) * xn(ijk)
622  ELSE
623  flux_out = flux_out - axy(ijk) * w(ijk) * rop(kp_of(ijk)) * xn(kp_of(ijk))
624  ENDIF
625 
626  CASE DEFAULT
627 ! IER = 1
628 ! CALL START_LOG
629 ! IF(DMP_LOG)WRITE (UNIT_LOG, '(A, A1)' ) 'From: Calc_mass_flux, Unknown Plane: ', Plane
630 ! CALL MFIX_EXIT(myPE)
631 
632  END SELECT
633  ENDDO
634  ENDDO
635  ENDDO
636 
637  call global_all_sum(flux_in)
638  call global_all_sum(flux_out)
639 
640  return
641  end SUBROUTINE calc_mass_flux_sp
642 
643 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
644 ! C
645 ! Module name: Calc_mass_fluxHR C
646 ! Purpose: Calculate the mass flux (g/s) across a plane. HR version. C
647 ! C
648 ! Author: M. Syamlal Date: 26-Jul-05 C
649 ! Reviewer: Date: C
650 ! C
651 ! Literature/Document References: C
652 ! C
653 ! Variables referenced: C
654 ! Variables modified: C
655 ! C
656 ! C
657 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
658 !
659  SUBROUTINE calc_mass_fluxhr(I1, I2, J1, J2, K1, K2, Plane, Flux_E, Flux_N, Flux_T, flux_in, flux_out, IER)
660  USE param
661  USE param1
662  IMPLICIT NONE
663 
664 !
665 ! Starting and ending indices (Make sure these represent a plane and NOT a volume)
666  INTEGER :: I1, I2, J1, J2, K1, K2
667 
668 ! For each (scalar) cell the plane (W, east, south, N, bottom, T) to be used for flux calc
669  Character :: Plane
670 
671 ! Components of flux
672  DOUBLE PRECISION :: Flux_E(dimension_3), Flux_N(dimension_3), Flux_T(dimension_3)
673 
674 ! Species Mass fraction
675  DOUBLE PRECISION :: Xn(dimension_3)
676 
677 ! flux across the plane (composed of many cell faces)
678 ! flux_in: flux from the scalar cell into the rest of the domain
679 ! flux_out: flux into the scalar cell from the rest of the domain
680  DOUBLE PRECISION :: flux_in, flux_out
681 
682  INTEGER :: IER
683 
684  xn = one
685  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)
686  return
687  end SUBROUTINE calc_mass_fluxhr
688 
689 
690 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
691 ! C
692 ! Module name: Calc_mass_flux_spHR C
693 ! Purpose: Calculate the species mass flux (g/s) across a plane. HR version C
694 ! C
695 ! Author: M. Syamlal Date: 26-Jul-05 C
696 ! Reviewer: Date: C
697 ! C
698 ! Literature/Document References: C
699 ! C
700 ! Variables referenced: C
701 ! Variables modified: C
702 ! C
703 ! C
704 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
705 !
706  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)
707  USE param
708  USE param1
709  USE geometry
710  USE physprop
711  USE indices
712  USE compar
713  USE mpi_utility
714  USE functions
715  USE xsi, only: calc_xsi
716  IMPLICIT NONE
717 
718 !
719 ! Starting and ending indices (Make sure these represent a plane and NOT a volume)
720  INTEGER :: I1, I2, J1, J2, K1, K2
721 
722 ! For each (scalar) cell the plane (W, east, south, N, bottom, T) to be used for flux calc
723  Character :: Plane
724 
725 ! Components of flux
726  DOUBLE PRECISION :: Flux_E(dimension_3), Flux_N(dimension_3), Flux_T(dimension_3)
727 
728 ! Species Mass fraction
729  DOUBLE PRECISION :: Xn(dimension_3)
730 
731 !
732 ! Discretizationindex
733  INTEGER Disc
734 
735 ! flux across the plane (composed of many cell faces)
736 ! flux_in: flux from the scalar cell into the rest of the domain
737 ! flux_out: flux into the scalar cell from the rest of the domain
738  DOUBLE PRECISION :: flux_in, flux_out
739 
740 ! Flux and Xn at the boundary plane
741  DOUBLE PRECISION :: Flux, PHI_HO
742 
743  INTEGER :: IER
744 !
745 ! Indices
746  INTEGER I, J, K, IJK
747 
748  DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t
749 
750 ! Unlike other calls to calc_xsi, the fluxes rather than the velocities are passed.
751 ! This should work fine because the velocities are used to determine the upwind bias,
752 ! for which fluxes should suffice. The flag incr is not needed and is set to 0.
753  CALL calc_xsi (disc, xn, flux_e, flux_n, flux_t, xsi_e, xsi_n, xsi_t, 0)
754 
755  flux_in = zero
756  flux_out = zero
757  ier = 0
758 
759  DO k = k1, k2
760  DO j = j1, j2
761  DO i = i1, i2
762  IF(.NOT.is_on_mype_owns(i, j, k)) cycle
763  IF(dead_cell_at(i,j,k)) cycle
764 
765  SELECT CASE (trim(plane))
766  CASE ('W')
767  ijk = im_of(funijk(i,j,k))
768  phi_ho = xsi_e(ijk)*xn(ip_of(ijk))+(1.0-xsi_e(ijk))*xn(ijk)
769  flux = flux_e(ijk)
770 
771  CASE ('E')
772  ijk = funijk(i,j,k)
773  phi_ho = xsi_e(ijk)*xn(ip_of(ijk))+(1.0-xsi_e(ijk))*xn(ijk)
774  flux = -flux_e(ijk)
775 
776  CASE ('S')
777  ijk = jm_of(funijk(i,j,k))
778  phi_ho = xsi_n(ijk)*xn(jp_of(ijk))+(1.0-xsi_n(ijk))*xn(ijk)
779  flux = flux_n(ijk)
780 
781  CASE ('N')
782  ijk = funijk(i,j,k)
783  phi_ho = xsi_n(ijk)*xn(jp_of(ijk))+(1.0-xsi_n(ijk))*xn(ijk)
784  flux = -flux_n(ijk)
785 
786  CASE ('B')
787  ijk = km_of(funijk(i,j,k))
788  phi_ho = xsi_t(ijk)*xn(kp_of(ijk))+(1.0-xsi_t(ijk))*xn(ijk)
789  flux = flux_t(ijk)
790 
791  CASE ('T')
792  ijk = funijk(i,j,k)
793  phi_ho = xsi_t(ijk)*xn(kp_of(ijk))+(1.0-xsi_t(ijk))*xn(ijk)
794  flux = -flux_t(ijk)
795 
796  CASE DEFAULT
797  phi_ho = zero
798  flux = zero
799 
800  END SELECT
801 
802  IF(flux > zero)THEN
803  flux_out = flux_out + flux * phi_ho
804  ELSE
805  flux_in = flux_in - flux * phi_ho
806  ENDIF
807 
808  ENDDO
809  ENDDO
810  ENDDO
811 
812  call global_all_sum(flux_in)
813  call global_all_sum(flux_out)
814 
815  return
816  end SUBROUTINE calc_mass_flux_sphr
817 
818 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
819 ! C
820 ! Module name: Accumulation(ro) C
821 ! Purpose: Intergrate density accumulation over the entire domain C
822 ! C
823 ! Author: M. Syamlal Date: 30-DEC-02 C
824 ! Local variables: None C
825 ! C
826 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
827 !
828  DOUBLE PRECISION FUNCTION accumulation(ro)
829 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98
830 !...Switches: -xf
831 !
832 ! Include param.inc file to specify parameter values
833 !
834 !-----------------------------------------------
835 ! M o d u l e s
836 !-----------------------------------------------
837  USE param
838  USE param1
839  IMPLICIT NONE
840 !-----------------------------------------------
841 ! G l o b a l P a r a m e t e r s
842 !-----------------------------------------------
843 !-----------------------------------------------
844 ! D u m m y A r g u m e n t s
845 !-----------------------------------------------
846 ! density distributiom
847  DOUBLE PRECISION :: RO(dimension_3)
848 ! Species Mass fraction
849  DOUBLE PRECISION :: Xn(dimension_3)
850 
851 !-----------------------------------------------
852 
853  xn = one
854  accumulation = accumulation_sp(ro, xn)
855 
856  RETURN
857  END FUNCTION accumulation
858 
859 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
860 ! C
861 ! Module name: Accumulation_sp(ro, Xn) C
862 ! Purpose: Intergrate density accumulation over the entire domain C
863 ! C
864 ! Author: M. Syamlal Date: 30-DEC-02 C
865 ! Local variables: None C
866 ! C
867 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
868 !
869  DOUBLE PRECISION FUNCTION accumulation_sp(ro, Xn)
870 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98
871 !...Switches: -xf
872 !
873 ! Include param.inc file to specify parameter values
874 !
875 !-----------------------------------------------
876 ! M o d u l e s
877 !-----------------------------------------------
878  USE param
879  USE param1
880  USE parallel
881  USE physprop
882  USE geometry
883  USE indices
884  USE compar
885  USE mpi_utility
886  USE functions
887  IMPLICIT NONE
888 !-----------------------------------------------
889 ! G l o b a l P a r a m e t e r s
890 !-----------------------------------------------
891 !-----------------------------------------------
892 ! D u m m y A r g u m e n t s
893 !-----------------------------------------------
894 ! Indices
895  INTEGER IJK
896 
897 ! density distributiom
898  DOUBLE PRECISION :: RO(dimension_3)
899 ! Species Mass fraction
900  DOUBLE PRECISION :: Xn(dimension_3)
901 
902  DOUBLE PRECISION SUM
903 !
904 !-----------------------------------------------
905 
906  sum = zero
907 !
908 !$omp parallel do default(none) private(IJK) shared(ijkstart3, ijkend3, i_of, j_of, k_of, ro, xn, vol) &
909 !$omp reduction(+:SUM)
910  DO ijk = ijkstart3, ijkend3
911  IF(.NOT.is_on_mype_wobnd(i_of(ijk), j_of(ijk), k_of(ijk))) cycle
912  IF (fluid_at(ijk)) sum = sum + ro(ijk) * xn(ijk) * vol(ijk)
913  END DO
914 
915  call global_all_sum(sum)
916  accumulation_sp = sum
917 
918  RETURN
919  END FUNCTION accumulation_sp
920 
921 
922 
923  DOUBLE PRECISION FUNCTION check_conservation (Phi, A_m, B_m, IJK)
924  USE param
925  USE param1
926  USE geometry
927  USE indices
928  USE compar
929  USE functions
930  IMPLICIT NONE
931 !
932 ! Septadiagonal matrix A_m
933  DOUBLE PRECISION A_m(dimension_3, -3:3, 0:dimension_m)
934 !
935 ! Vector b_m
936  DOUBLE PRECISION B_m(dimension_3, 0:dimension_m)
937 
938  DOUBLE PRECISION Phi(dimension_3)
939 
940 ! Indices
941  INTEGER IJK, IMJK, IJMK, IJKM, IPJK, IJPK, IJKP
942 
943  imjk = im_of(ijk)
944  ijmk = jm_of(ijk)
945  ijkm = km_of(ijk)
946  ipjk = ip_of(ijk)
947  ijpk = jp_of(ijk)
948  ijkp = kp_of(ijk)
949 
950  check_conservation = phi(ijk) * a_m(ijk, 0, 0) &
951  + phi(ipjk) * a_m(ijk, east, 0) &
952  + phi(imjk) * a_m(ijk, west, 0) &
953  + phi(ijpk) * a_m(ijk, north, 0) &
954  + phi(ijmk) * a_m(ijk, south, 0) &
955  + phi(ijkp) * a_m(ijk, top, 0) &
956  + phi(ijkm) * a_m(ijk, bottom, 0) &
957  - b_m(ijk,0)
958 
959  END FUNCTION check_conservation
960 
961  END MODULE check
integer, dimension(dimension_bc) bc_k_b
Definition: bc_mod.f:70
double precision accumulation_g
Definition: check_mod.f:25
double precision, dimension(:), allocatable flux_ge
Definition: mflux_mod.f:13
integer, parameter dim_n_g
Definition: param_mod.f:69
double precision, dimension(dimension_b!c) flux_in_g
Definition: check_mod.f:27
logical dmp_log
Definition: funits_mod.f:6
double precision function check_conservation(Phi, A_m, B_m, IJK)
Definition: check_mod.f:924
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
double precision report_time
Definition: check_mod.f:23
integer ijkend3
Definition: compar_mod.f:80
double precision dt_prev
Definition: run_mod.f:230
double precision, dimension(:), allocatable flux_gst
Definition: mflux_mod.f:32
double precision, dimension(:,:), allocatable flux_st
Definition: mflux_mod.f:24
double precision, parameter one
Definition: param1_mod.f:29
integer dimension_3
Definition: param_mod.f:11
subroutine check_mass_balance(init)
Definition: check_mod.f:38
double precision, dimension(:), allocatable axy
Definition: geometry_mod.f:210
double precision, dimension(dime!nsion_bc, dim_m) flux_in_s
Definition: check_mod.f:32
integer, dimension(dimension_bc) bc_i_w
Definition: bc_mod.f:54
logical added_mass
Definition: run_mod.f:91
integer, dimension(dimension_bc) bc_j_n
Definition: bc_mod.f:66
Definition: rxns_mod.f:1
logical, dimension(0:dim_m) species_eq
Definition: run_mod.f:115
double precision, dimension(:), allocatable flux_ssn
Definition: mflux_mod.f:31
subroutine calc_mass_fluxhr(I1, I2, J1, J2, K1, K2, Plane, Flux_E
Definition: check_mod.f:660
double precision start_time
Definition: check_mod.f:23
double precision integral_sum_r_g
Definition: check_mod.f:26
integer, parameter dim_m
Definition: param_mod.f:67
double precision flux_in_x_g
Definition: check_mod.f:28
double precision, dimension(:), allocatable sum_r_g
Definition: rxns_mod.f:28
double precision, dimension(:,:), allocatable sum_r_s
Definition: rxns_mod.f:35
double precision, dimension(dim_m) accumulation_s
Definition: check_mod.f:30
integer, parameter dimension_bc
Definition: param_mod.f:61
subroutine calc_xsi(DISCR, PHI, U, V, W, XSI_E, XSI_N, XSI_T, incr)
Definition: xsi_mod.f:23
double precision, dimension(:,:,:), allocatable rox_sc
Definition: rxns_mod.f:33
double precision, parameter undefined
Definition: param1_mod.f:18
integer east
Definition: param_mod.f:29
double precision, dimension(:), allocatable ayz
Definition: geometry_mod.f:206
subroutine calc_mass_flux_sphr(I1, I2, J1, J2, K1, K2, Plane, Flu
Definition: check_mod.f:707
Definition: xsi_mod.f:15
double precision report_mass_balance_dt
Definition: output_mod.f:27
character, dimension(dimension_bc) bc_plane
Definition: bc_mod.f:217
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
double precision, dimension(dim_m) integral_sum_r_s
Definition: check_mod.f:31
double precision, dimension(:), allocatable flux_gse
Definition: mflux_mod.f:28
integer, dimension(dimension_bc) bc_k_t
Definition: bc_mod.f:74
double precision function accumulation_sp(ro, Xn)
Definition: check_mod.f:870
double precision function accumulation(ro)
Definition: check_mod.f:829
integer mmax
Definition: physprop_mod.f:19
logical, dimension(:,:,:), allocatable dead_cell_at
Definition: compar_mod.f:127
double precision, dimension(:,:,:), allocatable x_s
Definition: fldvar_mod.f:78
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
double precision, dimension(:), allocatable flux_gn
Definition: mflux_mod.f:15
double precision, dimension(:,:), allocatable x_g
Definition: fldvar_mod.f:75
integer, dimension(dimension_bc) bc_j_s
Definition: bc_mod.f:62
integer m_am
Definition: run_mod.f:94
double precision, dimension(:,:,:), allocatable r_sp
Definition: rxns_mod.f:31
integer north
Definition: param_mod.f:37
logical, dimension(dimension_bc) bc_defined
Definition: bc_mod.f:207
integer south
Definition: param_mod.f:41
integer, parameter unit_log
Definition: funits_mod.f:21
Definition: run_mod.f:13
double precision, dimension(:), allocatable axz
Definition: geometry_mod.f:208
Definition: param_mod.f:2
double precision, dimension(dim_n_g) integral_r_g
Definition: check_mod.f:26
double precision, dimension(dim_m, di!m_n_s) integral_r_s
Definition: check_mod.f:31
double precision, dimension(:), allocatable flux_gt
Definition: mflux_mod.f:17
double precision, dimension(:,:), allocatable rox_gc
Definition: rxns_mod.f:26
integer, dimension(0:dim_m) nmax
Definition: physprop_mod.f:119
integer ijkstart3
Definition: compar_mod.f:80
integer west
Definition: param_mod.f:33
double precision, dimension(dim_n_g) accumulation_x_g
Definition: check_mod.f:25
integer, parameter dim_n_s
Definition: param_mod.f:71
integer top
Definition: param_mod.f:45
subroutine start_log
Definition: machine_mod.f:182
integer, dimension(dim_eqs) discretize
Definition: run_mod.f:67
double precision, dimension(:,:), allocatable rop_s
Definition: fldvar_mod.f:51
subroutine calc_mass_flux(I1, I2, J1, J2, K1, K2, Plane, U, V, W,
Definition: check_mod.f:480
double precision, dimension(:,:), allocatable flux_se
Definition: mflux_mod.f:22
double precision, dimension(:), allocatable flux_gsn
Definition: mflux_mod.f:30
double precision, dimension(dimension_bc, dim_n_g) flux_out_x_g
Definition: check_mod.f:28
double precision, dimension(dim_m,!dim_n_s) accumulation_x_s
Definition: check_mod.f:30
double precision time
Definition: run_mod.f:45
double precision, dimension(:), allocatable vol
Definition: geometry_mod.f:212
integer dimension_m
Definition: param_mod.f:18
subroutine calc_mass_flux_sp(I1, I2, J1, J2, K1, K2, Plane, U, V,
Definition: check_mod.f:530
double precision, dimension(:), allocatable rop_g
Definition: fldvar_mod.f:38
double precision, dimension(:,:), allocatable flux_sn
Definition: mflux_mod.f:20
integer bottom
Definition: param_mod.f:49
double precision, dimension(dimension_bc, dim_m, dim_n_s) flux_out_x_s
Definition: check_mod.f:33
integer, dimension(dimension_bc) bc_i_e
Definition: bc_mod.f:58
double precision, dimension(dimension_bc) flux_out_g
Definition: check_mod.f:27
double precision, parameter zero
Definition: param1_mod.f:27
subroutine end_log
Definition: machine_mod.f:208
double precision, dimension(:), allocatable flux_sst
Definition: mflux_mod.f:33
double precision, dimension(:,:), allocatable r_gp
Definition: rxns_mod.f:24
Definition: bc_mod.f:23
double precision, dimension(dimension_bc, dim_m) flux_out_s
Definition: check_mod.f:32
double precision, dimension(:), allocatable flux_sse
Definition: mflux_mod.f:29