File: N:\mfix\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
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
1319