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