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