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