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