File: RELATIVE:/../../../mfix.git/model/des/dif_phi_des.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: DIFFUSE_MEAN_FIELDS                                     !
4     !  Author: J.Musser                                   Date: 11-NOV-14  !
5     !                                                                      !
6     !  Purpose:                                                            !
7     !                                                                      !
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9           SUBROUTINE DIF_PHI_DES(M, DIF, A_M, B_M)
10     
11     !-----------------------------------------------
12     !   M o d u l e s
13     !-----------------------------------------------
14           USE param
15           USE param1
16           USE parallel
17           USE matrix
18           USE toleranc
19           USE run
20           USE geometry
21           USE compar
22           USE sendrecv
23           USE indices
24           USE fun_avg
25           USE functions
26           USE cutcell
27     
28           IMPLICIT NONE
29     
30     ! Phase index
31           INTEGER, INTENT(IN) :: M
32     
33     !  Septadiagonal matrix A_m
34           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
35     
36     ! Vector b_m
37           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
38     
39           DOUBLE PRECISION, INTENT(IN) :: DIF(DIMENSION_3)
40     
41     ! Fluid Cell indices
42           INTEGER :: I, J, K, IJK
43           INTEGER :: IMJK, IM, IJKW, IPJK, IJKE
44           INTEGER :: IJMK, JM, IJKS, IJPK, IJKN
45           INTEGER :: IJKM, KM, IJKB, IJKP, IJKT
46     !
47     ! Diffusion parameter
48           DOUBLE PRECISION :: D_f
49     
50     
51     !  Calculate convection-diffusion fluxes through each of the faces
52           DO IJK = ijkstart3, ijkend3
53     
54              I = I_OF(IJK)
55              J = J_OF(IJK)
56              K = K_OF(IJK)
57     
58              IF(.NOT.FLUID_AT(IJK)) CYCLE
59     
60              IPJK = IP_OF(IJK)
61              IJPK = JP_OF(IJK)
62     
63              IJKE = EAST_OF(IJK)
64              IJKN = NORTH_OF(IJK)
65     
66     ! East face (i+1/2, j, k)
67              D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*AYZ(IJK)
68              IF(CUT_TREATMENT_AT(IJK)) THEN
69                 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IPJK))) THEN
70                    D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*DY(J)*DZ(K)
71                 ENDIF
72              ENDIF
73     
74              A_M(IJK,E,M) = D_F
75              A_M(IPJK,W,M) = D_F
76     
77     ! West face (i-1/2, j, k)
78              IMJK = IM_OF(IJK)
79              IF (.NOT.FLUID_AT(IMJK)) THEN
80                 IM = IM1(I)
81                 IJKW = WEST_OF(IJK)
82                 D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*AYZ(IMJK)
83                 IF(CUT_TREATMENT_AT(IJK)) THEN
84                    IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IMJK))) THEN
85                       D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*DY(J)*DZ(K)
86                    ENDIF
87                 ENDIF
88                 A_M(IJK,W,M) = D_F
89              ENDIF
90     
91     
92     ! North face (i, j+1/2, k)
93              D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*AXZ(IJK)
94              IF(CUT_TREATMENT_AT(IJK)) THEN
95                 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJPK))) THEN
96                    D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*DX(I)*DZ(K)
97                 ENDIF
98              ENDIF
99              A_M(IJK,N,M) = D_F
100              A_M(IJPK,S,M) = D_F
101     
102     ! South face (i, j-1/2, k)
103              IJMK = JM_OF(IJK)
104              IF (.NOT.FLUID_AT(IJMK)) THEN
105                 JM = JM1(J)
106                 IJKS = SOUTH_OF(IJK)
107                 D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*AXZ(IJMK)
108                 IF(CUT_TREATMENT_AT(IJK)) THEN
109                    IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJMK))) THEN
110                       D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*DX(I)*DZ(K)
111                    ENDIF
112                 ENDIF
113                 A_M(IJK,S,M) = D_F
114              ENDIF
115     
116     
117     ! Top face (i, j, k+1/2)
118              IF (DO_K) THEN
119                 IJKP = KP_OF(IJK)
120                 IJKT = TOP_OF(IJK)
121                 D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*AXY(IJK)
122                 IF(CUT_TREATMENT_AT(IJK)) THEN
123                    IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKP))) THEN
124                       D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*DX(I)*DY(J)
125                    ENDIF
126                 ENDIF
127                 A_M(IJK,T,M) = D_F
128                 A_M(IJKP,B,M) = D_F
129     
130     
131     ! Bottom face (i, j, k-1/2)
132                 IJKM = KM_OF(IJK)
133                 IF (.NOT.FLUID_AT(IJKM)) THEN
134                    KM = KM1(K)
135                    IJKB = BOTTOM_OF(IJK)
136                    D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*AXY(IJKM)
137                    IF(CUT_TREATMENT_AT(IJK)) THEN
138                       IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKM))) THEN
139                          D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*DX(I)*DY(J)
140                       ENDIF
141                    ENDIF
142                    A_M(IJK,B,M) = D_F
143                 ENDIF
144              ENDIF
145           END DO
146     
147           CALL DIF_PHI_IS(DIF, A_M, M)
148     
149     
150           RETURN
151           END SUBROUTINE DIF_PHI_DES
152