File: /nfs/home/0/users/jenkins/mfix.git/model/conv_dif_phi.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29 SUBROUTINE CONV_DIF_PHI(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M,IER)
30
31
32
33
34
35
36
37
38 USE param
39 USE param1
40 USE run
41 USE geometry
42 USE compar
43 USE sendrecv
44 Use xsi_array
45 USE mpi_utility
46 USE indices
47 IMPLICIT NONE
48
49
50
51
52
53
54
55
56
57 DOUBLE PRECISION Phi(DIMENSION_3)
58
59
60 DOUBLE PRECISION Dif(DIMENSION_3)
61
62
63 INTEGER Disc
64
65
66 DOUBLE PRECISION Uf(DIMENSION_3), Vf(DIMENSION_3), Wf(DIMENSION_3)
67
68
69 DOUBLE PRECISION Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
70
71
72 INTEGER M
73
74
75 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
76
77
78 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
79
80
81 INTEGER IER
82
83
84
85
86
87
88
89
90 IF(DEF_COR)THEN
91 CALL CONV_DIF_PHI0(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M)
92 if (DISC > 1) CALL CONV_DIF_PHI_DC(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M,IER)
93 ELSE
94
95
96
97 IF (DISC == 0) THEN
98 CALL CONV_DIF_PHI0(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M)
99 ELSE
100 CALL CONV_DIF_PHI1(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M,IER)
101 ENDIF
102 ENDIF
103
104 CALL DIF_PHI_IS (DIF, A_M, B_M, M)
105
106 RETURN
107 END SUBROUTINE CONV_DIF_PHI
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132 SUBROUTINE CONV_DIF_PHI0(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M)
133
134
135
136
137
138
139
140
141 USE param
142 USE param1
143 USE parallel
144 USE matrix
145 USE toleranc
146 USE run
147 USE geometry
148 USE compar
149 USE sendrecv
150 USE indices
151 USE fun_avg
152 USE functions
153
154
155
156 USE cutcell
157
158
159
160 IMPLICIT NONE
161
162
163
164
165
166
167
168
169 DOUBLE PRECISION Phi(DIMENSION_3)
170
171
172 DOUBLE PRECISION Dif(DIMENSION_3)
173
174
175 INTEGER Disc
176
177
178 DOUBLE PRECISION Uf(DIMENSION_3), Vf(DIMENSION_3), Wf(DIMENSION_3)
179
180
181 DOUBLE PRECISION Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
182
183
184 INTEGER M
185
186
187 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
188
189
190 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
191
192
193 INTEGER I, J, K, IJK, IPJK, IJPK, IJKE, IJKN,&
194 IJKP, IJKT
195
196 INTEGER IMJK, IM, IJKW
197 INTEGER IJMK, JM, IJKS
198 INTEGER IJKM, KM, IJKB
199
200
201 DOUBLE PRECISION V_f
202
203
204 DOUBLE PRECISION D_f
205
206
207
208
209
210
211
212
213
214
215
216
217 DO IJK = ijkstart3, ijkend3
218
219 = I_OF(IJK)
220 J = J_OF(IJK)
221 K = K_OF(IJK)
222
223 IF (FLUID_AT(IJK)) THEN
224
225 = IP_OF(IJK)
226 IJPK = JP_OF(IJK)
227 IJKE = EAST_OF(IJK)
228 IJKN = NORTH_OF(IJK)
229
230
231
232 = UF(IJK)
233 D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*AYZ(IJK)
234
235
236
237 IF(CUT_TREATMENT_AT(IJK)) THEN
238 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IPJK))) THEN
239 D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*DY(J)*DZ(K)
240 ENDIF
241 ENDIF
242
243
244
245
246 IF (V_F >= ZERO) THEN
247 A_M(IJK,E,M) = D_F
248 A_M(IPJK,W,M) = D_F + FLUX_E(IJK)
249 ELSE
250 A_M(IJK,E,M) = D_F - FLUX_E(IJK)
251 A_M(IPJK,W,M) = D_F
252 ENDIF
253
254
255
256 = VF(IJK)
257 D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*AXZ(IJK)
258
259
260
261 IF(CUT_TREATMENT_AT(IJK)) THEN
262 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJPK))) THEN
263 D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*DX(I)*DZ(K)
264 ENDIF
265 ENDIF
266
267
268
269 IF (V_F >= ZERO) THEN
270 A_M(IJK,N,M) = D_F
271 A_M(IJPK,S,M) = D_F + FLUX_N(IJK)
272 ELSE
273 A_M(IJK,N,M) = D_F - FLUX_N(IJK)
274 A_M(IJPK,S,M) = D_F
275 ENDIF
276
277
278 IF (DO_K) THEN
279 IJKP = KP_OF(IJK)
280 IJKT = TOP_OF(IJK)
281 V_F = WF(IJK)
282 D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*AXY(IJK)
283
284
285
286 IF(CUT_TREATMENT_AT(IJK)) THEN
287 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKP))) THEN
288 D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*DX(I)*DY(J)
289 ENDIF
290 ENDIF
291
292
293
294 IF (V_F >= ZERO) THEN
295 A_M(IJK,T,M) = D_F
296 A_M(IJKP,B,M) = D_F + FLUX_T(IJK)
297 ELSE
298 A_M(IJK,T,M) = D_F - FLUX_T(IJK)
299 A_M(IJKP,B,M) = D_F
300 ENDIF
301 ENDIF
302
303
304
305 = IM_OF(IJK)
306 IF (.NOT.FLUID_AT(IMJK)) THEN
307 IM = IM1(I)
308 IJKW = WEST_OF(IJK)
309 V_F = UF(IMJK)
310 D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*AYZ(IMJK)
311
312
313
314 IF(CUT_TREATMENT_AT(IJK)) THEN
315 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IMJK))) THEN
316 D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*DY(J)*DZ(K)
317 ENDIF
318 ENDIF
319
320
321
322 IF (V_F >= ZERO) THEN
323 A_M(IJK,W,M) = D_F + FLUX_E(IMJK)
324 ELSE
325 A_M(IJK,W,M) = D_F
326 ENDIF
327 ENDIF
328
329
330 = JM_OF(IJK)
331 IF (.NOT.FLUID_AT(IJMK)) THEN
332 JM = JM1(J)
333 IJKS = SOUTH_OF(IJK)
334 V_F = VF(IJMK)
335 D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*AXZ(IJMK)
336
337
338
339 IF(CUT_TREATMENT_AT(IJK)) THEN
340 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJMK))) THEN
341 D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*DX(I)*DZ(K)
342 ENDIF
343 ENDIF
344
345
346
347 IF (V_F >= ZERO) THEN
348 A_M(IJK,S,M) = D_F + FLUX_N(IJMK)
349 ELSE
350 A_M(IJK,S,M) = D_F
351 ENDIF
352 ENDIF
353
354
355 IF (DO_K) THEN
356 IJKM = KM_OF(IJK)
357 IF (.NOT.FLUID_AT(IJKM)) THEN
358 KM = KM1(K)
359 IJKB = BOTTOM_OF(IJK)
360 V_F = WF(IJKM)
361 D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*AXY(&
362 IJKM)
363
364
365
366 IF(CUT_TREATMENT_AT(IJK)) THEN
367 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKM))) THEN
368 D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*DX(I)*DY(J)
369 ENDIF
370 ENDIF
371
372
373
374 IF (V_F >= ZERO) THEN
375 A_M(IJK,B,M) = D_F + FLUX_T(IJKM)
376 ELSE
377 A_M(IJK,B,M) = D_F
378 ENDIF
379 ENDIF
380 ENDIF
381
382 ENDIF
383 END DO
384
385 RETURN
386 END SUBROUTINE CONV_DIF_PHI0
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412 SUBROUTINE CONV_DIF_PHI_DC(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M,IER)
413
414
415
416
417
418
419
420
421 USE compar
422 USE discretization, ONLY: fpfoi_of
423 USE fun_avg
424 USE function3
425 USE functions
426 USE geometry
427 USE indices
428 USE matrix
429 USE parallel
430 USE param
431 USE param1
432 USE run
433 USE sendrecv
434 USE sendrecv3
435 USE toleranc
436 USE xsi
437 Use tmp_array
438 Use xsi_array
439 IMPLICIT NONE
440
441
442
443
444
445
446
447
448 DOUBLE PRECISION Phi(DIMENSION_3)
449
450
451 DOUBLE PRECISION Dif(DIMENSION_3)
452
453
454 INTEGER Disc
455
456
457 DOUBLE PRECISION Uf(DIMENSION_3), Vf(DIMENSION_3), Wf(DIMENSION_3)
458
459
460 DOUBLE PRECISION Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
461
462
463 INTEGER M
464
465
466 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
467
468
469 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
470
471
472 INTEGER IER
473
474
475 INTEGER I, J, K, IJK, IPJK, IJPK, IJKE, IJKN,&
476 IJKP, IJKT
477
478 INTEGER IMJK, IJKW
479 INTEGER IJMK, IJKS
480 INTEGER IJKM, IJKB
481 INTEGER IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
482 INTEGER IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
483
484
485 INTEGER incr
486
487
488
489 DOUBLE PRECISION V_F
490
491
492 DOUBLE PRECISION PHI_HO
493
494
495 DOUBLE PRECISION PHI_LO
496
497
498 DOUBLE PRECISION EAST_DC
499 DOUBLE PRECISION WEST_DC
500 DOUBLE PRECISION NORTH_DC
501 DOUBLE PRECISION SOUTH_DC
502 DOUBLE PRECISION TOP_DC
503 DOUBLE PRECISION BOTTOM_DC
504
505
506
507 call lock_xsi_array
508 call lock_tmp4_array
509
510
511
512
513
514 IF ( FPFOI ) THEN
515 Do IJK = ijkstart3, ijkend3
516 I = I_OF(IJK)
517 J = J_OF(IJK)
518 K = K_OF(IJK)
519 IJK4 = funijk3(I,J,K)
520 TMP4(IJK4) = PHI(IJK)
521 ENDDO
522 CALL send_recv3(tmp4)
523 ENDIF
524
525
526 =0
527
528
529 CALL CALC_XSI (DISC, PHI, UF, VF, WF, XSI_E, XSI_N, XSI_T,incr)
530
531
532
533
534
535
536
537
538
539
540
541
542
543 DO IJK = ijkstart3, ijkend3
544
545
546 = I_OF(IJK)
547 J = J_OF(IJK)
548 K = K_OF(IJK)
549
550 IF (FLUID_AT(IJK)) THEN
551
552
553 = IP_OF(IJK)
554 IMJK = IM_OF(IJK)
555 IJPK = JP_OF(IJK)
556 IJMK = JM_OF(IJK)
557 IJKP = KP_OF(IJK)
558 IJKM = KM_OF(IJK)
559 IJKE = EAST_OF(IJK)
560 IJKN = NORTH_OF(IJK)
561
562
563 = IP_OF(IP_OF(IPJK))
564 IPPP4 = funijk3(I_OF(IPPP), J_OF(IPPP), K_OF(IPPP))
565 IMMM = IM_OF(IM_OF(IMJK))
566 IMMM4 = funijk3(I_OF(IMMM), J_OF(IMMM), K_OF(IMMM))
567 JPPP = JP_OF(JP_OF(IJPK))
568 JPPP4 = funijk3(I_OF(JPPP), J_OF(JPPP), K_OF(JPPP))
569 JMMM = JM_OF(JM_OF(IJMK))
570 JMMM4 = funijk3(I_OF(JMMM), J_OF(JMMM), K_OF(JMMM))
571 KPPP = KP_OF(KP_OF(IJKP))
572 KPPP4 = funijk3(K_OF(IPPP), J_OF(KPPP), K_OF(KPPP))
573 KMMM = KM_OF(KM_OF(IJKM))
574 KMMM4 = funijk3(I_OF(KMMM), J_OF(KMMM), K_OF(KMMM))
575
576
577
578
579 = UF(IJK)
580 IF(V_F >= ZERO)THEN
581 PHI_LO = PHI(IJK)
582 IF ( FPFOI ) &
583 PHI_HO = FPFOI_OF(PHI(IPJK), PHI(IJK), &
584 PHI(IMJK), PHI(IM_OF(IMJK)))
585 ELSE
586 PHI_LO = PHI(IPJK)
587 IF ( FPFOI ) &
588 PHI_HO = FPFOI_OF(PHI(IJK), PHI(IPJK), &
589 PHI(IP_OF(IPJK)), TMP4(IPPP4))
590 ENDIF
591 IF (.NOT. FPFOI ) &
592 PHI_HO = XSI_E(IJK)*PHI(IPJK)+(1.0-XSI_E(IJK))*PHI(IJK)
593 EAST_DC = FLUX_E(IJK)*(PHI_LO - PHI_HO)
594
595
596
597
598 = VF(IJK)
599 IF(V_F >= ZERO)THEN
600 PHI_LO = PHI(IJK)
601 IF ( FPFOI ) &
602 PHI_HO = FPFOI_OF(PHI(IJPK), PHI(IJK), &
603 PHI(IJMK), PHI(JM_OF(IJMK)))
604 ELSE
605 PHI_LO = PHI(IJPK)
606 IF ( FPFOI ) &
607 PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJPK), &
608 PHI(JP_OF(IJPK)), TMP4(JPPP4))
609 ENDIF
610 IF (.NOT. FPFOI ) &
611 PHI_HO = XSI_N(IJK)*PHI(IJPK)+(1.0-XSI_N(IJK))*PHI(IJK)
612 NORTH_DC = FLUX_N(IJK)*(PHI_LO - PHI_HO)
613
614
615
616
617 IF (DO_K) THEN
618 IJKP = KP_OF(IJK)
619 IJKT = TOP_OF(IJK)
620 V_F = WF(IJK)
621 IF(V_F >= ZERO)THEN
622 PHI_LO = PHI(IJK)
623 IF ( FPFOI ) &
624 PHI_HO = FPFOI_OF(PHI(IJKP), PHI(IJK), &
625 PHI(IJKM), PHI(KM_OF(IJKM)))
626 ELSE
627 PHI_LO = PHI(IJKP)
628 IF ( FPFOI ) &
629 PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJKP), &
630 PHI(KP_OF(IJKP)), TMP4(KPPP4))
631 ENDIF
632 IF (.NOT. FPFOI ) &
633 PHI_HO = XSI_T(IJK)*PHI(IJKP)+(1.0-XSI_T(IJK))*PHI(IJK)
634 TOP_DC = FLUX_T(IJK)*(PHI_LO - PHI_HO)
635 ELSE
636 TOP_DC = ZERO
637 ENDIF
638
639
640
641
642 = IM_OF(IJK)
643 IJKW = WEST_OF(IJK)
644 V_F = UF(IMJK)
645 IF(V_F >= ZERO)THEN
646 PHI_LO = PHI(IMJK)
647 IF ( FPFOI ) &
648 PHI_HO = FPFOI_OF(PHI(IJK), PHI(IMJK), &
649 PHI(IM_OF(IMJK)), TMP4(IMMM4))
650 ELSE
651 PHI_LO = PHI(IJK)
652 IF ( FPFOI ) &
653 PHI_HO = FPFOI_OF(PHI(IMJK), PHI(IJK), &
654 PHI(IPJK), PHI(IP_OF(IPJK)))
655 ENDIF
656 IF (.NOT. FPFOI ) &
657 PHI_HO = XSI_E(IMJK)*PHI(IJK)+(ONE-XSI_E(IMJK))*PHI(IMJK)
658 WEST_DC = FLUX_E(IMJK)*(PHI_LO - PHI_HO)
659
660
661
662
663 = JM_OF(IJK)
664 IJKS = SOUTH_OF(IJK)
665 V_F = VF(IJMK)
666 IF(V_F >= ZERO)THEN
667 PHI_LO = PHI(IJMK)
668 IF ( FPFOI ) &
669 PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJMK), &
670 PHI(JM_OF(IJMK)), TMP4(JMMM4))
671 ELSE
672 PHI_LO = PHI(IJK)
673 IF ( FPFOI ) &
674 PHI_HO = FPFOI_OF(PHI(IJMK), PHI(IJK), &
675 PHI(IJPK), PHI(JP_OF(IJPK)))
676 ENDIF
677 IF (.NOT. FPFOI ) &
678 PHI_HO = XSI_N(IJMK)*PHI(IJK)+(ONE-XSI_N(IJMK))*PHI(IJMK)
679 SOUTH_DC = FLUX_N(IJMK)*(PHI_LO - PHI_HO)
680
681
682
683 IF (DO_K) THEN
684 IJKM = KM_OF(IJK)
685 IJKB = BOTTOM_OF(IJK)
686 V_F = WF(IJKM)
687 IF(V_F >= ZERO)THEN
688 PHI_LO = PHI(IJKM)
689 IF ( FPFOI ) &
690 PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJKM), &
691 PHI(KM_OF(IJKM)), TMP4(KMMM4))
692 ELSE
693 PHI_LO = PHI(IJK)
694 IF ( FPFOI ) &
695 PHI_HO = FPFOI_OF(PHI(IJKM), PHI(IJK), &
696 PHI(IJKP), PHI(KP_OF(IJKP)))
697 ENDIF
698 IF (.NOT. FPFOI ) &
699 PHI_HO = XSI_T(IJKM)*PHI(IJK)+(1.0-XSI_T(IJKM))*PHI(IJKM)
700 BOTTOM_DC = FLUX_T(IJKM)*(PHI_LO - PHI_HO)
701 ELSE
702 BOTTOM_DC = ZERO
703 ENDIF
704
705
706
707
708 (IJK,M) = B_M(IJK,M)+WEST_DC-EAST_DC+SOUTH_DC-NORTH_DC&
709 +BOTTOM_DC-TOP_DC
710
711 ENDIF
712 END DO
713 call unlock_tmp4_array
714 call unlock_xsi_array
715
716
717 RETURN
718 END SUBROUTINE CONV_DIF_PHI_DC
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746 SUBROUTINE CONV_DIF_PHI1(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M,IER)
747
748
749
750
751
752
753
754
755 USE param
756 USE param1
757 USE parallel
758 USE matrix
759 USE toleranc
760 USE run
761 USE geometry
762 USE compar
763 USE sendrecv
764 USE indices
765 USE vshear
766 USE xsi
767 USE xsi_array
768 USE fun_avg
769 USE functions
770
771
772
773 USE cutcell
774
775
776
777 IMPLICIT NONE
778
779
780
781
782
783
784
785
786 DOUBLE PRECISION Phi(DIMENSION_3)
787
788
789 DOUBLE PRECISION Dif(DIMENSION_3)
790
791
792 INTEGER Disc
793
794
795 DOUBLE PRECISION Uf(DIMENSION_3), Vf(DIMENSION_3), Wf(DIMENSION_3)
796
797
798 DOUBLE PRECISION Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
799
800
801 INTEGER M
802
803
804 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
805
806
807 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
808
809
810 INTEGER IER
811
812
813 INTEGER I, J, K, IJK, IPJK, IJPK, IJKE, IJKN,&
814 IJKP, IJKT
815
816 INTEGER IMJK, IM, IJKW
817 INTEGER IJMK, JM, IJKS
818 INTEGER IJKM, KM, IJKB
819
820 INTEGER incr
821
822
823
824
825 DOUBLE PRECISION D_f
826
827
828
829
830
831 call lock_xsi_array
832
833
834
835
836
837
838 =0
839
840
841 CALL CALC_XSI (DISC, PHI, UF, VF, WF, XSI_E, XSI_N, XSI_T,incr)
842
843
844
845
846 IF (SHEAR) THEN
847 DO IJK = ijkstart3, ijkend3
848 IF (FLUID_AT(IJK)) THEN
849 VF(IJK)=VF(IJK)+VSH(IJK)
850 END IF
851 END DO
852 END IF
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868 DO IJK = ijkstart3, ijkend3
869
870 = I_OF(IJK)
871 J = J_OF(IJK)
872 K = K_OF(IJK)
873
874 IF (FLUID_AT(IJK)) THEN
875
876 = IP_OF(IJK)
877 IJPK = JP_OF(IJK)
878 IJKE = EAST_OF(IJK)
879 IJKN = NORTH_OF(IJK)
880
881
882
883 = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*AYZ(IJK)
884
885
886
887 IF(CUT_TREATMENT_AT(IJK)) THEN
888 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IPJK))) THEN
889 D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*DY(J)*DZ(K)
890 ENDIF
891 ENDIF
892
893
894
895
896 (IJK,E,M) = D_F - XSI_E(IJK)*FLUX_E(IJK)
897
898 (IPJK,W,M) = D_F + (ONE - XSI_E(IJK))*FLUX_E(IJK)
899
900
901
902 = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*AXZ(IJK)
903
904
905
906 IF(CUT_TREATMENT_AT(IJK)) THEN
907 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJPK))) THEN
908 D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*DX(I)*DZ(K)
909 ENDIF
910 ENDIF
911
912
913
914
915 (IJK,N,M) = D_F - XSI_N(IJK)*FLUX_N(IJK)
916
917 (IJPK,S,M) = D_F + (ONE - XSI_N(IJK))*FLUX_N(IJK)
918
919
920
921 IF (DO_K) THEN
922 IJKP = KP_OF(IJK)
923 IJKT = TOP_OF(IJK)
924
925 = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*AXY(IJK)
926
927
928
929 IF(CUT_TREATMENT_AT(IJK)) THEN
930 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKP))) THEN
931 D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*DX(I)*DY(J)
932 ENDIF
933 ENDIF
934
935
936
937
938 (IJK,T,M) = D_F - XSI_T(IJK)*FLUX_T(IJK)
939
940 (IJKP,B,M)=D_F+(ONE-XSI_T(IJK))*FLUX_T(IJK)
941 ENDIF
942
943
944 = IM_OF(IJK)
945 IF (.NOT.FLUID_AT(IMJK)) THEN
946 IM = IM1(I)
947 IJKW = WEST_OF(IJK)
948
949 = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*AYZ(IMJK)
950
951
952
953 IF(CUT_TREATMENT_AT(IJK)) THEN
954 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IMJK))) THEN
955 D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*DY(J)*DZ(K)
956 ENDIF
957 ENDIF
958
959
960
961
962
963 (IJK,W,M) = D_F + (ONE - XSI_E(IMJK))*FLUX_E(IMJK)
964 ENDIF
965
966
967 = JM_OF(IJK)
968 IF (.NOT.FLUID_AT(IJMK)) THEN
969 JM = JM1(J)
970 IJKS = SOUTH_OF(IJK)
971
972 = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*AXZ(IJMK)
973
974
975
976 IF(CUT_TREATMENT_AT(IJK)) THEN
977 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJMK))) THEN
978 D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*DX(I)*DZ(K)
979 ENDIF
980 ENDIF
981
982
983
984
985
986 (IJK,S,M) = D_F + (ONE - XSI_N(IJMK))*FLUX_N(IJMK)
987 ENDIF
988
989
990 IF (DO_K) THEN
991 IJKM = KM_OF(IJK)
992 IF (.NOT.FLUID_AT(IJKM)) THEN
993 KM = KM1(K)
994 IJKB = BOTTOM_OF(IJK)
995
996 = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*AXY(&
997 IJKM)
998
999
1000
1001 IF(CUT_TREATMENT_AT(IJK)) THEN
1002 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKM))) THEN
1003 D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*DX(I)*DY(J)
1004 ENDIF
1005 ENDIF
1006
1007
1008
1009
1010 (IJK,B,M) = D_F + (ONE - XSI_T(IJKM))*FLUX_T(IJKM)
1011 ENDIF
1012 ENDIF
1013
1014 ENDIF
1015 END DO
1016
1017
1018 IF (SHEAR) THEN
1019 DO IJK = ijkstart3, ijkend3
1020 IF (FLUID_AT(IJK)) THEN
1021 VF(IJK)=VF(IJK)-VSH(IJK)
1022 END IF
1023 END DO
1024 END IF
1025
1026 call unlock_xsi_array
1027
1028
1029 RETURN
1030 END SUBROUTINE CONV_DIF_PHI1
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052 SUBROUTINE DIF_PHI_IS(DIF, A_M, B_M, M)
1053
1054
1055
1056
1057
1058
1059
1060
1061 USE param
1062 USE param1
1063 USE parallel
1064 USE matrix
1065 USE toleranc
1066 USE run
1067 USE geometry
1068 USE compar
1069 USE sendrecv
1070 USE indices
1071 USE scales
1072 USE constant
1073 USE physprop
1074 USE fldvar
1075 USE visc_s
1076 USE output
1077 USE is
1078 USE fun_avg
1079 USE function3
1080 USE functions
1081 IMPLICIT NONE
1082
1083
1084
1085
1086
1087
1088
1089
1090 INTEGER L
1091
1092
1093 INTEGER I, J, K, I1, I2, J1, J2, K1, K2, IJK,&
1094 IJKE, IJKN, IJKT, IPJK, IJPK, IJKP
1095
1096
1097 INTEGER M
1098
1099
1100 DOUBLE PRECISION Dif(DIMENSION_3)
1101
1102
1103 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
1104
1105
1106 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
1107
1108
1109 DOUBLE PRECISION D_f
1110
1111
1112
1113
1114 DO L = 1, DIMENSION_IS
1115 IF (IS_DEFINED(L)) THEN
1116 I1 = IS_I_W(L)
1117 I2 = IS_I_E(L)
1118 J1 = IS_J_S(L)
1119 J2 = IS_J_N(L)
1120 K1 = IS_K_B(L)
1121 K2 = IS_K_T(L)
1122
1123
1124
1125 IF(I1.LE.IEND2) I1 = MAX(I1, ISTART2)
1126
1127 IF(J1.LE.JEND2) J1 = MAX(J1, JSTART2)
1128
1129 IF(K1.LE.KEND2) K1 = MAX(K1, KSTART2)
1130
1131 IF(I2.GE.ISTART2) I2 = MIN(I2, IEND2)
1132
1133 IF(J2.GE.JSTART2) J2 = MIN(J2, JEND2)
1134
1135 IF(K2.GE.KSTART2) K2 = MIN(K2, KEND2)
1136
1137
1138 DO K = K1, K2
1139 DO J = J1, J2
1140 DO I = I1, I2
1141 IF (DEAD_CELL_AT(I,J,K)) CYCLE
1142 = FUNIJK(I,J,K)
1143
1144 SELECT CASE (TRIM(IS_PLANE(L)))
1145 CASE ('E')
1146 IJKE = EAST_OF(IJK)
1147 IPJK = IP_OF(IJK)
1148
1149 = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*AYZ(IJK)
1150
1151 (IJK,E,M) = A_M(IJK,E,M) - D_F
1152 A_M(IPJK,W,M) = A_M(IPJK,W,M) - D_F
1153
1154 CASE ('N')
1155 IJKN = NORTH_OF(IJK)
1156 IJPK = JP_OF(IJK)
1157
1158 = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*AXZ(IJK)
1159
1160 (IJK,N,M) = A_M(IJK,N,M) - D_F
1161 A_M(IJPK,S,M) = A_M(IJPK,S,M) - D_F
1162
1163 CASE ('T')
1164 IF (DO_K) THEN
1165 IJKT = TOP_OF(IJK)
1166 IJKP = KP_OF(IJK)
1167
1168 = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*&
1169 AXY(IJK)
1170
1171 (IJK,T,M) = A_M(IJK,T,M) - D_F
1172 A_M(IJKP,B,M) = A_M(IJKP,B,M) - D_F
1173
1174 ENDIF
1175 CASE DEFAULT
1176
1177 END SELECT
1178 END DO
1179 END DO
1180 END DO
1181 ENDIF
1182 END DO
1183
1184 RETURN
1185 END SUBROUTINE DIF_PHI_IS
1186