File: RELATIVE:/../../../mfix.git/model/des/dif_phi_bc_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_BC_DES(PHI, M, A_M, B_M)
10     
11           USE param
12           USE param1
13           USE matrix
14           USE geometry
15           USE indices
16           USE bc
17           USE compar
18           USE cutcell, only : CARTESIAN_GRID, CG_SAFE_MODE
19           USE fun_avg
20           USE functions
21     
22           IMPLICIT NONE
23     
24     !-----------------------------------------------
25     ! Dummy arguments
26     !-----------------------------------------------
27     
28     
29           DOUBLE PRECISION, INTENT(IN) :: PHI(DIMENSION_3)
30     ! Phase index
31           INTEGER, INTENT(IN) :: M
32     ! Septadiagonal matrix A_m
33           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
34     ! Vector b_m
35           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
36     !-----------------------------------------------
37     ! Local variables
38     !-----------------------------------------------
39     ! Boundary condition index
40           INTEGER :: L
41     ! Indices
42           INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK, &
43                      IM, JM, KM
44     !-----------------------------------------------
45     
46     ! Set up the default walls (i.e., bc_type='dummy' or undefined/default
47     ! boundaries) as non-conducting...
48     ! ---------------------------------------------------------------->>>
49     ! when setting up default walls do not use cutcells to avoid conflict
50     
51           IF(.NOT.CARTESIAN_GRID) THEN
52     
53     ! west yz plane
54              I1 = imin2
55              DO K1 = kmin3, kmax3
56              DO J1 = jmin3, jmax3
57                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
58                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
59                 IJK = FUNIJK(I1,J1,K1)
60                 IF (DEFAULT_WALL_AT(IJK)) THEN
61     ! Cutting the neighbor link between fluid cell and wall cell
62                    A_M(IP_OF(IJK),W,M) = ZERO
63     ! Setting the wall value equal to the adjacent fluid cell value
64                    A_M(IJK,:,M) = ZERO
65                    A_M(IJK,E,M) = ONE
66                    A_M(IJK,0,M) = -ONE
67                    B_M(IJK,M) = ZERO
68                 ENDIF
69              ENDDO
70              ENDDO
71     
72     ! east yz plane
73              I1 = IMAX2
74              DO K1 = kmin3, kmax3
75              DO J1 = jmin3, jmax3
76                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
77                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
78                 IJK = FUNIJK(I1,J1,K1)
79                 IF (DEFAULT_WALL_AT(IJK)) THEN
80     ! Cutting the neighbor link between fluid cell and wall cell
81                    A_M(IM_OF(IJK),E,M) = ZERO
82     ! Setting the wall value equal to the adjacent fluid cell value
83                    A_M(IJK,:,M) = ZERO
84                    A_M(IJK,W,M) = ONE
85                    A_M(IJK,0,M) = -ONE
86                    B_M(IJK,M) = ZERO
87                 ENDIF
88              ENDDO
89              ENDDO
90     
91     ! south xz plane
92              J1 = 1
93              DO K1 = kmin3, kmax3
94              DO I1 = imin3, imax3
95                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
96                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
97                 IJK = FUNIJK(I1,J1,K1)
98                 IF (DEFAULT_WALL_AT(IJK)) THEN
99     ! Cutting the neighbor link between fluid cell and wall cell
100                    A_M(JP_OF(IJK),S,M) = ZERO
101     ! Setting the wall value equal to the adjacent fluid cell value
102                    A_M(IJK,:,M) = ZERO
103                    A_M(IJK,N,M) = ONE
104                    A_M(IJK,0,M) = -ONE
105                    B_M(IJK,M) = ZERO
106                 ENDIF
107              ENDDO
108              ENDDO
109     
110     ! north xz plane
111              J1 = JMAX2
112              DO K1 = kmin3, kmax3
113              DO I1 = imin3, imax3
114                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
115                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
116                 IJK = FUNIJK(I1,J1,K1)
117                 IF (DEFAULT_WALL_AT(IJK)) THEN
118     ! Cutting the neighbor link between fluid cell and wall cell
119                    A_M(JM_OF(IJK),N,M) = ZERO
120     ! Setting the wall value equal to the adjacent fluid cell value
121                    A_M(IJK,:,M) = ZERO
122                    A_M(IJK,S,M) = ONE
123                    A_M(IJK,0,M) = -ONE
124                    B_M(IJK,M) = ZERO
125                 ENDIF
126              ENDDO
127              ENDDO
128     
129              IF (DO_K) THEN
130     ! bottom xy plane
131                 K1 = 1
132                 DO J1 = jmin3, jmax3
133                 DO I1 = imin3, imax3
134                    IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
135                    IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
136                    IJK = FUNIJK(I1,J1,K1)
137     
138                    IF (DEFAULT_WALL_AT(IJK)) THEN
139     ! Cutting the neighbor link between fluid cell and wall cell
140                       A_M(KP_OF(IJK),B,M) = ZERO
141     ! Setting the wall value equal to the adjacent fluid cell value (set
142     ! the boundary cell value equal to adjacent fluid cell value)
143                       A_M(IJK,:,M) = ZERO
144                       A_M(IJK,T,M) = ONE
145                       A_M(IJK,0,M) = -ONE
146                       B_M(IJK,M) = ZERO
147                    ENDIF
148                 ENDDO
149                 ENDDO
150     
151     ! top xy plane
152                 K1 = KMAX2
153                 DO J1 = jmin3, jmax3
154                 DO I1 = imin3, imax3
155                    IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
156                    IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
157                    IJK = FUNIJK(I1,J1,K1)
158                    IF (DEFAULT_WALL_AT(IJK)) THEN
159     ! Cutting the neighbor link between fluid cell and wall cell
160                       A_M(KM_OF(IJK),T,M) = ZERO
161     ! Setting the wall value equal to the adjacent fluid cell value
162                       A_M(IJK,:,M) = ZERO
163                       A_M(IJK,B,M) = ONE
164                       A_M(IJK,0,M) = -ONE
165                       B_M(IJK,M) = ZERO
166                    ENDIF
167                 ENDDO
168                 ENDDO
169              ENDIF
170     
171           ENDIF
172     
173     ! Set user defined wall boundary conditions
174     ! ---------------------------------------------------------------->>>
175           DO L = 1, DIMENSION_BC
176              IF (BC_DEFINED(L)) THEN
177     
178                 I1 = BC_I_W(L)
179                 I2 = BC_I_E(L)
180                 J1 = BC_J_S(L)
181                 J2 = BC_J_N(L)
182                 K1 = BC_K_B(L)
183                 K2 = BC_K_T(L)
184     
185                 DO K = K1, K2
186                 DO J = J1, J2
187                 DO I = I1, I2
188     
189                    IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
190                    IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
191     
192                    IJK = FUNIJK(I,J,K)
193     
194     
195     ! Set the boundary cell equal to the known value in that
196     ! cell
197                    A_M(IJK,:,M) = ZERO
198                    B_M(IJK,M) = ZERO
199     ! second modify the matrix equation according to the user specified
200     ! boundary condition
201                    IF (FLUID_AT(EAST_OF(IJK))) THEN
202                       A_M(IJK,0,M) = -ODX_E(I)
203                       A_M(IJK,E,M) =  ODX_E(I)
204     
205                    ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
206                       IM = IM1(I)
207                       A_M(IJK,W,M) =  ODX_E(IM)
208                       A_M(IJK,0,M) = -ODX_E(IM)
209     
210                    ELSEIF (FLUID_AT(NORTH_OF(IJK))) THEN
211                       A_M(IJK,0,M) = -ODY_N(J)
212                       A_M(IJK,N,M) =  ODY_N(J)
213     
214                    ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
215                       JM = JM1(J)
216                       A_M(IJK,S,M) =  ODY_N(JM)
217                       A_M(IJK,0,M) = -ODY_N(JM)
218     
219                    ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
220                       A_M(IJK,0,M) = -OX(I)*ODZ_T(K)
221                       A_M(IJK,T,M) =  OX(I)*ODZ_T(K)
222     
223                    ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
224                       KM = KM1(K)
225                       A_M(IJK,B,M) =  OX(I)*ODZ_T(KM)
226                       A_M(IJK,0,M) = -OX(I)*ODZ_T(KM)
227     
228                    ENDIF
229                 ENDDO
230                 ENDDO
231                 ENDDO
232              ENDIF   ! end if (bc_defined)
233           ENDDO   ! end L do loop over dimension_bc
234     
235     ! modifications for cartesian grid implementation
236           IF(CARTESIAN_GRID .AND. .NOT.(CG_SAFE_MODE(1)==1)) &
237              CALL DIF_PHI_BC_DES_CG(PHI, 0, A_M, B_M)
238     
239           RETURN
240           END SUBROUTINE DIF_PHI_BC_DES
241     
242     
243     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
244     !                                                                      C
245     !  Subroutine: BC_PHI_CG                                               C
246     !  Purpose: Modify boundary conditions for cartesian grid cut-cell     C
247     !           implementation                                             C
248     !                                                                      C
249     !  Author: Jeff Dietiker                                               C
250     !                                                                      C
251     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
252           SUBROUTINE DIF_PHI_BC_DES_CG(PHI, M, A_M, B_M)
253     
254     !-----------------------------------------------
255     ! Modules
256     !-----------------------------------------------
257           USE param
258           USE param1
259           USE matrix
260           USE geometry
261           USE indices
262           USE bc
263           USE compar
264           USE cutcell
265           USE fun_avg
266           USE functions
267     
268           IMPLICIT NONE
269     
270     !-----------------------------------------------
271     ! Dummy arguments
272     !-----------------------------------------------
273     ! The field variable being solved for:
274     !     e.g., T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
275     !     e_Turb_G
276           DOUBLE PRECISION, INTENT(IN) :: PHI(DIMENSION_3)
277     
278     ! Phase index
279           INTEGER, INTENT(IN) :: M
280     ! Septadiagonal matrix A_m
281           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
282     ! Vector b_m
283           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
284     !-----------------------------------------------
285     ! Local variables
286     !-----------------------------------------------
287     ! Indices
288           INTEGER :: I, J, K, IJK, IM, JM, KM
289     
290           LOGICAL :: ALONG_GLOBAL_GHOST_LAYER
291     
292     !-----------------------------------------------
293     
294           DO IJK = ijkstart3, ijkend3
295     
296              I = I_OF(IJK)
297              J = J_OF(IJK)
298              K = K_OF(IJK)
299     
300              ALONG_GLOBAL_GHOST_LAYER = (I < IMIN1) .OR. ( I>IMAX1 )       &
301                 .OR. (J < JMIN1) .OR. (J > JMAX1)
302     
303              IF(DO_K) ALONG_GLOBAL_GHOST_LAYER = ALONG_GLOBAL_GHOST_LAYER  &
304                 .OR. (K < KMIN1) .OR. (K > KMAX1)
305     
306     ! Setting the value in the boundary cell equal to what is known
307              IF(BLOCKED_CELL_AT(IJK)) THEN
308                 A_M(IJK,:,M) = ZERO
309                 A_M(IJK,0,M) = -ONE
310                 B_M(IJK,M) = -PHI(IJK)
311              ENDIF
312     
313              IF(BLOCKED_CELL_AT(IJK).OR.ALONG_GLOBAL_GHOST_LAYER) THEN
314     
315     
316                 IF (CUT_CELL_AT(IP_OF(IJK))) THEN
317                    IF( BC_ID(IP_OF(IJK)) > 0 ) THEN
318                       A_M(IJK,0,M) = -ODX_E(I)
319                       A_M(IJK,E,M) =  ODX_E(I)
320                       B_M(IJK,M) =  ZERO
321                    ENDIF
322     
323                 ELSEIF (CUT_CELL_AT(IM_OF(IJK))) THEN
324                    IF(BC_ID(IM_OF(IJK)) > 0 ) THEN
325                       IM = IM1(I)
326                       A_M(IJK,W,M) =  ODX_E(IM)
327                       A_M(IJK,0,M) = -ODX_E(IM)
328                       B_M(IJK,M) = ZERO
329                    ENDIF
330     
331                 ELSEIF (CUT_CELL_AT(JP_OF(IJK))) THEN
332                    IF(BC_ID(JP_OF(IJK)) > 0 ) THEN
333                       A_M(IJK,0,M) = -ODY_N(J)
334                       A_M(IJK,N,M) =  ODY_N(J)
335                       B_M(IJK,M) = ZERO
336                    ENDIF
337     
338                 ELSEIF (CUT_CELL_AT(JM_OF(IJK))) THEN
339                    IF(BC_ID(JM_OF(IJK)) > 0 ) THEN
340                       JM = JM1(J)
341                       A_M(IJK,S,M) =  ODY_N(JM)
342                       A_M(IJK,0,M) = -ODY_N(JM)
343                       B_M(IJK,M) = ZERO
344                    ENDIF
345     
346                 ELSEIF (CUT_CELL_AT(KP_OF(IJK))) THEN
347                    IF(BC_ID(KP_OF(IJK)) > 0 ) THEN
348                       A_M(IJK,0,M) = -OX(I)*ODZ_T(K)
349                       A_M(IJK,T,M) =  OX(I)*ODZ_T(K)
350                       B_M(IJK,M)  = ZERO
351                    ENDIF
352     
353                 ELSEIF (CUT_CELL_AT(KM_OF(IJK))) THEN
354                    IF(BC_ID(KM_OF(IJK)) > 0 ) THEN
355                       KM = KM1(K)
356                       A_M(IJK,B,M) =  OX(I)*ODZ_T(KM)
357                       A_M(IJK,0,M) = -OX(I)*ODZ_T(KM)
358                       B_M(IJK,M) = ZERO
359                    ENDIF
360                 ENDIF
361              ENDIF
362           ENDDO
363     
364           RETURN
365     
366           END SUBROUTINE DIF_PHI_BC_DES_CG
367     
368     
369