File: /nfs/home/0/users/jenkins/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, IER)
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           INTEGER :: IER
42     
43     ! Fluid Cell indices
44           INTEGER :: I, J, K, IJK
45           INTEGER :: IMJK, IM, IJKW, IPJK, IJKE
46           INTEGER :: IJMK, JM, IJKS, IJPK, IJKN
47           INTEGER :: IJKM, KM, IJKB, IJKP, IJKT
48     !
49     ! Difusion parameter
50           DOUBLE PRECISION :: D_f
51     
52     
53     !  Calculate convection-diffusion fluxes through each of the faces
54           DO IJK = ijkstart3, ijkend3
55     
56              I = I_OF(IJK)
57              J = J_OF(IJK)
58              K = K_OF(IJK)
59     
60              IF(.NOT.FLUID_AT(IJK)) CYCLE
61     
62              IPJK = IP_OF(IJK)
63              IJPK = JP_OF(IJK)
64     
65              IJKE = EAST_OF(IJK)
66              IJKN = NORTH_OF(IJK)
67     
68     ! East face (i+1/2, j, k)
69              D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*AYZ(IJK)
70              IF(CUT_TREATMENT_AT(IJK)) THEN
71                 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IPJK))) THEN
72                    D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*DY(J)*DZ(K)
73                 ENDIF
74              ENDIF
75     
76              A_M(IJK,E,M) = D_F
77              A_M(IPJK,W,M) = D_F
78     
79     ! West face (i-1/2, j, k)
80              IMJK = IM_OF(IJK)
81              IF (.NOT.FLUID_AT(IMJK)) THEN
82                 IM = IM1(I)
83                 IJKW = WEST_OF(IJK)
84                 D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*AYZ(IMJK)
85                 IF(CUT_TREATMENT_AT(IJK)) THEN
86                    IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IMJK))) THEN
87                       D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*DY(J)*DZ(K)
88                    ENDIF
89                 ENDIF
90                 A_M(IJK,W,M) = D_F
91              ENDIF
92     
93     
94     ! North face (i, j+1/2, k)
95              D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*AXZ(IJK)
96              IF(CUT_TREATMENT_AT(IJK)) THEN
97                 IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJPK))) THEN
98                    D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*DX(I)*DZ(K)
99                 ENDIF
100              ENDIF
101              A_M(IJK,N,M) = D_F
102              A_M(IJPK,S,M) = D_F
103     
104     ! South face (i, j-1/2, k)
105              IJMK = JM_OF(IJK)
106              IF (.NOT.FLUID_AT(IJMK)) THEN
107                 JM = JM1(J)
108                 IJKS = SOUTH_OF(IJK)
109                 D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*AXZ(IJMK)
110                 IF(CUT_TREATMENT_AT(IJK)) THEN
111                    IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJMK))) THEN
112                       D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*DX(I)*DZ(K)
113                    ENDIF
114                 ENDIF
115                 A_M(IJK,S,M) = D_F
116              ENDIF
117     
118     
119     ! Top face (i, j, k+1/2)
120              IF (DO_K) THEN
121                 IJKP = KP_OF(IJK)
122                 IJKT = TOP_OF(IJK)
123                 D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*AXY(IJK)
124                 IF(CUT_TREATMENT_AT(IJK)) THEN
125                    IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKP))) THEN
126                       D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*DX(I)*DY(J)
127                    ENDIF
128                 ENDIF
129                 A_M(IJK,T,M) = D_F
130                 A_M(IJKP,B,M) = D_F
131     
132     
133     ! Bottom face (i, j, k-1/2)
134                 IJKM = KM_OF(IJK)
135                 IF (.NOT.FLUID_AT(IJKM)) THEN
136                    KM = KM1(K)
137                    IJKB = BOTTOM_OF(IJK)
138                    D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*AXY(IJKM)
139                    IF(CUT_TREATMENT_AT(IJK)) THEN
140                       IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKM))) THEN
141                          D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*DX(I)*DY(J)
142                       ENDIF
143                    ENDIF
144                    A_M(IJK,B,M) = D_F
145                 ENDIF
146              ENDIF
147           END DO
148     
149           CALL DIF_PHI_IS(DIF, A_M, B_M, M, IER)
150     
151     
152           RETURN
153           END SUBROUTINE DIF_PHI_DES
154