File: N:\mfix\model\des\pic\interpolate_pic.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !  Subroutine: CALC_PS_GRAD_PIC                                        !
3     !  Author: R. Garg                                                     !
4     !                                                                      !
5     !  Purpose:                                                            !
6     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
7           SUBROUTINE INTERPOLATE_PIC
8     
9           use particle_filter, only: DES_INTERP_SCHEME_ENUM
10           use particle_filter, only: DES_INTERP_GARG
11     
12           IMPLICIT NONE
13     
14           SELECT CASE(DES_INTERP_SCHEME_ENUM)
15           CASE(DES_INTERP_GARG) ; CALL INTERPOLATE_PIC0
16           CASE DEFAULT; CALL INTERPOLATE_PIC1
17           END SELECT
18     
19           RETURN
20           END SUBROUTINE INTERPOLATE_PIC
21     
22     
23     
24     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
25     !  Subroutine: INTERPOLATE_PIC_NONE                                    !
26     !  Author: R. Garg                                                     !
27     !                                                                      !
28     !  Purpose:                                                            !
29     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
30           SUBROUTINE INTERPOLATE_PIC1
31     
32           use discretelement, only: MAX_PIP, PIJK
33           USE geometry, only: DO_K
34           use fldvar, only: EP_G, u_s, v_s, w_s
35           use functions, only: FLUID_AT
36           use functions, only: IS_NONEXISTENT
37           use mfix_pic, only: AVGVEL => AVGSOLVEL_P
38           use mfix_pic, only: PS_GRAD, PS_FORCE_PIC
39           use mfix_pic, only: EPG_P
40     ! Flag to use interpolation
41           use particle_filter, only: DES_INTERP_ON
42     ! Interpolation cells and weights
43           use particle_filter, only: FILTER_CELL, FILTER_WEIGHT
44           use physprop, only: mmax
45     
46           implicit none
47     
48     ! Local Variables
49     !---------------------------------------------------------------------//
50     ! Interpolation weight
51           DOUBLE PRECISION :: WEIGHT
52     ! temporary variables used to calculate pressure at scalar cell edge
53           INTEGER :: IJK, NP, LC
54     ! Loop bound for filter
55           INTEGER :: LP_BND
56     !......................................................................!
57     
58     ! Loop bounds for interpolation.
59           LP_BND = merge(27,9,DO_K)
60     
61     
62           DO NP=1,MAX_PIP
63     
64     ! Skip parcels that don't exist or are in a non-fluid cell.
65              IF(IS_NONEXISTENT(NP) .OR. .NOT.FLUID_AT(PIJK(NP,4))) THEN
66                 AVGVEL(:,NP) = 0.0d0
67                 PS_GRAD(:,NP) = 0.0d0
68                 EPG_P(NP) = 1.0d0
69     
70     ! Interpolated solids pressure and average solids velocity.
71              ELSEIF(DES_INTERP_ON) THEN
72                 AVGVEL(:,NP) = 0.0d0
73                 PS_GRAD(:,NP) = 0.0d0
74                 EPG_P(NP) = 0.0d0
75                 DO LC=1,LP_BND
76                    IJK = FILTER_CELL(LC,NP)
77                    WEIGHT = FILTER_WEIGHT(LC,NP)
78     ! Gas phase velocity.
79                    AVGVEL(1,NP) = AVGVEL(1,NP) + U_s(IJK,MMAX+1)*WEIGHT
80                    AVGVEL(2,NP) = AVGVEL(2,NP) + V_s(IJK,MMAX+1)*WEIGHT
81                    AVGVEL(3,NP) = AVGVEL(3,NP) + W_s(IJK,MMAX+1)*WEIGHT
82     ! Particle normal stress.
83                    PS_GRAD(:,NP) = PS_GRAD(:,NP)+PS_FORCE_PIC(:,IJK)*WEIGHT
84     ! TO BE REMOVED
85                    EPG_P(NP) = EPG_P(NP) + EP_G(IJK)*WEIGHT
86                 ENDDO
87     
88     ! Centroid method. 
89             ELSE
90                 IJK = PIJK(NP,4)
91     ! Gas phase velocity.
92                 AVGVEL(1,NP) = U_s(IJK,MMAX+1)
93                 AVGVEL(2,NP) = V_s(IJK,MMAX+1)
94                 AVGVEL(3,NP) = W_s(IJK,MMAX+1)
95     ! Particle normal stress.
96                 PS_GRAD(:,NP) = PS_FORCE_PIC(:,IJK)
97     ! TO BE REMOVED
98                 EPG_P(NP) = EP_G(IJK)
99              ENDIF
100           ENDDO
101     
102           RETURN
103           END SUBROUTINE INTERPOLATE_PIC1
104     
105     
106     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
107     !  Subroutine: CALC_PS_GRAD_PIC                                        !
108     !  Author: R. Garg                                                     !
109     !                                                                      !
110     !  Purpose:                                                            !
111     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
112           SUBROUTINE INTERPOLATE_PIC0
113     
114           USE bc
115           USE compar
116           USE constant
117           USE cutcell
118           USE derived_types, only: pic
119           USE discretelement
120           USE fldvar, only: ep_g, u_g, v_g, w_g
121           USE fun_avg
122           USE functions
123           USE geometry
124           USE indices
125           USE interpolation
126           USE mfix_pic
127           USE parallel
128           USE param
129           USE param1
130           USE physprop
131           USE run
132           USE sendrecv
133     
134           implicit none
135     
136           ! general i, j, k indices
137           INTEGER I, J, K, IJK, IPJK, IJPK, IJKP, IDIM, M
138     
139           ! temporary variables used to calculate pressure at scalar cell edge
140           DOUBLE PRECISION avg_factor, VOL_TOT_VEC(3), VOL_TOT_SCAL
141     
142           integer :: korder, iw,ie,js,jn,kb,ktp, onew, pcell(3), cur_ijk, NP, nindx
143     
144           integer :: ii,jj,kk, ipjpk, ijpkp, ipjkp, ipjpkp
145     
146           double precision :: vol_ijk, vol_ipjk, vol_ijpk, vol_ipjpk
147           double precision :: vol_ijkp, vol_ipjkp, vol_ijpkp, vol_ipjpkp
148     
149     
150           CALL SET_INTERPOLATION_SCHEME(2)
151     
152           KORDER = merge ( 1, 2, no_k) !1+(DIMN-2)
153     
154     ! avg_factor=0.25 (in 3D) or =0.5 (in 2D)
155           !avg_factor = 0.25d0*(dimn-2) + 0.5d0*(3-dimn)
156           AVG_FACTOR = merge(0.5d0, 0.25D0, NO_K)
157     
158           do ijk = ijkstart3,ijkend3
159     
160              if(.not.fluid_at(ijk) .or. pinc(ijk).eq.0) cycle
161              i = i_of(ijk)
162              j = j_of(ijk)
163              k = k_of(ijk)
164     
165              pcell(1) = i-1
166              pcell(2) = j-1
167              pcell(3) = merge(1, k-1, no_k) ! =k-1 (in 3d) or =1 (in 2d)
168     
169              CALL SET_INTERPOLATION_STENCIL(PCELL,IW,IE,JS,JN,KB,          &
170                 KTP,INTERP_SCHEME,MERGE(2,3,NO_K),ORDERNEW = ONEW)
171     
172     ! Compute velocity at grid nodes and set the geometric stencil
173              DO K = 1, merge(1, ONEW, NO_K)!  (3-DIMN)*1+(DIMN-2)*ONEW
174              DO J = 1,ONEW
175              DO I = 1,ONEW
176     
177                 II = IW + I-1
178                 JJ = JS + J-1
179                 KK = KB + K-1
180     
181                 CUR_IJK = FUNIJK_MAP_C(II,JJ,KK)
182     
183                 IPJK    = FUNIJK_MAP_C(II+1,JJ,KK)
184                 IJPK    = FUNIJK_MAP_C(II,JJ+1,KK)
185                 IPJPK   = FUNIJK_MAP_C(II+1,JJ+1,KK)
186     
187                 VOL_IJK = ZERO
188                 VOL_IPJK = ZERO
189                 VOL_IJPK = ZERO
190                 VOL_IPJPK = ZERO
191     
192                 VOL_IJKP = ZERO
193                 VOL_IPJKP = ZERO
194                 VOL_IJPKP = ZERO
195                 VOL_IPJPKP = ZERO
196     
197                 IF(FLUID_AT(CUR_IJK)) VOL_IJK   = VOL(CUR_IJK)
198                 IF(FLUID_AT(IPJK))    VOL_IPJK  = VOL(IPJK)
199                 IF(FLUID_AT(IJPK))    VOL_IJPK  = VOL(IJPK)
200                 IF(FLUID_AT(IPJPK))   VOL_IPJPK = VOL(IPJPK)
201     
202                 IF(DO_K) THEN
203                    IJKP   = FUNIJK_MAP_C(II,JJ,KK+1)
204                    IJPKP  = FUNIJK_MAP_C(II,JJ+1,KK+1)
205                    IPJKP  = FUNIJK_MAP_C(II+1,JJ,KK+1)
206                    IPJPKP = FUNIJK_MAP_C(II+1,JJ+1,KK+1)
207     
208                    IF(FLUID_AT(IJKP))   VOL_IJKP   = VOL(IJKP)
209                    IF(FLUID_AT(IPJKP))  VOL_IPJKP  = VOL(IPJKP)
210                    IF(FLUID_AT(IJPKP))  VOL_IJPKP  = VOL(IJPKP)
211                    IF(FLUID_AT(IPJPKP)) VOL_IPJPKP = VOL(IPJPKP)
212                 ENDIF
213     
214                 GSTENCIL(I,J,K,1) = XE(II)
215                 GSTENCIL(I,J,K,2) = YN(JJ)
216                 GSTENCIL(I,J,K,3) = merge(DZ(1), ZT(KK), NO_K)
217     
218                 VOL_TOT_SCAL = ZERO
219     
220                 VOL_TOT_SCAL = VOL_IJK + VOL_IPJK + VOL_IJPK + VOL_IPJPK + &
221                    VOL_IJKP + VOL_IPJKP + VOL_IJPKP + VOL_IPJPKP
222     
223                 VOL_TOT_VEC = ZERO
224     
225                 VOL_TOT_VEC(1) = VOL(CUR_IJK) + VOL(IJPK)
226                 VOL_TOT_VEC(2) = VOL(CUR_IJK) + VOL(IPJK)
227     
228                 DO M = MMAX+1,DES_MMAX+MMAX
229                    VEL_SOL_STENCIL(I,J,K,1,M) = &
230                       PIC_U_S(CUR_IJK,M)*VOL(CUR_IJK) + &
231                       PIC_U_S(IJPK,M)*VOL(IJPK)
232     
233                    VEL_SOL_STENCIL(I,J,K,2,M) = &
234                       PIC_V_S(CUR_IJK,M)*VOL(CUR_IJK) + &
235                       PIC_V_S(IPJK,M)*VOL(IPJK)
236                 ENDDO
237     
238                 SSTENCIL(I,J,K) = EP_G(CUR_IJK)*VOL_IJK + &
239                    EP_G(IPJK)*VOL_IPJK + EP_G(IJPK)*VOL_IJPK + &
240                    EP_G(IPJPK)*VOL_IPJPK
241     
242                 PSGRADSTENCIL(I,J,K,1) = &
243                    PS_FORCE_PIC(1,CUR_IJK)*VOL(CUR_IJK) + &
244                    PS_FORCE_PIC(1,IJPK)*VOL(IJPK)
245     
246                 PSGRADSTENCIL(I,J,K,2) = &
247                    PS_FORCE_PIC(2,CUR_IJK)*VOL(CUR_IJK) + &
248                    PS_FORCE_PIC(2,IPJK)*VOL(IPJK)
249     
250                 VSTENCIL(I,J,K,1) = &
251                    U_G(CUR_IJK)*VOL(CUR_IJK) + U_G(IJPK)*VOL(IJPK)
252     
253                 VSTENCIL(I,J,K,2) = &
254                    V_G(CUR_IJK)*VOL(CUR_IJK) + V_G(IPJK)*VOL(IPJK)
255     
256                 IF(DO_K) THEN
257                    VOL_TOT_VEC(1) = VOL_TOT_VEC(1) + VOL(IJKP) + VOL(IJPKP)
258                    VOL_TOT_VEC(2) = VOL_TOT_VEC(2) + VOL(IJKP) + VOL(IPJKP)
259                    VOL_TOT_VEC(3) = &
260                       VOL(CUR_IJK) + VOL(IPJK) + VOL(IJPK) + VOL(IPJPK)
261     
262                    SSTENCIL(I,J,K) = SSTENCIL(I,J,K) + &
263                       EP_G(IJKP)*VOL_IJKP + EP_G(IPJKP)*VOL_IPJKP + &
264                       EP_G(IJPKP)*VOL_IJPKP + EP_G(IPJPKP)*VOL_IPJPKP
265     
266                    PSGRADSTENCIL(I,J,K,1) = PSGRADSTENCIL(I,J,K,1) + &
267                       PS_FORCE_PIC(1,IJKP)*VOL(IJKP) + &
268                       PS_FORCE_PIC(1,IJPKP)*VOL(IJPKP)
269     
270                    PSGRADSTENCIL(I,J,K,2) = PSGRADSTENCIL(I,J,K,2) + &
271                       PS_FORCE_PIC(2,IJKP)*VOL(IJKP) + &
272                       PS_FORCE_PIC(2,IPJKP)*VOL(IPJKP)
273     
274                    PSGRADSTENCIL(I,J,K,3) = &
275                       PS_FORCE_PIC(3,CUR_IJK)*VOL(CUR_IJK)+ &
276                       PS_FORCE_PIC(3,IJPK)*VOL(IJPK)+       &
277                       PS_FORCE_PIC(3,IPJK)*VOL(IPJK)+       &
278                       PS_FORCE_PIC(3,IPJPK)*VOL(IPJPK)
279     
280                    VSTENCIL(I,J,K,1) = VSTENCIL(I,J,K,1) + &
281                       U_G(IJKP)*VOL(IJKP) + U_G(IJPKP)*VOL(IJPKP)
282     
283                    VSTENCIL(I,J,K,2) = VSTENCIL(I,J,K,2) + &
284                       V_G(IJKP)*VOL(IJKP) + V_G(IPJKP)*VOL(IPJKP)
285     
286                    VSTENCIL(I,J,K,3) = W_G(CUR_IJK)*VOL(CUR_IJK)+&
287                       W_G(IJPK)*VOL(IJPK) + W_G(IPJK)*VOL(IPJK) + &
288                       W_G(IPJPK)*VOL(IPJPK)
289     
290                    DO M = MMAX+1, DES_MMAX+MMAX
291                       VEL_SOL_STENCIL(I,J,K,1,M) = &
292                          VEL_SOL_STENCIL(I,J,K,1,M) + &
293                          PIC_U_S(IJKP,M)*VOL(IJKP) + & 
294                          PIC_U_S(IJPKP,M)*VOL(IJPKP)
295     
296                       VEL_SOL_STENCIL(I,J,K,2, M) = &
297                          VEL_SOL_STENCIL(I,J,K,2,M) + &
298                          PIC_V_S(IJKP,M)*VOL(IJKP) +&
299                          PIC_V_S(IPJKP,M)*VOL(IPJKP)
300     
301                       VEL_SOL_STENCIL(I,J,K,3, M) = &
302                          PIC_W_S(CUR_IJK,M)*VOL(CUR_IJK) + &
303                          PIC_W_S(IJPK,M)*VOL(IJPK) + &
304                          PIC_W_S(IPJK,M)*VOL(IPJK) + &
305                          PIC_W_S(IPJPK,M)*VOL(IPJPK)
306                    ENDDO
307                 ELSE
308                    PSGRADSTENCIL(I,J,K,3) = 0.D0
309                    VEL_SOL_STENCIL(I,J,K,3, MMAX+1:DES_MMAX+MMAX) = 0.D0
310                    VSTENCIL(I,J,K,3) = 0.D0
311                 ENDIF
312     
313     
314     
315     
316                 DO IDIM = 1, merge(2,3,NO_K)
317                    IF(VOL_TOT_VEC(IDIM).GT.ZERO)  THEN
318                       psgradstencil(i,j,k,idim) = &
319                          psgradstencil(i,j,k,idim)/VOL_TOT_VEC(idim)
320     
321                       VEL_SOL_STENCIL(i,j,k,idim, MMAX+1:DES_MMAX+MMAX) = &
322                          VEL_SOL_STENCIL(i,j,k,idim, MMAX+1:DES_MMAX+MMAX)/&
323                          VOL_TOT_VEC(idim)
324     
325                       vstencil(i,j,k,idim) = &
326                          vstencil(i,j,k,idim)/VOL_TOT_VEC(idim)
327                    ENDIF
328                 ENDDO
329     
330     
331                 IF(VOL_TOT_SCAL.GT.ZERO) &
332                    SSTENCIL(I,J,K) = SSTENCIL(I,J,K)/VOL_TOT_SCAL
333     
334              enddo
335              enddo
336              enddo
337     
338     
339     
340     ! Loop through particles in the cell
341              do nindx = 1,pinc(ijk)
342                 np = pic(ijk)%p(nindx)
343                 m = pijk(np,5)
344     
345     ! Interpolate particle stress to particle position
346                 if(NO_K) then !2-D
347                    call interpolator(gstencil(1:onew,1:onew,1,1:dimn), &
348                    psgradstencil(1:onew,1:onew,1,1:dimn), &
349                    des_pos_new(np,1:dimn),PS_GRAD(1:dimn,np),  &
350                    onew,interp_scheme,weightp)
351                 else !3-D, diff in psgradstencil size
352                    call interpolator(gstencil(1:onew,1:onew,1:onew,1:dimn), &
353                    psgradstencil(1:onew,1:onew,1:onew,1:dimn), &
354                    des_pos_new(np,1:dimn),PS_GRAD(1:dimn,np),  &
355                    onew,interp_scheme,weightp)
356                 endif
357     
358                 do idim = 1,  merge(2,3,NO_K)
359                    AVGSOLVEL_P(IDIM,NP) = ARRAY_DOT_PRODUCT(               &
360                       VEL_SOL_STENCIL(:,:,:,IDIM,M),WEIGHTP(:,:,:))
361                 ENDDO
362     
363                 EPG_P(NP) = ARRAY_DOT_PRODUCT(SSTENCIL(:,:,:),WEIGHTP(:,:,:))
364     
365              ENDDO
366           ENDDO
367     
368     
369           END SUBROUTINE INTERPOLATE_PIC0
370