MFIX  2016-1
calc_resid.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: CALC_RESID_C C
4 ! Purpose: Calculate residuals for continuity equations C
5 ! C
6 ! Author: M. Syamlal Date: 21-MAY-96 C
7 ! Reviewer: Date: C
8 ! C
9 ! Literature/Document References: C
10 ! Variables referenced: C
11 ! Variables modified: C
12 ! Local variables: C
13 ! C
14 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
15 
16  SUBROUTINE calc_resid_c(VAR, A_M, B_M, M, NUM, DEN, &
17  resid, max_resid, ijk_resid)
18 
19 !-----------------------------------------------
20 ! Modules
21 !-----------------------------------------------
22  USE param
23  USE param1, ONLY: zero, one, undefined
24  USE parallel
25  USE geometry
26  USE indices
27  USE compar
28  USE mpi_utility
29  USE run
30  IMPLICIT NONE
31 !-----------------------------------------------
32 ! Dummy arguments
33 !-----------------------------------------------
34 ! Variable being evaluated
35  DOUBLE PRECISION, INTENT(IN) :: Var(dimension_3)
36 ! Septadiagonal matrix A_m
37  DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3, 0:dimension_m)
38 ! Vector b_m
39  DOUBLE PRECISION, INTENT(IN) :: B_m(dimension_3, 0:dimension_m)
40 ! Phase index
41  INTEGER, INTENT(IN) :: M
42 ! Numerator and denominator
43  DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
44 ! Average value of Residual
45  DOUBLE PRECISION, INTENT(OUT) :: RESID
46 ! Maximum value of Residual
47  DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
48 ! IJK of Maximum value of Residual
49  INTEGER, INTENT(OUT) :: IJK_RESID
50 !-----------------------------------------------
51 ! Local variables
52 !-----------------------------------------------
53 ! Indices
54  INTEGER :: IJK, IJKW, IJKS, IJKB, IJKE, IJKN, IJKT
55 ! Numerators and denominators
56  DOUBLE PRECISION :: NUM1, DEN1
57 ! Number of fluid cells
58  INTEGER :: NCELLS
59 ! New local variables for DMP version
60  DOUBLE PRECISION, DIMENSION(ijksize3_all(myPE)) :: RESID_IJK
61  DOUBLE PRECISION :: MAX_RESID_GL(0:numpes-1), MAX_RESID_L(0:numpes-1)
62  INTEGER :: IJK_RESID_GL(0:numpes-1), IJK_RESID_L(0:numpes-1)
63  INTEGER :: nproc
64 !-----------------------------------------------
65 
66 ! initializing values
67  num = zero
68  den = zero
69  max_resid = -one
70  ncells = 0
71 
72 !!$omp parallel do private( IJK )
73  DO ijk = ijkstart3, ijkend3
74  resid_ijk(ijk) = zero
75  ENDDO
76 
77 !!$omp parallel do private( IJK, IJKW, IJKS, IJKB, IJKE, IJKN, IJKT, &
78 !!$omp& NUM1, DEN1) &
79 !!$omp& REDUCTION(+:NUM,DEN,NCELLS)
80  DO ijk = ijkstart3, ijkend3
81  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
82  IF (fluid_at(ijk)) THEN
83  ijkw = west_of(ijk)
84  ijks = south_of(ijk)
85  ijke = east_of(ijk)
86  ijkn = north_of(ijk)
87 
88 ! evaluating the residual at cell ijk:
89 ! RESp = B-sum(Anb*VARnb)-Ap*VARp
90 ! (where nb = neighbor cells and p = center/0 cell)
91  num1 = b_m(ijk,m) - (a_m(ijk,0,m)*var(ijk)+a_m(ijk,east,m)*var(ijke)+&
92  a_m(ijk,west,m)*var(ijkw)+a_m(ijk,north,m)*var(ijkn)+a_m(ijk,south,m)*var(&
93  ijks))
94 
95  IF (do_k) THEN
96  ijkb = bottom_of(ijk)
97  ijkt = top_of(ijk)
98  num1 = num1 - (a_m(ijk,top,m)*var(ijkt)+a_m(ijk,bottom,m)*var(ijkb))
99  ENDIF
100 
101  num1 = abs(num1)
102  den1 = abs(a_m(ijk,0,m)*var(ijk))
103 ! storing value of residual at each ijk location
104  resid_ijk(ijk) = num1
105 
106 ! adding to terms that are accumulated
107  ncells = ncells + 1
108  num = num + num1
109  den = den + den1
110  ENDIF
111  ENDDO
112 
113  IF(.not.debug_resid) RETURN
114 
115 ! Collecting all the information among all the procesors -
116 ! determining the global sum
117  call global_all_sum(num)
118  call global_all_sum(den)
119  call global_all_sum(ncells)
120 
121  ijk_resid = 1
122  max_resid = resid_ijk( ijk_resid )
123  DO ijk = ijkstart3, ijkend3
124  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
125  IF (resid_ijk(ijk) > max_resid) THEN
126  ijk_resid = ijk
127  max_resid = resid_ijk( ijk_resid )
128  ENDIF
129  ENDDO
130 
131 ! Determining the max residual
132  do nproc=0,numpes-1
133  if(nproc.eq.mype) then
134  max_resid_l(nproc) = max_resid
135  ijk_resid_l(nproc) = funijk_gl(i_of(ijk_resid), j_of(ijk_resid), k_of(ijk_resid))
136  else
137  max_resid_l(nproc) = 0.0
138  ijk_resid_l(nproc) = 0
139  endif
140  enddo
141 
142 ! Determining the maximum among all the procesors
143  call global_all_max(max_resid)
144 
145 ! Collecting all the information among all the procesors
146  call global_all_sum(max_resid_l, max_resid_gl)
147  call global_all_sum(ijk_resid_l, ijk_resid_gl)
148 
149 ! Calling to determine the global IJK location w.r.t. serial version
150  ijk_resid = ijkmax2
151  do nproc=0,numpes-1
152  if(max_resid_gl(nproc).eq.max_resid.and.ijk_resid_gl(nproc).lt.ijk_resid) then
153  ijk_resid = ijk_resid_gl(nproc)
154  endif
155  enddo
156 
157 ! Normalizing the residual
158  IF (den > zero) THEN
159  resid = num/den
160  max_resid = ncells*max_resid/den
161  ELSEIF (num == zero) THEN
162  resid = zero
163  max_resid = zero
164  ijk_resid = 0
165  ELSE
166  resid = undefined
167  max_resid = undefined
168 ! WRITE (LINE, *) 'Warning: All center coefficients are zero.'
169 ! CALL WRITE_ERROR ('CALC_RESID_C', LINE, 1)
170  ENDIF
171 
172  RETURN
173 
174  CONTAINS
175 
176  include 'functions.inc'
177 
178  END SUBROUTINE calc_resid_c
179 
180 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
181 ! C
182 ! Subroutine: CALC_RESID_S C
183 ! Purpose: Calculate residuals for scalar equations C
184 ! C
185 ! Author: M. Syamlal Date: 21-MAY-96 C
186 ! Reviewer: Date: C
187 ! C
188 ! Literature/Document References: C
189 ! Variables referenced: C
190 ! Variables modified: C
191 ! Local variables: C
192 ! C
193 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
194 
195  SUBROUTINE calc_resid_s(VAR, A_M, B_M, M, NUM, DEN, &
196  resid, max_resid, ijk_resid, tol)
197 
198 !-----------------------------------------------
199 ! Modules
200 !-----------------------------------------------
201  USE param
202  USE param1
203  USE parallel
204  USE geometry
205  USE indices
206  USE compar
207  USE mpi_utility
208  USE run
209 
210  USE fldvar
211  USE physprop
212  USE toleranc
213 
214  IMPLICIT NONE
215 !-----------------------------------------------
216 ! Dummy arguments
217 !-----------------------------------------------
218 ! Variable being evaluated
219  DOUBLE PRECISION, INTENT(IN) :: Var(dimension_3)
220 ! Septadiagonal matrix A_m
221  DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3, 0:dimension_m)
222 ! Vector b_m
223  DOUBLE PRECISION, INTENT(IN) :: B_m(dimension_3, 0:dimension_m)
224 ! Phase index
225  INTEGER, INTENT(IN) :: M
226 ! Numerator and denominator
227  DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
228 ! Average value of Residual
229  DOUBLE PRECISION, INTENT(OUT) :: RESID
230 ! Maximum value of Residual
231  DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
232 ! IJK of Maximum value of Residual
233  INTEGER, INTENT(OUT) :: IJK_RESID
234 ! Ignore residual calculation for scalar values below this
235  DOUBLE PRECISION, INTENT(IN) :: TOL
236 !-----------------------------------------------
237 ! Local variables
238 !-----------------------------------------------
239 ! Indices
240  INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
241 ! Numerators and denominators
242  DOUBLE PRECISION :: NUM1, DEN1
243 ! Number of fluid cells
244  INTEGER :: NCELLS
245 ! New local variables for DMP version
246  DOUBLE PRECISION, DIMENSION(ijksize3_all(myPE)) :: RESID_IJK
247  DOUBLE PRECISION :: MAX_RESID_GL(0:numpes-1), MAX_RESID_L(0:numpes-1)
248  INTEGER :: IJK_RESID_GL(0:numpes-1), IJK_RESID_L(0:numpes-1)
249  INTEGER :: nproc
250 !-----------------------------------------------
251 
252 ! initializing
253  num = zero
254  den = zero
255  max_resid = -one
256  ncells = 0
257 
258 !!$omp parallel do private( IJK )
259  DO ijk = ijkstart3, ijkend3
260  resid_ijk(ijk) = zero
261  ENDDO
262 
263 !!$omp parallel do &
264 !!$omp& private( IJK, IMJK, IJMK, IPJK, IJPK, IJKM, IJKP, &
265 !!$omp& NUM1, DEN1) &
266 !!$omp& REDUCTION(+:NUM, DEN,NCELLS)
267 
268  DO ijk = ijkstart3, ijkend3
269  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
270 
271  IF (fluid_at(ijk) .AND. abs(var(ijk)) > tol) THEN
272  imjk = im_of(ijk)
273  ijmk = jm_of(ijk)
274  ipjk = ip_of(ijk)
275  ijpk = jp_of(ijk)
276 
277  if(m/=0) then
278  if(ep_s(ijk,m) <= dil_ep_s) cycle
279  endif
280 
281 
282 ! evaluating the residual at cell ijk:
283 ! RESp = B-sum(Anb*VARnb)-Ap*VARp
284 ! (where nb = neighbor cells and p = center/0 cell)
285  num1 = b_m(ijk,m) - (a_m(ijk,0,m)*var(ijk)+a_m(ijk,east,m)*var(ipjk)+&
286  a_m(ijk,west,m)*var(imjk)+a_m(ijk,north,m)*var(ijpk)+a_m(ijk,south,m)*var(&
287  ijmk))
288  IF (do_k) THEN
289  ijkm = km_of(ijk)
290  ijkp = kp_of(ijk)
291  num1 = num1 - (a_m(ijk,top,m)*var(ijkp)+a_m(ijk,bottom,m)*var(ijkm))
292  ENDIF
293 
294  num1 = abs(num1)
295  den1 = abs(a_m(ijk,0,m)*var(ijk))
296 ! storing value of residual at each ijk location
297  resid_ijk(ijk) = num1
298 
299 ! adding to terms that are accumulated
300  ncells = ncells + 1
301  num = num + num1
302  den = den + den1
303  ENDIF
304  ENDDO
305 
306  IF(.not.debug_resid) RETURN
307 
308 ! Collecting all the information among all the procesors -
309 ! determining the global sum
310  call global_all_sum(num)
311  call global_all_sum(den)
312  call global_all_sum(ncells)
313 
314  ijk_resid = 1
315  max_resid = resid_ijk( ijk_resid )
316  DO ijk = ijkstart3, ijkend3
317  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
318  IF (resid_ijk(ijk) > max_resid) THEN
319  ijk_resid = ijk
320  max_resid = resid_ijk( ijk_resid )
321  ENDIF
322  ENDDO
323 
324 ! Determining the max residual
325  do nproc=0,numpes-1
326  if(nproc.eq.mype) then
327  max_resid_l(nproc) = max_resid
328  ijk_resid_l(nproc) = funijk_gl(i_of(ijk_resid), j_of(ijk_resid), k_of(ijk_resid))
329  else
330  max_resid_l(nproc) = 0.0
331  ijk_resid_l(nproc) = 0
332  endif
333  enddo
334 
335 ! Determining the maximum among all the procesors
336  call global_all_max(max_resid)
337 
338 ! Collecting all the information among all the procesors
339  call global_all_sum(max_resid_l, max_resid_gl)
340  call global_all_sum(ijk_resid_l, ijk_resid_gl)
341 
342 ! Determining the global IJK location w.r.t. serial version
343  ijk_resid = ijkmax2
344  do nproc=0,numpes-1
345  if(max_resid_gl(nproc).eq.max_resid.and.ijk_resid_gl(nproc).lt.ijk_resid) then
346  ijk_resid = ijk_resid_gl(nproc)
347  endif
348  enddo
349 
350 ! Normalizing the residual
351  IF (den > zero) THEN
352  resid = num/den
353  max_resid = ncells*max_resid/den
354  ELSEIF (num == zero) THEN
355  resid = zero
356  max_resid = zero
357  ijk_resid = 0
358  ELSE
359  resid = undefined
360  max_resid = undefined
361 ! WRITE(LINE,*)'Message: All center coefficients are zero.'
362 ! CALL WRITE_ERROR('CALC_RESID_S', LINE, 1)
363  ENDIF
364 
365  RETURN
366 
367  CONTAINS
368 
369  include 'functions.inc'
370 
371  END SUBROUTINE calc_resid_s
372 
373 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
374 ! C
375 ! Subroutine: CALC_RESID_pp C
376 ! Purpose: Calculate residuals for pressure correction equation C
377 ! C
378 ! Author: M. Syamlal Date: 21-JUN-96 C
379 ! Reviewer: Date: C
380 ! C
381 ! Comments: C
382 ! for correction equations the convergence for the corrections must C
383 ! go to zero, therefore the vector b must go to zero. this value C
384 ! cannot be normalized as the other equations are since the C
385 ! denominator here will vanish. thus the residual is normalized C
386 ! based on its value in the first iteration C
387 ! C
388 ! C
389 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
390 
391  SUBROUTINE calc_resid_pp(B_M, NORM, NUM, DEN, RESID, MAX_RESID, &
392  ijk_resid)
393 
394 !-----------------------------------------------
395 ! Modules
396 !-----------------------------------------------
397  USE param
398  USE param1
399  USE parallel
400  USE geometry
401  USE indices
402  USE compar
403  USE mpi_utility
404  USE run
405  IMPLICIT NONE
406 !-----------------------------------------------
407 ! Dummy arguments
408 !-----------------------------------------------
409 ! Vector b_m
410  DOUBLE PRECISION, INTENT(IN) :: B_m(dimension_3, 0:dimension_m)
411 ! Normalization factor
412  DOUBLE PRECISION, INTENT(IN) :: NORM
413 ! Numerator and denominator
414  DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
415 ! Average value of Residual
416  DOUBLE PRECISION, INTENT(OUT) :: RESID
417 ! Maximum value of Residual
418  DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
419 ! IJK of Maximum value of Residual
420  INTEGER, INTENT(OUT) :: IJK_RESID
421 !-----------------------------------------------
422 ! Local variables
423 !-----------------------------------------------
424 ! Indices
425  INTEGER :: IJK
426 ! Number of fluid cells
427  INTEGER :: NCELLS
428 ! Numerators and denominators
429  DOUBLE PRECISION :: NUM1, DEN1
430 ! New local variables for DMP version
431  DOUBLE PRECISION :: MAX_RESID_GL(0:numpes-1), MAX_RESID_L(0:numpes-1)
432  INTEGER :: IJK_RESID_GL(0:numpes-1), IJK_RESID_L(0:numpes-1)
433  INTEGER :: nproc
434 !-----------------------------------------------
435 
436 ! initializing values
437  num = zero
438  den = zero
439  max_resid = -one
440  ncells = 0
441  den1 = one
442  ijk_resid = 1
443 
444  DO ijk = ijkstart3, ijkend3
445  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
446  IF (fluid_at(ijk)) THEN
447 
448 ! evaluating the residual at cell ijk:
449  num1 = abs(b_m(ijk,0))
450  IF (num1 > max_resid) THEN
451  max_resid = num1
452  ijk_resid = ijk
453  ENDIF
454 
455 ! adding to terms that are accumulated
456  ncells = ncells + 1
457  num = num + num1
458  den = den + den1
459  ENDIF
460  ENDDO
461 
462 
463  IF(.not.debug_resid) THEN
464 ! Collecting all the information among all the procesors
465  call global_all_sum(num)
466  call global_all_sum(den)
467 
468 ! Normalizing the residual
469  IF (den*norm > zero) THEN
470 ! if norm=1 then this simply becomes an unscaled 'average' residual
471  resid = num/(den*norm)
472  ELSEIF (num == zero) THEN
473  resid = zero
474  ELSE
475  resid = large_number
476  ENDIF
477  ELSE ! if(debug_resid) branch
478 
479 ! Collecting all the information among all the procesors -
480 ! determining the global sum
481  call global_all_sum(num)
482  call global_all_sum(den)
483  call global_all_sum(ncells)
484 
485 ! Determining the max residual
486  do nproc=0,numpes-1
487  if(nproc.eq.mype) then
488  max_resid_l(nproc) = max_resid
489  ijk_resid_l(nproc) = funijk_gl(i_of(ijk_resid), j_of(ijk_resid), k_of(ijk_resid))
490  else
491  max_resid_l(nproc) = 0.0
492  ijk_resid_l(nproc) = 0
493  endif
494  enddo
495 
496 ! Determining the maximum among all the procesors
497  call global_all_max(max_resid)
498 
499 ! Collecting all the information among all the procesors
500  call global_all_sum(max_resid_l, max_resid_gl)
501  call global_all_sum(ijk_resid_l, ijk_resid_gl)
502 
503 ! Determining the global IJK location w.r.t. serial version
504  ijk_resid = ijkmax2
505  do nproc=0,numpes-1
506  if(max_resid_gl(nproc).eq.max_resid.and.ijk_resid_gl(nproc).lt.ijk_resid) then
507  ijk_resid = ijk_resid_gl(nproc)
508  endif
509  enddo
510 
511 ! Normalizing the residual
512  IF (den*norm > zero) THEN
513  resid = num/(den*norm)
514  max_resid = ncells*max_resid/(den*norm)
515  ELSEIF (num == zero) THEN
516  resid = zero
517  max_resid = zero
518  ijk_resid = 0
519  ELSE
520  resid = large_number
521  max_resid = large_number
522 ! WRITE (LINE, *) 'Warning: All center coefficients are zero.'
523 ! CALL WRITE_ERROR ('CALC_RESID_pp', LINE, 1)
524  ENDIF
525 
526  ENDIF ! end if/else debug_resid branch
527 
528  RETURN
529 
530  CONTAINS
531 
532  include 'functions.inc'
533 
534  END SUBROUTINE calc_resid_pp
535 
536 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
537 ! C
538 ! Subroutine: CALC_RESID_mb C
539 ! Purpose: Calculate overall mass balance error C
540 ! C
541 ! C
542 ! Author: M. Syamlal Date: 9-DEC-02 C
543 ! Reviewer: Date: C
544 ! C
545 ! C
546 ! C
547 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
548 
549  SUBROUTINE calc_resid_mb(INIT, ErrorPercent)
550 
551 !-----------------------------------------------
552 ! Modules
553 !-----------------------------------------------
554  USE bc
555  USE check, only: calc_mass_fluxhr, accumulation
556  USE compar
557  USE constant
558  USE fldvar
559  USE geometry
560  USE indices
561  USE mflux
562  USE mpi_utility
563  USE parallel
564  USE param
565  USE param1
566  USE physprop
567  USE residual
568  USE run, only: added_mass, dt, m_am, steady_state
569  USE rxns
570  IMPLICIT NONE
571 !-----------------------------------------------
572 ! Dummy arguments
573 !-----------------------------------------------
574 ! Flag to check whether this is an initialization call
575 ! 0 -> initialize old accumulation; 1 -> calc residual
576  INTEGER, INTENT(IN) :: init
577 ! Total mass balance error as a % of inflow
578  DOUBLE PRECISION, INTENT(OUT) :: ErrorPercent(0:mmax)
579 !-----------------------------------------------
580 ! Local variables
581 !-----------------------------------------------
582 ! phase index
583  INTEGER :: M
584 ! boundary index
585  INTEGER :: L
586 ! locally define dt, so that this routine works when dt is not defined
587  DOUBLE PRECISION :: dt_local
588 ! error index
589  INTEGER :: IER
590  DOUBLE PRECISION :: flux_in, flux_out, fin, fout
591  DOUBLE PRECISION :: err, accum_new, denom
592 !-----------------------------------------------
593 
594  if(steady_state)then
595  dt_local = one
596  else
597  dt_local = dt
598  endif
599 
600  IF(init == 0) THEN
601 ! Initialize this routine
602 ! ---------------------------------------------------------------->>>
603 
604 ! Accumulation
605  if(steady_state)then
606  accum_resid_g = zero
607  else
608  accum_resid_g = accumulation(rop_g)
609  endif
610  DO m=1, mmax
611  if(steady_state)then
612  accum_resid_s(m) = zero
613  else
614  accum_resid_s(m) = accumulation(rop_s(1,m))
615  endif
616  ENDDO
617  RETURN
618 ! end initialization
619 ! ----------------------------------------------------------------<<<
620 
621  ELSE
622 
623 ! Calculate residual
624 ! ---------------------------------------------------------------->>>
625  if(steady_state)then
626  accum_new = - accumulation(sum_r_g) * dt_local
627  else
628  accum_new = accumulation(rop_g) - accumulation(sum_r_g) * dt_local
629  endif
630 
631  flux_out = zero
632  flux_in = zero
633  DO l = 1, dimension_bc
634  IF (bc_defined(l)) THEN
635 ! call Calc_mass_flux(BC_I_W(L), BC_I_E(L), BC_J_S(L), &
636 ! BC_J_N(L), BC_K_B(L), BC_K_T(L), BC_PLANE(L), U_g, V_g, W_g, &
637 ! ROP_g, fin, fout, IER)
638  IF(.NOT.added_mass) THEN
639  call calc_mass_fluxhr(bc_i_w(l), bc_i_e(l), bc_j_s(l), &
640  bc_j_n(l), bc_k_b(l), bc_k_t(l), bc_plane(l), &
641  flux_ge, flux_gn, flux_gt, fin, fout, ier)
642  ELSE
643  call calc_mass_fluxhr(bc_i_w(l), bc_i_e(l), bc_j_s(l), &
644  bc_j_n(l), bc_k_b(l), bc_k_t(l), bc_plane(l), flux_gse,&
645  flux_gsn, flux_gst, fin, fout, ier)
646  ENDIF
647  flux_out = flux_out + fout * dt_local
648  flux_in = flux_in + fin * dt_local
649  ENDIF
650  END DO
651 
652  err = (accum_new - accum_resid_g) - (flux_in - flux_out)
653  denom = max(abs(accum_new), abs(accum_resid_g), abs(flux_in), abs(flux_out))
654  IF (denom /= zero) THEN
655  errorpercent(0) = err*100./denom
656  ELSE
657  errorpercent(0) = err*100./small_number
658  ENDIF
659 
660  DO m =1, mmax
661  if(steady_state)then
662  accum_new = - accumulation(sum_r_s(1,m)) * dt_local
663  else
664  accum_new = accumulation(rop_s(1,m)) - accumulation(sum_r_s(1,m)) * dt_local
665  endif
666 
667  flux_out = zero
668  flux_in = zero
669  DO l = 1, dimension_bc
670  IF (bc_defined(l)) THEN
671 ! call Calc_mass_flux(BC_I_W(L), BC_I_E(L), BC_J_S(L), BC_J_N(L), &
672 ! BC_K_B(L), BC_K_T(L), BC_PLANE(L), U_s(1,M), V_s(1,M), W_s(1,M), &
673 ! ROP_s(1,M), fin, fout, IER)
674  IF(.NOT.added_mass .OR. m /= m_am) THEN
675  call calc_mass_fluxhr(bc_i_w(l), bc_i_e(l), bc_j_s(l), &
676  bc_j_n(l), bc_k_b(l), bc_k_t(l), bc_plane(l), &
677  flux_se(1,m), flux_sn(1,m), flux_st(1,m), fin, fout, ier)
678  ELSE
679  call calc_mass_fluxhr(bc_i_w(l), bc_i_e(l), bc_j_s(l), &
680  bc_j_n(l), bc_k_b(l), bc_k_t(l), bc_plane(l), &
681  flux_sse, flux_ssn, flux_sst, fin, fout, ier)
682  ENDIF
683  flux_out = flux_out + fout * dt_local
684  flux_in = flux_in + fin * dt_local
685  ENDIF
686  ENDDO
687 
688  err = (accum_new - accum_resid_s(m)) - (flux_in - flux_out)
689  denom = max(abs(accum_new), abs(accum_resid_s(m)), abs(flux_in), abs(flux_out))
690  if(denom /= zero) THEN
691  errorpercent(m) = err*100./denom
692  else
693  errorpercent(m) = err*100./small_number
694  endif
695  ENDDO ! end do m=1,mmax
696 ! end calculate residual
697 ! ----------------------------------------------------------------<<<
698 
699  ENDIF ! end if/else (init==0)
700 
701  RETURN
702 
703  CONTAINS
704 
705  include 'functions.inc'
706 
707  END SUBROUTINE calc_resid_mb
708 
709 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
710 ! C
711 ! Subroutine: CALC_RESID_U C
712 ! Purpose: Calculate residuals for u-momentum equations C
713 ! C
714 ! Author: M. Syamlal Date: 21-MAY-96 C
715 ! Reviewer: Date: C
716 ! C
717 ! C
718 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
719 
720  SUBROUTINE calc_resid_u(U_M, V_M, W_M, A_M, B_M, M, NUM, DEN, &
721  resid, max_resid, ijk_resid)
722 
723 !-----------------------------------------------
724 ! Modules
725 !-----------------------------------------------
726  USE param
727  USE param1
728  USE parallel
729  USE geometry
730  USE indices
731  USE compar
732  USE mpi_utility
733  USE run
734  USE fldvar
735  USE physprop
736  USE toleranc
737  USE fun_avg
738 
739  IMPLICIT NONE
740 !-----------------------------------------------
741 ! Dummy arguments
742 !-----------------------------------------------
743 ! U velocity (x-dir)
744  DOUBLE PRECISION, INTENT(IN) :: U_m(dimension_3)
745 ! V velocity (y-dir), used here for scaling
746  DOUBLE PRECISION, INTENT(IN) :: V_m(dimension_3)
747 ! W velocity (z-dir), used here for scaling
748  DOUBLE PRECISION, INTENT(IN) :: W_m(dimension_3)
749 ! Septadiagonal matrix A_m
750  DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3, 0:dimension_m)
751 ! Vector b_m
752  DOUBLE PRECISION, INTENT(IN) :: B_m(dimension_3, 0:dimension_m)
753 ! Phase index
754  INTEGER, INTENT(IN) :: M
755 ! Numerator and denominator
756  DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
757 ! Average value of Residual
758  DOUBLE PRECISION, INTENT(OUT) :: RESID
759 ! Maximum value of Residual
760  DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
761 ! IJK of Maximum value of Residual
762  INTEGER, INTENT(OUT) :: IJK_RESID
763 !-----------------------------------------------
764 ! Local variables
765 !-----------------------------------------------
766 ! Velocity magnitude
767  DOUBLE PRECISION :: VEL
768 ! Indices
769  INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
770 ! Numerators and denominators
771  DOUBLE PRECISION :: NUM1, DEN1
772 ! Number of fluid cells
773  INTEGER :: NCELLS
774 ! New local variables for DMP version
775  DOUBLE PRECISION, DIMENSION(ijksize3_all(myPE)) :: RESID_IJK
776  DOUBLE PRECISION :: MAX_RESID_GL(0:numpes-1), MAX_RESID_L(0:numpes-1)
777  INTEGER :: IJK_RESID_GL(0:numpes-1), IJK_RESID_L(0:numpes-1)
778  INTEGER :: nproc
779 
780 ! Solids volume fraction at face
781  DOUBLE PRECISION :: EPSA
782 
783 !-----------------------------------------------
784 
785 ! initializing
786  num = zero
787  den = zero
788  max_resid = -one
789  ncells = 0
790 
791 !$omp parallel default(none) &
792 !$omp private( IMJK, IPJK, IJMK, IJPK, IJKM, IJKP, &
793 !$omp VEL, NUM1, DEN1,EPSA) &
794 !$omp shared (ijkstart3, ijkend3, resid_ijk, i_of, j_of, k_of,m,a_m,b_m,w_m,do_k,u_m,v_m,num,den,ncells)
795 !$omp do reduction(+:num, DEN, NCELLS )
796  DO ijk = ijkstart3, ijkend3
797  resid_ijk(ijk) = zero
798 
799  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
800 
801 ! Skip walls where some values are undefined.
802  IF(wall_at(ijk)) cycle
803 
804  if(m/=0) then
805  epsa = avg_x(ep_s(ijk,m),ep_s(east_of(ijk),m),i_of(ijk))
806  if(epsa <= dil_ep_s) cycle
807  endif
808 
809  IF (.NOT.ip_at_e(ijk)) THEN
810  imjk = im_of(ijk)
811  ijmk = jm_of(ijk)
812  ipjk = ip_of(ijk)
813  ijpk = jp_of(ijk)
814 
815 ! evaluating the residual at cell ijk:
816 ! RESp = B-sum(Anb*VARnb)-Ap*VARp
817 ! (where nb = neighbor cells and p = center/0 cell)
818  num1 = b_m(ijk,m) - (a_m(ijk,0,m)*u_m(ijk)+&
819  a_m(ijk,east,m)*u_m(ipjk)+a_m(ijk,west,m)*u_m(imjk)+&
820  a_m(ijk,north,m)*u_m(ijpk)+a_m(ijk,south,m)*u_m(ijmk))
821  IF (do_k) THEN
822  ijkm = km_of(ijk)
823  ijkp = kp_of(ijk)
824  num1 = num1 - (a_m(ijk,top,m)*u_m(ijkp)+a_m(ijk,bottom,m)*u_m(ijkm))
825  ENDIF
826 
827 ! Ignore momentum residual in stagnant regions. Need an alternative
828 ! criteria for residual scaling for such cases.
829  vel = sqrt(u_m(ijk)**2+v_m(ijk)**2+w_m(ijk)**2)
830  IF (vel > small_number) THEN
831  num1 = abs(num1)
832  den1 = abs(a_m(ijk,0,m)*vel)
833 ! storing value of residual at each ijk location
834  resid_ijk(ijk) = num1
835 ! adding to terms that are accumulated
836  ncells = ncells + 1
837  num = num + num1
838  den = den + den1
839  ENDIF
840  ENDIF
841  ENDDO
842 !$omp end parallel
843 
844  IF(.not.debug_resid) RETURN
845 
846 ! Collecting all the information among all the procesors -
847 ! determining the global sum
848  call global_all_sum(num)
849  call global_all_sum(den)
850  call global_all_sum(ncells)
851 
852  ijk_resid = 1
853  max_resid = resid_ijk( ijk_resid )
854  DO ijk = ijkstart3, ijkend3
855  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
856  IF (resid_ijk( ijk ) > max_resid) THEN
857  ijk_resid = ijk
858  max_resid = resid_ijk( ijk_resid )
859  ENDIF
860  ENDDO
861 
862 ! Determining the max residual
863  do nproc=0,numpes-1
864  if(nproc.eq.mype) then
865  max_resid_l(nproc) = max_resid
866  ijk_resid_l(nproc) = funijk_gl(i_of(ijk_resid), j_of(ijk_resid), k_of(ijk_resid))
867  else
868  max_resid_l(nproc) = 0.0
869  ijk_resid_l(nproc) = 0
870  endif
871  enddo
872 
873 ! Determining the maximum among all the procesors
874  call global_all_max(max_resid)
875 
876 ! Collecting all the information among all the procesors
877  call global_all_sum(max_resid_l, max_resid_gl)
878  call global_all_sum(ijk_resid_l, ijk_resid_gl)
879 
880 ! Determining the global IJK location w.r.t. serial version
881  ijk_resid = ijkmax2
882  do nproc=0,numpes-1
883  IF(max_resid_gl(nproc).eq.max_resid.and.&
884  ijk_resid_gl(nproc).lt.ijk_resid) THEN
885  ijk_resid = ijk_resid_gl(nproc)
886  ENDIF
887  ENDDO
888 
889 ! Normalizing the residual
890  IF (den > zero) THEN
891  resid = num/den
892  max_resid = ncells*max_resid/den
893  ELSEIF (num == zero) THEN
894  resid = zero
895  max_resid = zero
896  ijk_resid = 0
897  ELSE
898  resid = large_number
899  max_resid = large_number
900 ! WRITE (LINE, *) 'Warning: All center coefficients are zero.'
901 ! CALL WRITE_ERROR ('CALC_RESID_U', LINE, 1)
902  ENDIF
903 
904  RETURN
905 
906  CONTAINS
907 
908  include 'functions.inc'
909 
910  END SUBROUTINE calc_resid_u
911 
912 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
913 ! C
914 ! Subroutine: CALC_RESID_V C
915 ! Purpose: Calculate residuals for v-momentum equations C
916 ! C
917 ! Author: M. Syamlal Date: 21-MAY-96 C
918 ! Reviewer: Date: C
919 ! C
920 ! C
921 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
922 
923  SUBROUTINE calc_resid_v(U_M, V_M, W_M, A_M, B_M, M, NUM, DEN, &
924  resid, max_resid, ijk_resid)
925 
926 !-----------------------------------------------
927 ! Modules
928 !-----------------------------------------------
929  USE param
930  USE param1
931  USE parallel
932  USE geometry
933  USE indices
934  USE compar
935  USE mpi_utility
936  USE run
937  USE fldvar
938  USE physprop
939  USE toleranc
940  USE fun_avg
941 
942  IMPLICIT NONE
943 !-----------------------------------------------
944 ! Dummy arguments
945 !-----------------------------------------------
946 ! U velocity (x-dir)
947  DOUBLE PRECISION, INTENT(IN) :: U_m(dimension_3)
948 ! V velocity (y-dir), used here for scaling
949  DOUBLE PRECISION, INTENT(IN) :: V_m(dimension_3)
950 ! W velocity (z-dir), used here for scaling
951  DOUBLE PRECISION, INTENT(IN) :: W_m(dimension_3)
952 ! Septadiagonal matrix A_m
953  DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3, 0:dimension_m)
954 ! Vector b_m
955  DOUBLE PRECISION, INTENT(IN) :: B_m(dimension_3, 0:dimension_m)
956 ! Phase index
957  INTEGER, INTENT(IN) :: M
958 ! Numerator and denominator
959  DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
960 ! Average value of Residual
961  DOUBLE PRECISION, INTENT(OUT) :: RESID
962 ! Maximum value of Residual
963  DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
964 ! IJK of Maximum value of Residual
965  INTEGER, INTENT(OUT) :: IJK_RESID
966 !-----------------------------------------------
967 ! Local variables
968 !-----------------------------------------------
969 ! Velocity magnitude
970  DOUBLE PRECISION :: VEL
971 ! Indices
972  INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
973 ! Numerators and denominators
974  DOUBLE PRECISION :: NUM1, DEN1
975 ! Number of fluid cells
976  INTEGER :: NCELLS
977 ! New local variables for DMP version
978  DOUBLE PRECISION, DIMENSION(ijksize3_all(myPE)) :: RESID_IJK
979  DOUBLE PRECISION :: MAX_RESID_GL(0:numpes-1), MAX_RESID_L(0:numpes-1)
980  INTEGER :: IJK_RESID_GL(0:numpes-1), IJK_RESID_L(0:numpes-1)
981  INTEGER :: nproc
982 ! Solids volume fraction at face
983  DOUBLE PRECISION :: EPSA
984 
985 !-----------------------------------------------
986 
987 ! initializing
988  num = zero
989  den = zero
990  max_resid = -one
991  ncells = 0
992 
993 !$omp parallel default(none) &
994 !$omp private( IMJK, IPJK, IJMK, IJPK, IJKM, IJKP, &
995 !$omp VEL, NUM1, DEN1,EPSA) &
996 !$omp shared (ijkstart3, ijkend3, resid_ijk, i_of, j_of, k_of,m,a_m,b_m,w_m,do_k,u_m,v_m,num,den,ncells)
997 !$omp do reduction(+:num, DEN, NCELLS )
998  DO ijk = ijkstart3, ijkend3
999  resid_ijk(ijk) = zero
1000 
1001  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
1002 
1003 ! Skip walls where some values are undefined.
1004  IF(wall_at(ijk)) cycle
1005 
1006  if(m/=0) then
1007  epsa = avg_y(ep_s(ijk,m),ep_s(north_of(ijk),m),j_of(ijk))
1008  if(epsa <= dil_ep_s) cycle
1009  endif
1010 
1011  IF (.NOT.ip_at_n(ijk)) THEN
1012  imjk = im_of(ijk)
1013  ijmk = jm_of(ijk)
1014  ipjk = ip_of(ijk)
1015  ijpk = jp_of(ijk)
1016 
1017 ! evaluating the residual at cell ijk:
1018 ! RESp = B-sum(Anb*VARnb)-Ap*VARp
1019 ! (where nb = neighbor cells and p = center/0 cell)
1020  num1 = b_m(ijk,m) - (a_m(ijk,0,m)*v_m(ijk)+&
1021  a_m(ijk,east,m)*v_m(ipjk)+a_m(ijk,west,m)*v_m(imjk)+&
1022  a_m(ijk,north,m)*v_m(ijpk)+a_m(ijk,south,m)*v_m(ijmk))
1023  IF (do_k) THEN
1024  ijkm = km_of(ijk)
1025  ijkp = kp_of(ijk)
1026  num1 = num1 - (a_m(ijk,top,m)*v_m(ijkp)+a_m(ijk,bottom,m)*v_m(ijkm))
1027  ENDIF
1028 
1029 ! Ignore momentum residual in stagnant regions. Need an alternative
1030 ! criteria for residual scaling for such cases.
1031  vel = sqrt(u_m(ijk)**2+v_m(ijk)**2+w_m(ijk)**2)
1032  IF (vel > small_number) THEN
1033  num1 = abs(num1)
1034  den1 = abs(a_m(ijk,0,m)*vel)
1035 ! storing value of residual at each ijk location
1036  resid_ijk(ijk) = num1
1037 ! adding to terms that are accumulated
1038  ncells = ncells + 1
1039  num = num + num1
1040  den = den + den1
1041  ENDIF
1042  ENDIF
1043  ENDDO
1044 !$omp end parallel
1045 
1046  if(.not.debug_resid) return
1047 
1048 ! Collecting all the information among all the procesors -
1049 ! determining the global sum
1050  call global_all_sum(num)
1051  call global_all_sum(den)
1052  call global_all_sum(ncells)
1053 
1054  ijk_resid = 1
1055  max_resid = resid_ijk( ijk_resid )
1056  DO ijk = ijkstart3, ijkend3
1057  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
1058  IF (resid_ijk( ijk ) > max_resid) THEN
1059  ijk_resid = ijk
1060  max_resid = resid_ijk( ijk_resid )
1061  ENDIF
1062  ENDDO
1063 
1064 ! Determining the max residual
1065  do nproc=0,numpes-1
1066  if(nproc.eq.mype) then
1067  max_resid_l(nproc) = max_resid
1068  ijk_resid_l(nproc) = funijk_gl(i_of(ijk_resid), j_of(ijk_resid), k_of(ijk_resid))
1069  else
1070  max_resid_l(nproc) = 0.0
1071  ijk_resid_l(nproc) = 0
1072  endif
1073  enddo
1074 
1075 ! Determine the maximum among all the procesors
1076  call global_all_max(max_resid)
1077 
1078 ! Collecting all the information among all the procesors
1079  call global_all_sum(max_resid_l, max_resid_gl)
1080  call global_all_sum(ijk_resid_l, ijk_resid_gl)
1081 
1082 ! Determining the global IJK location w.r.t. serial version
1083  ijk_resid = ijkmax2
1084  do nproc=0,numpes-1
1085  if(max_resid_gl(nproc).eq.max_resid.and.ijk_resid_gl(nproc).lt.ijk_resid) then
1086  ijk_resid = ijk_resid_gl(nproc)
1087  endif
1088  enddo
1089 
1090 ! Normalizing the residual
1091  IF (den > zero) THEN
1092  resid = num/den
1093  max_resid = ncells*max_resid/den
1094  ELSEIF (num == zero) THEN
1095  resid = zero
1096  max_resid = zero
1097  ijk_resid = 0
1098  ELSE
1099  resid = large_number
1100  max_resid = large_number
1101 ! WRITE (LINE, *) 'Warning: All center coefficients are zero.'
1102 ! CALL WRITE_ERROR ('CALC_RESID_V', LINE, 1)
1103  ENDIF
1104 
1105  RETURN
1106 
1107  CONTAINS
1108 
1109  include 'functions.inc'
1110 
1111  END SUBROUTINE calc_resid_v
1112 
1113 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1114 ! C
1115 ! Subroutine: CALC_RESID_W C
1116 ! Purpose: Calculate residuals for w-momentum equations C
1117 ! C
1118 ! Author: M. Syamlal Date: 21-MAY-96 C
1119 ! Reviewer: Date: C
1120 ! C
1121 ! Literature/Document References: C
1122 ! Variables referenced: C
1123 ! Variables modified: C
1124 ! Local variables: C
1125 ! C
1126 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1127 
1128  SUBROUTINE calc_resid_w(U_M, V_M, W_M, A_M, B_M, M, NUM, DEN, &
1129  resid, max_resid, ijk_resid)
1130 
1131 !-----------------------------------------------
1132 ! Modules
1133 !-----------------------------------------------
1134  USE param
1135  USE param1
1136  USE parallel
1137  USE geometry
1138  USE indices
1139  USE compar
1140  USE mpi_utility
1141  USE run
1142  USE fldvar
1143  USE physprop
1144  USE toleranc
1145  USE fun_avg
1146 
1147  IMPLICIT NONE
1148 !-----------------------------------------------
1149 ! Dummy arguments
1150 !-----------------------------------------------
1151 ! U velocity (x-dir)
1152  DOUBLE PRECISION, INTENT(IN) :: U_m(dimension_3)
1153 ! V velocity (y-dir), used here for scaling
1154  DOUBLE PRECISION, INTENT(IN) :: V_m(dimension_3)
1155 ! W velocity (z-dir), used here for scaling
1156  DOUBLE PRECISION, INTENT(IN) :: W_m(dimension_3)
1157 ! Septadiagonal matrix A_m
1158  DOUBLE PRECISION, INTENT(IN) :: A_m(dimension_3, -3:3, 0:dimension_m)
1159 ! Vector b_m
1160  DOUBLE PRECISION, INTENT(IN) :: B_m(dimension_3, 0:dimension_m)
1161 ! Phase index
1162  INTEGER, INTENT(IN) :: M
1163 ! Numerator and denominator
1164  DOUBLE PRECISION, INTENT(OUT) :: NUM, DEN
1165 ! Average value of Residual
1166  DOUBLE PRECISION, INTENT(OUT) :: RESID
1167 ! Maximum value of Residual
1168  DOUBLE PRECISION, INTENT(OUT) :: MAX_RESID
1169 ! IJK of Maximum value of Residual
1170  INTEGER, INTENT(OUT) :: IJK_RESID
1171 !-----------------------------------------------
1172 ! Local variables
1173 !-----------------------------------------------
1174 ! Velocity magnitude
1175  DOUBLE PRECISION :: VEL
1176 ! Indices
1177  INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
1178 ! Numerators and denominators
1179  DOUBLE PRECISION :: NUM1, DEN1
1180 ! Number of fluid cells
1181  INTEGER :: NCELLS
1182 ! New local variables for DMP version
1183  DOUBLE PRECISION, DIMENSION(ijksize3_all(myPE)) :: RESID_IJK
1184  DOUBLE PRECISION :: MAX_RESID_GL(0:numpes-1), MAX_RESID_L(0:numpes-1)
1185  INTEGER :: IJK_RESID_GL(0:numpes-1), IJK_RESID_L(0:numpes-1)
1186  INTEGER :: nproc
1187 ! Solids volume fraction at face
1188  DOUBLE PRECISION :: EPSA
1189 
1190 !-----------------------------------------------
1191 
1192 ! initializing
1193  num = zero
1194  den = zero
1195  max_resid = -one
1196  ncells = 0
1197 
1198 !$omp parallel default(none) &
1199 !$omp private( IMJK, IPJK, IJMK, IJPK, IJKM, IJKP, &
1200 !$omp VEL, NUM1, DEN1,EPSA) &
1201 !$omp shared (ijkstart3, ijkend3, resid_ijk, i_of, j_of, k_of,m,a_m,b_m,w_m,do_k,u_m,v_m,num,den,ncells)
1202 !$omp do reduction(+:num, DEN, NCELLS )
1203  DO ijk = ijkstart3, ijkend3
1204  resid_ijk(ijk) = zero
1205 
1206  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
1207 
1208 ! Skip walls where some values are undefined.
1209  IF(wall_at(ijk)) cycle
1210 
1211  if(m/=0) then
1212  epsa = avg_z(ep_s(ijk,m),ep_s(top_of(ijk),m),k_of(ijk))
1213  if(epsa <= dil_ep_s) cycle
1214  endif
1215 
1216  IF (.NOT.ip_at_t(ijk)) THEN
1217  imjk = im_of(ijk)
1218  ijmk = jm_of(ijk)
1219  ipjk = ip_of(ijk)
1220  ijpk = jp_of(ijk)
1221 
1222 ! evaluating the residual at cell ijk:
1223 ! RESp = B-sum(Anb*VARnb)-Ap*VARp
1224 ! (where nb = neighbor cells and p = center/0 cell)
1225  num1 = b_m(ijk,m) - (a_m(ijk,0,m)*w_m(ijk)+&
1226  a_m(ijk,east,m)*w_m(ipjk)+a_m(ijk,west,m)*w_m(imjk)+&
1227  a_m(ijk,north,m)*w_m(ijpk)+a_m(ijk,south,m)*w_m(ijmk))
1228  IF (do_k) THEN
1229  ijkm = km_of(ijk)
1230  ijkp = kp_of(ijk)
1231  num1 = num1 - (a_m(ijk,top,m)*w_m(ijkp)+a_m(ijk,bottom,m)*w_m(ijkm))
1232  ENDIF
1233 
1234 ! Ignore momentum residual in stagnant regions. Need an alternative
1235 ! criteria for residual scaling for such cases.
1236  vel = sqrt(u_m(ijk)**2+v_m(ijk)**2+w_m(ijk)**2)
1237  IF (vel > small_number) THEN
1238  num1 = abs(num1)
1239  den1 = abs(a_m(ijk,0,m)*vel)
1240 ! storing value of residual at each ijk location
1241  resid_ijk(ijk) = num1
1242 ! adding to terms that are accumulated
1243  ncells = ncells + 1
1244  num = num + num1
1245  den = den + den1
1246  ENDIF
1247  ENDIF
1248  ENDDO
1249 !$omp end parallel
1250 
1251  if(.not.debug_resid) return
1252 
1253 ! Collecting all the information among all the procesors -
1254 ! determining the global sum
1255  call global_all_sum(num)
1256  call global_all_sum(den)
1257  call global_all_sum(ncells)
1258 
1259  ijk_resid = 1
1260  max_resid = resid_ijk( ijk_resid )
1261  DO ijk = ijkstart3, ijkend3
1262  IF(.NOT.is_on_mype_wobnd(i_of(ijk),j_of(ijk), k_of(ijk))) cycle
1263 
1264  IF (resid_ijk( ijk ) > max_resid) THEN
1265  ijk_resid = ijk
1266  max_resid = resid_ijk( ijk_resid )
1267  ENDIF
1268  ENDDO
1269 
1270 ! Determining the max residual
1271  do nproc=0,numpes-1
1272  if(nproc.eq.mype) then
1273  max_resid_l(nproc) = max_resid
1274  ijk_resid_l(nproc) = funijk_gl(i_of(ijk_resid), j_of(ijk_resid), k_of(ijk_resid))
1275  else
1276  max_resid_l(nproc) = 0.0
1277  ijk_resid_l(nproc) = 0
1278  endif
1279  enddo
1280 
1281 ! Determining the maximum among all the procesors
1282  call global_all_max(max_resid)
1283 
1284 ! Collecting all the information among all the procesors
1285  call global_all_sum(max_resid_l, max_resid_gl)
1286  call global_all_sum(ijk_resid_l, ijk_resid_gl)
1287 
1288 ! Determining the global IJK location w.r.t. serial version
1289  ijk_resid = ijkmax2
1290 
1291  do nproc=0,numpes-1
1292  if(max_resid_gl(nproc).eq.max_resid.and.ijk_resid_gl(nproc).lt.ijk_resid) then
1293  ijk_resid = ijk_resid_gl(nproc)
1294  endif
1295  enddo
1296 
1297 ! Normalizing the residual
1298  IF (den > zero) THEN
1299  resid = num/den
1300  max_resid = ncells*max_resid/den
1301  ELSE IF (num == zero) THEN
1302  resid = zero
1303  max_resid = zero
1304  ijk_resid = 0
1305  ELSE
1306  resid = large_number
1307  max_resid = large_number
1308 ! WRITE (LINE, *) 'Warning: All center coefficients are zero.'
1309 ! CALL WRITE_ERROR ('CALC_RESID_W', LINE, 1)
1310  ENDIF
1311 
1312  RETURN
1313 
1314  CONTAINS
1315 
1316  include 'functions.inc'
1317 
1318  END SUBROUTINE calc_resid_w
calc_resid_u
subroutine calc_resid_u(U_M, V_M, W_M, A_M, B_M, M, NUM, DEN, RESID, MAX_RESID, IJK_RESID)
Definition: calc_resid.f:722
bc::bc_k_b
integer, dimension(dimension_bc) bc_k_b
Definition: bc_mod.f:70
mpi_utility
Definition: mpi_utility_mod.f:6
mflux::flux_ge
double precision, dimension(:), allocatable flux_ge
Definition: mflux_mod.f:13
calc_resid_mb
subroutine calc_resid_mb(INIT, ErrorPercent)
Definition: calc_resid.f:550
run::steady_state
logical steady_state
Definition: run_mod.f:57
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
residual::accum_resid_s
double precision, dimension(dim_m) accum_resid_s
Definition: residual_mod.f:70
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
mflux::flux_gst
double precision, dimension(:), allocatable flux_gst
Definition: mflux_mod.f:32
mflux::flux_st
double precision, dimension(:,:), allocatable flux_st
Definition: mflux_mod.f:24
param1::one
double precision, parameter one
Definition: param1_mod.f:29
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
geometry::ijkmax2
integer ijkmax2
Definition: geometry_mod.f:80
bc::bc_i_w
integer, dimension(dimension_bc) bc_i_w
Definition: bc_mod.f:54
run::added_mass
logical added_mass
Definition: run_mod.f:91
bc::bc_j_n
integer, dimension(dimension_bc) bc_j_n
Definition: bc_mod.f:66
rxns
Definition: rxns_mod.f:1
mflux::flux_ssn
double precision, dimension(:), allocatable flux_ssn
Definition: mflux_mod.f:31
check::calc_mass_fluxhr
subroutine calc_mass_fluxhr(I1, I2, J1, J2, K1, K2, Plane, Flux_E
Definition: check_mod.f:660
rxns::sum_r_g
double precision, dimension(:), allocatable sum_r_g
Definition: rxns_mod.f:28
rxns::sum_r_s
double precision, dimension(:,:), allocatable sum_r_s
Definition: rxns_mod.f:35
run::dt
double precision dt
Definition: run_mod.f:51
param::dimension_bc
integer, parameter dimension_bc
Definition: param_mod.f:61
calc_resid_s
subroutine calc_resid_s(VAR, A_M, B_M, M, NUM, DEN, RESID, MAX_RESID, IJK_RESID, TOL)
Definition: calc_resid.f:197
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
param::east
integer east
Definition: param_mod.f:29
compar::numpes
integer numpes
Definition: compar_mod.f:24
bc::bc_plane
character, dimension(dimension_bc) bc_plane
Definition: bc_mod.f:217
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
mflux::flux_gse
double precision, dimension(:), allocatable flux_gse
Definition: mflux_mod.f:28
fun_avg
Definition: fun_avg_mod.f:1
bc::bc_k_t
integer, dimension(dimension_bc) bc_k_t
Definition: bc_mod.f:74
check::accumulation
double precision function accumulation(ro)
Definition: check_mod.f:829
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
mpi_utility::global_all_max
Definition: mpi_utility_mod.f:110
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
mflux::flux_gn
double precision, dimension(:), allocatable flux_gn
Definition: mflux_mod.f:15
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
bc::bc_j_s
integer, dimension(dimension_bc) bc_j_s
Definition: bc_mod.f:62
run::m_am
integer m_am
Definition: run_mod.f:94
param::north
integer north
Definition: param_mod.f:37
bc::bc_defined
logical, dimension(dimension_bc) bc_defined
Definition: bc_mod.f:207
calc_resid_pp
subroutine calc_resid_pp(B_M, NORM, NUM, DEN, RESID, MAX_RESID, IJK_RESID)
Definition: calc_resid.f:393
param::south
integer south
Definition: param_mod.f:41
mpi_utility::global_all_sum
Definition: mpi_utility_mod.f:66
run
Definition: run_mod.f:13
param1::large_number
double precision, parameter large_number
Definition: param1_mod.f:23
param
Definition: param_mod.f:2
toleranc::dil_ep_s
double precision, parameter dil_ep_s
Definition: toleranc_mod.f:24
mflux::flux_gt
double precision, dimension(:), allocatable flux_gt
Definition: mflux_mod.f:17
residual::accum_resid_g
double precision accum_resid_g
Definition: residual_mod.f:70
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
calc_resid_c
subroutine calc_resid_c(VAR, A_M, B_M, M, NUM, DEN, RESID, MAX_RESID, IJK_RESID)
Definition: calc_resid.f:18
compar::mype
integer mype
Definition: compar_mod.f:24
calc_resid_w
subroutine calc_resid_w(U_M, V_M, W_M, A_M, B_M, M, NUM, DEN, RESID, MAX_RESID, IJK_RESID)
Definition: calc_resid.f:1130
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
param::west
integer west
Definition: param_mod.f:33
calc_resid_v
subroutine calc_resid_v(U_M, V_M, W_M, A_M, B_M, M, NUM, DEN, RESID, MAX_RESID, IJK_RESID)
Definition: calc_resid.f:925
fldvar::ep_s
double precision function ep_s(IJK, xxM)
Definition: fldvar_mod.f:178
param::top
integer top
Definition: param_mod.f:45
fldvar::rop_s
double precision, dimension(:,:), allocatable rop_s
Definition: fldvar_mod.f:51
mflux::flux_se
double precision, dimension(:,:), allocatable flux_se
Definition: mflux_mod.f:22
mflux::flux_gsn
double precision, dimension(:), allocatable flux_gsn
Definition: mflux_mod.f:30
param::dimension_m
integer dimension_m
Definition: param_mod.f:18
fldvar::rop_g
double precision, dimension(:), allocatable rop_g
Definition: fldvar_mod.f:38
mflux::flux_sn
double precision, dimension(:,:), allocatable flux_sn
Definition: mflux_mod.f:20
param::bottom
integer bottom
Definition: param_mod.f:49
bc::bc_i_e
integer, dimension(dimension_bc) bc_i_e
Definition: bc_mod.f:58
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
mflux::flux_sst
double precision, dimension(:), allocatable flux_sst
Definition: mflux_mod.f:33
bc
Definition: bc_mod.f:23
mflux::flux_sse
double precision, dimension(:), allocatable flux_sse
Definition: mflux_mod.f:29
run::debug_resid
logical debug_resid
Definition: run_mod.f:244