File: RELATIVE:/../../../mfix.git/model/calc_resid.f

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