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