File: N:\mfix\model\des\gas_drag.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: GAS_DRAG_W                                              !
4     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
5     !                                                                      !
6     !  Purpose: Account for the equal and opposite drag force on the gas   !
7     !           phase due to particles by introducing the drag as a        !
8     !           source term.  Face centered.                               !
9     !                                                                      !
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
11           SUBROUTINE GAS_DRAG_U(A_M, B_M, IER)
12     
13     ! Global Variables:
14     !---------------------------------------------------------------------//
15     ! Runtime Flag: Interpolate DES values
16           use particle_filter, only: DES_INTERP_SCHEME_ENUM, DES_INTERP_GARG
17     ! Flag: Gas sees the effect of particles in gas/solids flows.
18           use discretelement, only: DES_ONEWAY_COUPLED
19     ! Coefficient at cell corners added to the gas momentum A matix.
20           use discretelement, only: DRAG_AM, F_GDS
21     ! Coefficient at cell corners added to gas momentum B vector.
22           use discretelement, only: DRAG_BM
23     ! Volume of X-momentum cell
24           use geometry, only: VOL_U
25     ! Flag to calculate Z direction
26           use geometry, only: DO_K
27     
28     ! Global Parameters:
29     !---------------------------------------------------------------------//
30     ! Size of IJK arrays and solids phase..
31           use param, only: DIMENSION_3, DIMENSION_M
32     ! Fluid grid loop bounds.
33           use compar, only: IJKStart3, IJKEnd3
34     ! The I, J, and K values that comprise an IJK
35           use indices, only: I_OF, J_OF, K_OF
36     ! Flag: Fluid exists at indexed cell
37           use functions, only: FLUID_AT
38     ! IJK of cell to east.
39           use functions, only: EAST_OF
40     ! IJK function for I,J,K that includes mapped indices.
41           use compar, only: FUNIJK_MAP_C
42     ! Function for averaging to a scalar cell's east face.
43           use fun_avg, only: AVG_X
44     ! Domain index bounds.
45           use compar, only: ISTART2, JSTART2, KSTART2
46           use compar, only: IEND2, JEND2, KEND2
47     
48           IMPLICIT NONE
49     
50     ! Dummy Arguments:
51     !---------------------------------------------------------------------//
52     ! Septadiagonal matrix A_m
53           DOUBLE PRECISION, INTENT(INOUT) :: A_M(DIMENSION_3,-3:3,0:DIMENSION_M)
54     ! Vector b_m
55           DOUBLE PRECISION, INTENT(INOUT) :: B_M(DIMENSION_3, 0:DIMENSION_M)
56     ! Error index
57           INTEGER, INTENT(INOUT) :: IER
58     
59     
60     ! Local Variables:
61     !---------------------------------------------------------------------//
62     ! Grid cell indices
63           INTEGER :: I, J, K, IJK, IJMK, IJKM, IJMKM, IJKE
64     ! temporary variables for matrix A_M and vector B_M
65           DOUBLE PRECISION :: tmp_A, tmp_B
66     ! Averaging factor
67           DOUBLE PRECISION :: AVG_FACTOR
68     !......................................................................!
69     
70     ! Initialize error flag.
71           IER = 0
72     
73     ! Skip this routine if the gas/solids are only one-way coupled.
74           IF(DES_ONEWAY_COUPLED) RETURN
75     
76     ! Average the interpoalted drag force from the cell corners to the cell face.
77           IF(DES_INTERP_SCHEME_ENUM == DES_INTERP_GARG)THEN
78     
79              AVG_FACTOR = merge(0.25d0, 0.5d0, DO_K)
80     
81     !!$omp parallel do schedule(guided,50) default(none)                    &
82     !!$omp shared(IJKSTART3, IJKEND3, FUNIJK_MAP_C, I_OF, J_OF,             &
83     !!$omp    K_OF, DO_K, AVG_FACTOR, DRAG_AM, DRAG_BM, A_M, B_M, VOL_U)    &
84     !!$omp private(IJK, I, J, K, IJMK, IJKM, IJMKM, tmp_A, tmp_B)
85              DO IJK = IJKSTART3, IJKEND3
86     
87                 IF(.NOT.FLUID_AT(IJK)) CYCLE
88     
89                 I = I_OF(IJK)
90                 J = J_OF(IJK)
91                 K = K_OF(IJK)
92     
93                 IF (I.LT.ISTART2 .OR. I.GT.IEND2) CYCLE
94                 IF (J.LT.JSTART2 .OR. J.GT.JEND2) CYCLE
95                 IF (K.LT.KSTART2 .OR. K.GT.KEND2) CYCLE
96     
97                 IJMK = FUNIJK_MAP_C(I, J-1, K)
98     
99                 tmp_A = -AVG_FACTOR*(DRAG_AM(IJK) + DRAG_AM(IJMK))
100                 tmp_B = -AVG_FACTOR*(DRAG_BM(IJK,1) + DRAG_BM(IJMK,1))
101     
102                 IF(DO_K) THEN
103                    IJKM = FUNIJK_MAP_C(I, J, K-1)
104                    IJMKM = FUNIJK_MAP_C(I, J-1, K-1)
105                    tmp_A = tmp_A - AVG_FACTOR*                             &
106                       (DRAG_AM(IJKM) + DRAG_AM(IJMKM))
107                    tmp_B = tmp_B - AVG_FACTOR*                             &
108                       (DRAG_BM(IJKM,1) + DRAG_BM(IJMKM,1))
109                 ENDIF
110     
111                 A_M(IJK,0,0) = A_M(IJK,0,0) + tmp_A*VOL_U(IJK)
112                 B_M(IJK,0) = B_M(IJK,0) + tmp_B*VOL_U(IJK)
113     
114              ENDDO
115     !!$omp end parallel do
116     
117           ELSE
118     
119     !$omp parallel do default(none) &
120     !$omp private(IJK, I, IJKE, tmp_A, tmp_B) &
121     !$omp shared(IJKSTART3,IJKEND3,I_OF, F_GDS, DRAG_BM, A_M, B_M, VOL_U)
122              DO IJK = IJKSTART3, IJKEND3
123                 IF(FLUID_AT(IJK)) THEN
124                    I = I_OF(IJK)
125                    IJKE = EAST_OF(IJK)
126     
127                    tmp_A = AVG_X(F_GDS(IJK), F_GDS(IJKE), I)
128                    tmp_B = AVG_X(DRAG_BM(IJK,1), DRAG_BM(IJKE,1), I)
129     
130                    A_M(IJK,0,0) = A_M(IJK,0,0) - VOL_U(IJK) * tmp_A
131                    B_M(IJK,0) = B_M(IJK,0) - VOL_U(IJK) * tmp_B
132                 ENDIF
133              ENDDO
134     !$omp end parallel do
135           ENDIF
136     
137           END SUBROUTINE GAS_DRAG_U
138     
139     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
140     !                                                                      !
141     !  Subroutine: GAS_DRAG_V                                              !
142     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
143     !                                                                      !
144     !  Purpose: Account for the equal and opposite drag force on the gas   !
145     !           phase due to particles by introducing the drag as a        !
146     !           source term.  Face centered.                               !
147     !                                                                      !
148     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
149           SUBROUTINE GAS_DRAG_V(A_M, B_M, IER)
150     
151     
152     ! Global Variables:
153     !---------------------------------------------------------------------//
154     ! Runtime Flag: Interpolate DES values
155           use particle_filter, only: DES_INTERP_SCHEME_ENUM, DES_INTERP_GARG
156     ! Flag: Gas sees the effect of particles in gas/solids flows.
157           use discretelement, only: DES_ONEWAY_COUPLED
158     ! Coefficient at cell corners added to the gas momentum A matix.
159           use discretelement, only: DRAG_AM, F_GDS
160     ! Coefficient at cell corners added to gas momentum B vector.
161           use discretelement, only: DRAG_BM
162     ! Volume of Y-momentum cell
163           use geometry, only: VOL_V
164     ! Flag to calculate Z direction
165           use geometry, only: DO_K
166     
167     ! Global Parameters:
168     !---------------------------------------------------------------------//
169     ! Size of IJK arrays and solids phase..
170           use param, only: DIMENSION_3, DIMENSION_M
171     ! Fluid grid loop bounds.
172           use compar, only: IJKStart3, IJKEnd3
173     ! The I, J, and K values that comprise an IJK
174           use indices, only: I_OF, J_OF, K_OF
175     ! Flag: Fluid exists at indexed cell
176           use functions, only: FLUID_AT
177     ! IJK of cell to north.
178           use functions, only: NORTH_OF
179     ! IJK function for I,J,K that includes mapped indices.
180           use compar, only: FUNIJK_MAP_C
181     ! Function for averaging to a scalar cell's north face.
182           use fun_avg, only: AVG_Y
183     ! Domain index bounds.
184           use compar, only: ISTART2, JSTART2, KSTART2
185           use compar, only: IEND2, JEND2, KEND2
186     
187     
188           IMPLICIT NONE
189     
190     ! Dummy Arguments:
191     !---------------------------------------------------------------------//
192     ! Septadiagonal matrix A_m
193           DOUBLE PRECISION, INTENT(INOUT) :: A_M(DIMENSION_3, -3:3, 0:DIMENSION_M)
194     ! Vector b_m
195           DOUBLE PRECISION, INTENT(INOUT) :: B_M(DIMENSION_3, 0:DIMENSION_M)
196     ! Error index
197           INTEGER, INTENT(INOUT) :: IER
198     
199     ! Local Variables:
200     !---------------------------------------------------------------------//
201     ! Grid cell indices
202           INTEGER :: I, J, K, IJK, IMJK, IJKM, IMJKM, IJKN
203     ! temporary variables for matrix A_M and vector B_M
204           DOUBLE PRECISION tmp_A, tmp_B
205     ! Averaging factor
206           DOUBLE PRECISION :: AVG_FACTOR
207     !......................................................................!
208     
209     ! Initialize error flag.
210           IER = 0
211     
212     ! Skip this routine if the gas/solids are only one-way coupled.
213           IF(DES_ONEWAY_COUPLED) RETURN
214     
215           IF(DES_INTERP_SCHEME_ENUM == DES_INTERP_GARG)THEN
216     
217              AVG_FACTOR = merge(0.25d0, 0.5d0, DO_K)
218     
219     !!$omp parallel do schedule (guided,50) default(none)                   &
220     !!$omp shared(IJKSTART3, IJKEND3, FUNIJK_MAP_C, I_OF, J_OF,             &
221     !!$omp    K_OF, DO_K, AVG_FACTOR, DRAG_AM, DRAG_BM, A_M, B_M, VOL_V)    &
222     !!$omp private(IJK, I, J, K, IMJK, IJKM, IMJKM, tmp_A, tmp_B)
223              DO IJK = IJKSTART3, IJKEND3
224                 IF(.NOT.FLUID_AT(IJK)) CYCLE
225     
226                 I = I_OF(IJK)
227                 J = J_OF(IJK)
228                 K = K_OF(IJK)
229     
230                 IF (I.LT.ISTART2 .OR. I.GT.IEND2) CYCLE
231                 IF (J.LT.JSTART2 .OR. J.GT.JEND2) CYCLE
232                 IF (K.LT.KSTART2 .OR. K.GT.KEND2) CYCLE
233     
234                 IMJK = FUNIJK_MAP_C(I-1,J,K)
235     
236                 tmp_A = -AVG_FACTOR*(DRAG_AM(IJK) + DRAG_AM(IMJK))
237                 tmp_B = -AVG_FACTOR*(DRAG_BM(IJK,2) + DRAG_BM(IMJK,2))
238     
239                 IF(DO_K) THEN
240     
241                    IJKM = FUNIJK_MAP_C(I,J,K-1)
242                    IMJKM = FUNIJK_MAP_C(I-1,J,K-1)
243     
244                    tmp_A = tmp_A - AVG_FACTOR*                             &
245                       (DRAG_AM(IJKM) + DRAG_AM(IMJKM))
246                    tmp_B = tmp_B - AVG_FACTOR*                             &
247                       (DRAG_BM(IJKM,2) + DRAG_BM(IMJKM,2))
248                 ENDIF
249     
250                 A_M(IJK,0,0) = A_M(IJK,0,0) + tmp_A*VOL_V(IJK)
251                 B_M(IJK,0) = B_M(IJK,0) + tmp_B*VOL_V(IJK)
252     
253              ENDDO
254     !!$omp end parallel do
255     
256           ELSE
257     
258     !$omp parallel do default(none) &
259     !$omp private(IJK, J, IJKN, tmp_a, tmp_B) &
260     !$omp shared(IJKSTART3, IJKEND3, J_OF, F_GDS, DRAG_AM, DRAG_BM, A_M, B_M, VOL_V)
261              DO IJK = IJKSTART3, IJKEND3
262                 IF(FLUID_AT(IJK)) THEN
263                    J = J_OF(IJK)
264                    IJKN = NORTH_OF(IJK)
265     
266                    tmp_A = AVG_Y(F_GDS(IJK), F_GDS(IJKN), J)
267                    tmp_B = AVG_Y(DRAG_BM(IJK,2), DRAG_BM(IJKN,2), J)
268     
269                    A_M(IJK,0,0) = A_M(IJK,0,0) - VOL_V(IJK) * tmp_A
270                    B_M(IJK,0) = B_M(IJK,0) - VOL_V(IJK) * tmp_B
271                 ENDIF
272              ENDDO
273     !$omp end parallel do
274           ENDIF
275     
276           END SUBROUTINE GAS_DRAG_V
277     
278     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
279     !                                                                      !
280     !  Subroutine: GAS_DRAG_W                                              !
281     !  Author: Jay Boyalakuntla                           Date: 12-Jun-04  !
282     !                                                                      !
283     !  Purpose: Account for the equal and opposite drag force on the gas   !
284     !           phase due to particles by introducing the drag as a        !
285     !           source term.  Face centered.                               !
286     !                                                                      !
287     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
288           SUBROUTINE GAS_DRAG_W(A_M, B_M, IER)
289     
290     ! Global Variables:
291     !---------------------------------------------------------------------//
292     ! Runtime Flag: Interpolate DES values
293           use particle_filter, only: DES_INTERP_SCHEME_ENUM, DES_INTERP_GARG
294     ! Flag: Gas sees the effect of particles in gas/solids flows.
295           use discretelement, only: DES_ONEWAY_COUPLED
296     ! Coefficient at cell corners added to the gas momentum A matix.
297           use discretelement, only: DRAG_AM, F_GDS
298     ! Coefficient at cell corners added to gas momentum B vector.
299           use discretelement, only: DRAG_BM
300     ! Volume of Z-momentum cell
301           use geometry, only: VOL_W
302     
303     ! Global Parameters:
304     !---------------------------------------------------------------------//
305     ! Size of IJK arrays and solids phase..
306           use param, only: DIMENSION_3, DIMENSION_M
307     ! Fluid grid loop bounds.
308           use compar, only: IJKStart3, IJKEnd3
309     ! The I, J, and K values that comprise an IJK
310           use indices, only: I_OF, J_OF, K_OF
311     ! Flag: Fluid exists at indexed cell
312           use functions, only: FLUID_AT
313     ! IJK of cell to top.
314           use functions, only: TOP_OF
315     ! IJK function for I,J,K that includes mapped indices.
316           use compar, only: FUNIJK_MAP_C
317     ! Function for averaging to a scalar cell's north face.
318           use fun_avg, only: AVG_Z
319     ! Domain index bounds.
320           use compar, only: ISTART2, JSTART2, KSTART2
321           use compar, only: IEND2, JEND2, KEND2
322     
323           IMPLICIT NONE
324     
325     ! Dummy Arguments:
326     !---------------------------------------------------------------------//
327     ! Septadiagonal matrix A_m
328           DOUBLE PRECISION, INTENT(INOUT) :: A_M(DIMENSION_3,-3:3,0:DIMENSION_M)
329     ! Vector b_m
330           DOUBLE PRECISION, INTENT(INOUT) :: B_M(DIMENSION_3, 0:DIMENSION_M)
331     ! Error index
332           INTEGER, INTENT(INOUT) :: IER
333     
334     ! Local Variables:
335     !---------------------------------------------------------------------//
336     ! Grid cell indices
337           INTEGER :: I, J, K, IJK, IMJK, IJMK, IMJMK, IJKT
338     ! temporary variables for matrix A_M and vector B_M
339           DOUBLE PRECISION tmp_A, tmp_B
340     ! Averaging factor
341     ! (=0.25 in 3D and =0.5 in 2D)
342           DOUBLE PRECISION :: AVG_FACTOR
343     !......................................................................!
344     
345     ! Initialize error flag.
346           IER = 0
347     
348     ! Skip this routine if the gas/solids are only one-way coupled.
349           IF(DES_ONEWAY_COUPLED) RETURN
350     
351           IF(DES_INTERP_SCHEME_ENUM == DES_INTERP_GARG)THEN
352     
353              AVG_FACTOR = 0.25d0
354     
355     !!$omp parallel do schedule (guided,50) default(none)                   &
356     !!$omp shared(IJKSTART3, IJKEND3, FUNIJK_MAP_C, I_OF, J_OF,             &
357     !!$omp    K_OF, AVG_FACTOR, DRAG_AM, DRAG_BM, A_M, B_M, VOL_W)          &
358     !!$omp private(IJK, I, J, K, IMJK, IJMK, IMJMK, tmp_A, tmp_B)
359              DO IJK = IJKSTART3, IJKEND3
360                 IF(.NOT.FLUID_AT(IJK)) CYCLE
361     
362                 I = I_OF(IJK)
363                 J = J_OF(IJK)
364                 K = K_OF(IJK)
365     
366                 IF (I.LT.ISTART2 .OR. I.GT.IEND2) CYCLE
367                 IF (J.LT.JSTART2 .OR. J.GT.JEND2) CYCLE
368                 IF (K.LT.KSTART2 .OR. K.GT.KEND2) CYCLE
369     
370                 IMJK = FUNIJK_MAP_C(I-1,J,K)
371                 IJMK = FUNIJK_MAP_C(I,J-1,K)
372                 IMJMK = FUNIJK_MAP_C(I-1,J-1,K)
373     
374                 tmp_A = -AVG_FACTOR*(DRAG_AM(IJK) + DRAG_AM(IMJK) +        &
375                    DRAG_AM(IJMK) + DRAG_AM(IMJMK))
376     
377                 tmp_B = -AVG_FACTOR*(DRAG_BM(IJK,3) + DRAG_BM(IMJK,3) +    &
378                    DRAG_BM(IJMK,3) + DRAG_BM(IMJMK,3))
379     
380                 A_M(IJK,0,0) = A_M(IJK,0,0) + tmp_A*VOL_W(IJK)
381                 B_M(IJK,0) = B_M(IJK,0) + tmp_B*VOL_W(IJK)
382     
383              ENDDO
384     !!$omp end parallel do
385     
386           ELSE
387     
388     !$omp parallel do default(none) &
389     !$omp private(IJK, K, IJKT, tmp_A, tmp_B) &
390     !$omp shared(IJKSTART3,IJKEND3,K_OF, F_GDS, DRAG_BM, A_M, B_M, VOL_W)
391              DO IJK = IJKSTART3, IJKEND3
392                 IF(FLUID_AT(IJK)) THEN
393                    K = K_OF(IJK)
394                    IJKT = TOP_OF(IJK)
395     
396                    tmp_A = AVG_Z(F_GDS(IJK), F_GDS(IJKT), K)
397                    tmp_B = AVG_Z(DRAG_BM(IJK,3), DRAG_BM(IJKT,3), K)
398     
399                    A_M(IJK,0,0) = A_M(IJK,0,0) - VOL_W(IJK) * tmp_A
400                    B_M(IJK,0) = B_M(IJK,0) - VOL_W(IJK) * tmp_B
401                 ENDIF
402              ENDDO
403     !$omp end parallel do
404     
405           ENDIF
406     
407           RETURN
408           END SUBROUTINE GAS_DRAG_W
409