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