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