File: RELATIVE:/../../../mfix.git/model/des/comp_mean_fields0.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
3           SUBROUTINE COMP_MEAN_FIELDS0
4     
5     !-----------------------------------------------
6     ! Modules
7     !-----------------------------------------------
8           USE param
9           USE param1
10           USE parallel
11           USE constant
12           USE physprop
13           USE fldvar
14           USE run
15           USE geometry
16           USE indices
17           USE bc
18           USE compar
19           USE sendrecv
20           USE discretelement
21           USE drag
22           USE interpolation
23           use desmpi
24           USE cutcell
25           USE mfix_pic
26           USE mpi_utility
27     
28     
29           use mpi_node_des, only: des_addnodevalues_mean_fields
30           use particle_filter, only: DES_REPORT_MASS_INTERP
31     
32     
33           use functions, only: FLUID_AT
34           use functions, only: FUNIJK
35           use functions, only: IS_ON_myPE_wobnd
36     
37           IMPLICIT NONE
38     !-----------------------------------------------
39     ! Local variables
40     !-----------------------------------------------
41     ! general i, j, k indices
42           INTEGER :: I, J, K, IJK, &
43                      II, JJ, KK
44           INTEGER :: I1, I2, J1, J2, K1, K2
45           INTEGER :: IDIM, IJK2
46           INTEGER :: CUR_IJK
47     ! indices used for interpolation stencil (unclear why IE, JN, KTP are
48     ! needed)
49           INTEGER :: IW, IE, JS, JN, KB, KTP
50     ! i,j,k indices of the fluid cell the particle resides in minus 1
51     ! (e.g., shifted 1 in west, south, bottom direction)
52           INTEGER, DIMENSION(3):: PCELL
53     ! order of interpolation set in the call to set_interpolation_scheme
54     ! unless it is re/set later through the call to
55     ! set_interpolation_stencil
56           INTEGER :: ONEW
57     ! constant whose value depends on dimension of system
58     ! avg_factor=0.250 (in 3D) or =0.50 (in 2D)
59           DOUBLE PRECISION :: AVG_FACTOR
60     ! index of solid phase that particle NP belongs to
61           INTEGER :: M
62     ! particle number index, used for looping
63           INTEGER :: NP, NINDX
64     ! Statistical weight of the particle. Equal to one for DEM
65           DOUBLE PRECISION :: WTP
66     
67           DOUBLE PRECISION :: MASS_SOL1, MASS_SOL2
68     ! sum of mass_sol1 and mass_sol2 across all processors
69           DOUBLE PRECISION :: MASS_SOL1_ALL, MASS_SOL2_ALL
70     
71           DOUBLE PRECISION :: TEMP1, TEMP2
72     
73           DOUBLE PRECISION, DIMENSION(3) :: DES_VEL_DENSITY
74           DOUBLE PRECISION :: DES_ROP_DENSITY
75     
76           INTEGER :: COUNT_NODES_OUTSIDE, COUNT_NODES_INSIDE, &
77                      COUNT_NODES_INSIDE_MAX
78           double precision :: RESID_ROPS(DES_MMAX), &
79                               RESID_VEL(3, DES_MMAX)
80           double precision :: NORM_FACTOR
81     !Handan Liu added on Jan 17 2013
82               DOUBLE PRECISION, DIMENSION(2,2,2,3) :: gst_tmp
83               DOUBLE PRECISION, DIMENSION(2,2,2) :: weight_ft
84     !-----------------------------------------------
85     
86     ! initializing
87           MASS_SOL1 = ZERO
88           MASS_SOL2 = ZERO
89           MASS_SOL1_ALL = ZERO
90           MASS_SOL2_ALL = ZERO
91     ! avg_factor=0.25 (in 3D) or =0.5 (in 2D)
92           AVG_FACTOR = merge(0.5d0, 0.25D0, NO_K)
93     
94     ! cartesian_grid related quantities
95           COUNT_NODES_INSIDE_MAX = merge(4, 8, NO_K)
96     
97     
98     ! Initialize entire arrays to zero
99           DES_VEL_NODE = ZERO
100           DES_ROPS_NODE = ZERO
101           DES_ROP_S = zero
102           DES_U_S = ZERO
103           DES_V_S = ZERO
104           IF(DO_K) DES_W_S = ZERO
105     
106     
107     ! sets several quantities including interp_scheme, scheme, and
108     ! order and allocates arrays necessary for interpolation
109           CALL SET_INTERPOLATION_SCHEME(2)
110     
111     !$omp parallel default(shared)                                             &
112     !$omp private(IJK, I, J, K, PCELL, IW, IE, JS, JN, KB, KTP, ONEW, GST_TMP, &
113     !$omp    COUNT_NODES_INSIDE, II, JJ, KK, CUR_IJK, NINDX, NP, WTP, M,       &
114     !$omp    WEIGHT_FT, I1, I2, J1, J2, K1, K2, IDIM,                          &
115     !$omp    IJK2, NORM_FACTOR, RESID_ROPS, RESID_VEL,COUNT_NODES_OUTSIDE, TEMP1)
116     !$omp do reduction(+:MASS_SOL1) reduction(+:DES_ROPS_NODE,DES_VEL_NODE)
117           DO IJK = IJKSTART3,IJKEND3
118     
119     ! Cycle this cell if not in the fluid domain or if it contains no
120     ! particle/parcel
121              IF(.NOT.FLUID_AT(IJK)) CYCLE
122              IF( PINC(IJK) == 0) CYCLE
123     
124              PCELL(1) = I_OF(IJK)-1
125              PCELL(2) = J_OF(IJK)-1
126              PCELL(3) = merge(K_OF(IJK)-1, 1, DO_K)
127     
128     ! setup the stencil based on the order of interpolation and factoring in
129     ! whether the system has any periodic boundaries. sets onew to order.
130              CALL SET_INTERPOLATION_STENCIL(PCELL, IW, IE, JS, JN, KB, KTP,&
131                 INTERP_SCHEME, DIMN, ORDERNEW=ONEW)
132     
133              COUNT_NODES_OUTSIDE = 0
134     ! Computing/setting the geometric stencil
135              DO K=1, merge(1, ONEW, NO_K)
136              DO J=1, ONEW
137              DO I=1, ONEW
138     
139                 II = IW + I-1
140                 JJ = JS + J-1
141                 KK = KB + K-1
142                 CUR_IJK = funijk_map_c(II,JJ,KK)
143     
144                 GST_TMP(I,J,K,1) = XE(II)
145                 GST_TMP(I,J,K,2) = YN(JJ)
146                 GST_TMP(I,J,K,3) = merge(DZ(1), ZT(KK), NO_K)
147     
148                 IF(CARTESIAN_GRID) THEN
149                    IF(SCALAR_NODE_ATWALL(CUR_IJK))                   &
150                       COUNT_NODES_OUTSIDE = COUNT_NODES_OUTSIDE + 1
151                 ENDIF
152     
153              ENDDO
154              ENDDO
155              ENDDO
156     
157     
158     ! Calculate des_rops_node so des_rop_s, and in turn, ep_g can be updated
159     !----------------------------------------------------------------->>>
160     
161     ! looping through particles in the cell
162              DO NINDX=1, PINC(IJK)
163                 NP = PIC(IJK)%P(NINDX)
164     
165                 call DRAG_WEIGHTFACTOR(gst_tmp,des_pos_new(:,np),weight_ft)
166     
167                 M = PIJK(NP,5)
168                 WTP = ONE
169                 IF(MPPIC) WTP = DES_STAT_WT(NP)
170     
171                 MASS_SOL1 = MASS_SOL1 + PMASS(NP)*WTP
172     
173                 TEMP2 = DES_RO_S(M)*PVOL(NP)*WTP
174     
175                 DO K = 1, merge(1, ONEW, NO_K)
176                 DO J = 1, ONEW
177                 DO I = 1, ONEW
178     ! shift loop index to new variables for manipulation
179                    II = IW + I-1
180                    JJ = JS + J-1
181                    KK = KB + K-1
182     
183                    CUR_IJK = FUNIJK_MAP_C(II,JJ,KK)
184     
185                    TEMP1 = WEIGHT_FT(I,J,K)*TEMP2
186     
187                    DES_ROPS_NODE(CUR_IJK,M) = DES_ROPS_NODE(CUR_IJK,M) +   &
188                       TEMP1
189     
190                    DES_VEL_NODE(CUR_IJK,:,M) = DES_VEL_NODE(CUR_IJK,:,M) + &
191                       TEMP1*DES_VEL_NEW(:,NP)
192                 ENDDO
193                 ENDDO
194                 ENDDO
195              ENDDO   ! end do (nindx=1,pinc(ijk))
196     !-----------------------------------------------------------------<<<
197     
198     
199     ! Only for cutcell cases may count_nodes_inside become less than its
200     ! original set value. In such an event, the contribution of scalar nodes
201     ! that do not reside in the domain is added to a residual array. This
202     ! array is then redistribited equally to the nodes that are in the fluid
203     ! domain. These steps are done to conserve mass.
204     !----------------------------------------------------------------->>>
205              IF (CARTESIAN_GRID) THEN
206     
207     ! only for cartesian_grid will count_nodes_outside be modified from zero
208                 COUNT_NODES_INSIDE = &
209                    COUNT_NODES_INSIDE_MAX - COUNT_NODES_OUTSIDE
210     
211                 IF(COUNT_NODES_INSIDE.LT.COUNT_NODES_INSIDE_MAX) THEN
212     
213     ! initializing
214                    RESID_ROPS(1:DES_MMAX) = ZERO
215                    RESID_VEL(:, 1:DES_MMAX) = ZERO
216     
217     ! Convention used to number node numbers
218     ! i=1, j=2           i=2, j=2
219     !   _____________________
220     !   |                   |
221     !   |  I = 2, J = 2     |
222     !   |___________________|
223     ! i=1, j=1           i=2, j=1
224     ! setting indices based on convention
225                    I = I_OF(IJK)
226                    J = J_OF(IJK)
227                    K = K_OF(IJK)
228                    I1 = I-1
229                    I2 = I
230                    J1 = J-1
231                    J2 = J
232                    K1 = merge(K, K-1, NO_K)
233                    K2 = K
234     ! first calculate the residual des_rops_node and des_vel_node that was
235     ! computed on nodes that do not belong to the domain
236     
237                    DO KK = K1, K2
238                    DO JJ = J1, J2
239                    DO II = I1, I2
240     
241                       IJK2 = funijk(II, JJ, KK)
242     
243                       IF(SCALAR_NODE_ATWALL(IJK2)) THEN
244                          RESID_ROPS(1:DES_MMAX) = RESID_ROPS(1:DES_MMAX) + &
245                             DES_ROPS_NODE(IJK2,1:DES_MMAX)
246     
247                          DES_ROPS_NODE(IJK2,1:DES_MMAX) = ZERO
248                          DO IDIM = 1, merge(2,3,NO_K)
249                             RESID_VEL(IDIM, 1:DES_MMAX) =                  &
250                                 RESID_VEL(IDIM, 1:DES_MMAX) +              &
251                                DES_VEL_NODE(IJK2,IDIM, 1:DES_MMAX)
252                             DES_VEL_NODE(IJK2,IDIM, 1:DES_MMAX) = ZERO
253                          ENDDO
254                       ENDIF
255                    ENDDO
256                    ENDDO
257                    ENDDO
258     
259     ! now add this residual equally to the remaining nodes
260                    NORM_FACTOR = ONE/REAL(COUNT_NODES_INSIDE)
261                    DO KK = K1, K2
262                    DO JJ = J1, J2
263                    DO II = I1, I2
264                       IJK2 = funijk(II, JJ, KK)
265     
266                       IF(.NOT.SCALAR_NODE_ATWALL(IJK2)) THEN
267                          DES_ROPS_NODE(IJK2,1:DES_MMAX) =                  &
268                             DES_ROPS_NODE(IJK2,1:DES_MMAX) +               &
269                             RESID_ROPS(1:DES_MMAX)*NORM_FACTOR
270                          DO IDIM = 1, merge(2,3,NO_K)
271                             DES_VEL_NODE(IJK2,IDIM, 1:DES_MMAX) =          &
272                                DES_VEL_NODE(IJK2,IDIM, 1:DES_MMAX) +       &
273                                RESID_VEL(IDIM, 1:DES_MMAX)*NORM_FACTOR
274                          ENDDO
275                       ENDIF
276     
277                    ENDDO
278                    ENDDO
279                    ENDDO
280                 ENDIF
281              ENDIF   ! end if (cartesian_grid)
282           ENDDO
283     !$omp end parallel
284     
285     
286     ! At the interface des_rops_node has to be added since particles
287     ! across the processors will contribute to the same scalar node.
288     ! sendrecv will be called and the node values will be added
289     ! at the junction. des_rops_node is altered by the routine when
290     ! periodic boundaries are invoked
291           CALL DES_ADDNODEVALUES_MEAN_FIELDS
292     
293     
294     ! Now go from node to scalar center. Same convention as sketched
295     ! earlier
296     !----------------------------------------------------------------->>>
297     ! Explanation by RG: 08/17/2012
298     ! the approach used here is to make it general enough for cutcells to be
299     ! included as well. The new changes do not alter earlier calculations
300     ! but make the technique general as to include cartesian grid (cut-cell)
301     ! simulations.
302     ! Previously, the volume of the node (by array des_vol_node) was used to
303     ! first scale the nodal values. Subsequently, these nodal values were
304     ! equally added to compute the cell centered values for the scalar cell.
305     
306     ! Consider an internal node next to an edge node (a node adjacent to a
307     ! boundary). In 2D, the volume of an edge node will be half that of an
308     ! internal node. And an edge node will contribute double compared to
309     ! an internal node to the value of the scalar cell they share. These
310     ! calculations were previously accomplished via the variable volume of
311     ! node.  Now this is accomplished by the ratio vol(ijk2)/vol_sur, where
312     ! vol(ijk2) is the volume of the scalar cell in consideration and
313     ! vol_sur is the sum of all the scalar cell volumes that have this node
314     ! as the common node.
315     !---------------------------------------------------------------------//
316     !$omp parallel do default(none) collapse (3)                           &
317     !$omp shared(KSTART2, KEND1, JSTART2, JEND1, ISTART2, IEND1, DO_K, VOL,&
318     !$omp   DEAD_CELL_AT, FUNIJK_MAP_C, VOL_SURR, DES_MMAX, DES_ROPS_NODE, &
319     !$omp   DES_VEL_NODE)                                                  &
320     !$omp private(I, J, K, IJK, M, II, JJ, KK, IJK2, DES_ROP_DENSITY,      &
321     !$omp   DES_VEL_DENSITY)                                               &
322     !$omp reduction(+:DES_ROP_S, DES_U_S, DES_V_S, DES_W_S)
323           DO K = KSTART2, KEND1
324           DO J = JSTART2, JEND1
325           DO I = ISTART2, IEND1
326              IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
327              IJK = funijk(I,J,K)
328              if (VOL_SURR(IJK).eq.ZERO) CYCLE ! no FLUID_AT any of the stencil points have
329     
330     ! looping over stencil points (NODE VALUES)
331              DO M = 1, DES_MMAX
332     
333                 DES_ROP_DENSITY = DES_ROPS_NODE(IJK, M)/VOL_SURR(IJK)
334                 DES_VEL_DENSITY(:) = DES_VEL_NODE(IJK, :, M)/VOL_SURR(IJK)
335     
336                 DO KK = K, merge(K+1, K, DO_K)
337                 DO JJ = J, J+1
338                 DO II = I, I+1
339                    IF (DEAD_CELL_AT(II,JJ,KK)) CYCLE  ! skip dead cells
340     
341                    IJK2 = funijk_map_c(II, JJ, KK)
342                    IF(FLUID_AT(IJK2).and.(IS_ON_myPE_wobnd(II, JJ, KK))) THEN
343     ! Since the data in the ghost cells is spurious anyway and overwritten during
344     ! subsequent send receives, do not compute any value here as this will
345     ! mess up the total mass value that is computed below to ensure mass conservation
346     ! between Lagrangian and continuum representations
347                       DES_ROP_S(IJK2, M) = DES_ROP_S(IJK2, M) + DES_ROP_DENSITY*VOL(IJK2)
348                       DES_U_S(IJK2, M) = DES_U_S(IJK2, M) + DES_VEL_DENSITY(1)*VOL(IJK2)
349                       DES_V_S(IJK2, M) = DES_V_S(IJK2, M) + DES_VEL_DENSITY(2)*VOL(IJK2)
350                       IF(DO_K) DES_W_S(IJK2, M) = DES_W_S(IJK2, M) + DES_VEL_DENSITY(3)*VOL(IJK2)
351                    ENDIF
352                 ENDDO  ! end do (ii=i1,i2)
353                 ENDDO  ! end do (jj=j1,j2)
354                 ENDDO  ! end do (kk=k1,k2)
355              ENDDO
356           ENDDO   ! end do (i=istart2,iend1)
357           ENDDO   ! end do (j=jstart2,jend1)
358           ENDDO   ! end do (k=kstart2,kend1)
359     !omp end parallel do
360     
361     !-----------------------------------------------------------------<<<
362     
363     
364     !$omp parallel do default(none) private(IJK, M)                        &
365     !$omp shared(IJKSTART3, IJKEND3, DO_K, DES_MMAX, DES_ROP_s, DES_U_S,   &
366     !$omp   DES_V_S, DES_W_S, VOL)
367           DO IJK = IJKSTART3, IJKEND3
368              IF(.NOT.FLUID_AT(IJK)) CYCLE
369     
370              DO M = 1, DES_MMAX
371                 IF(DES_ROP_S(IJK, M).GT.ZERO) THEN
372                    DES_U_S(IJK, M) = DES_U_S(IJK,M)/DES_ROP_S(IJK, M)
373                    DES_V_S(IJK, M) = DES_V_S(IJK,M)/DES_ROP_S(IJK, M)
374                    IF(DO_K) DES_W_S(IJK, M) = DES_W_S(IJK,M)/DES_ROP_S(IJK, M)
375     
376     ! Divide by scalar cell volume to obtain the bulk density
377                    DES_ROP_S(IJK, M) = DES_ROP_S(IJK, M)/VOL(IJK)
378     
379                 ENDIF
380              ENDDO   ! end loop over M=1,DES_MMAX
381           ENDDO  ! end loop over IJK=ijkstart3,ijkend3
382     !omp end parallel do
383     
384     
385           IF (MPPIC) CALL SEND_RECV(DES_ROP_S,2)
386     
387           CALL CALC_EPG_DES
388     
389     
390           IF(MPPIC) THEN
391     
392     ! Now calculate Eulerian mean velocity fields like U_S, V_S, and W_S.
393              CALL SEND_RECV(DES_U_S,2)
394              CALL SEND_RECV(DES_V_S,2)
395              IF(DO_K) CALL SEND_RECV(DES_W_S,2)
396     
397     ! The Eulerian velocity field is used to set up the stencil to interpolate
398     ! mean solid velocity at the parcel's location. DES_U_S could have also been
399     ! used, but that also would have require the communication at this stage.
400     ! The final interpolated value does not change if the stencil is formed by
401     ! first obtaining face centered Eulerian velocities (U_S, etc.)
402     ! and then computing the node velocities from them or directly computing
403     ! the node velocities from cell centered average velocity field (DES_U_S,
404     ! etc.). We are using the first approach as it is more natural to set
405     ! BC's on solid velocity field in the face centered represenation (U_S,
406     ! etc.)
407     
408              IF(.NOT.CARTESIAN_GRID) THEN
409                 CALL MPPIC_COMP_EULERIAN_VELS_NON_CG
410              ELSE
411                 CALL MPPIC_COMP_EULERIAN_VELS_CG
412              ENDIF
413           ENDIF   ! end if (.not.mppic)
414     
415     ! turn on the below statements to check if the mass is conserved
416     ! between discrete and continuum representations. Should be turned to
417     ! false for any production runs.
418           IF(DES_REPORT_MASS_INTERP) THEN
419     
420     
421              DO IJK = IJKSTART3, IJKEND3
422                 IF(.NOT.FLUID_AT(IJK)) CYCLE
423     
424                 I = I_OF(IJK)
425                 J = J_OF(IJK)
426                 K = K_OF(IJK)
427     
428     ! It is important to check both FLUID_AT and IS_ON_MYPE_WOBND.
429                 IF(IS_ON_myPE_wobnd(I,J,K)) MASS_SOL2 = MASS_SOL2 +        &
430                    sum(DES_ROP_S(IJK,1:DES_MMAX))*VOL(IJK)
431              ENDDO
432     
433     
434              CALL GLOBAL_SUM(MASS_SOL1, MASS_SOL1_ALL)
435              CALL GLOBAL_SUM(MASS_SOL2, MASS_SOL2_ALL)
436              if(myPE.eq.pe_IO) THEN
437                 WRITE(*,'(/,5x,A,4(2x,g17.8),/)') &
438                      'SOLIDS MASS DISCRETE AND CONTINUUM =  ', &
439                      MASS_SOL1_ALL, MASS_SOL2_ALL
440              ENDIF
441           ENDIF
442           END SUBROUTINE COMP_MEAN_FIELDS0
443     
444     
445     
446     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
447     !  Subroutine: DRAG_WEIGHTFACTOR                                        C
448     !  Purpose: DES - Calculate the fluid velocity interpolated at the      C
449     !           particle's location and weights. Replace 'interpolator'     C
450     !                       interface for OpenMP implementation.            C
451     !                                                                       C
452     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
453     
454           SUBROUTINE DRAG_WEIGHTFACTOR(GSTEN,DESPOS,WEIGHTFACTOR)
455     
456           use geometry, only: NO_K
457     
458             IMPLICIT NONE
459     
460     !-----------------------------------------------
461     ! Local Variables
462     !-----------------------------------------------
463             DOUBLE PRECISION, DIMENSION(2,2,2,3), INTENT(IN):: GSTEN
464             DOUBLE PRECISION, DIMENSION(3), INTENT(IN):: DESPOS
465             DOUBLE PRECISION, DIMENSION(2,2,2), INTENT(OUT) :: WEIGHTFACTOR
466             INTEGER :: II, JJ, KK
467     
468             DOUBLE PRECISION, DIMENSION(2) :: XXVAL, YYVAL, ZZVAL
469             DOUBLE PRECISION :: DXX, DYY, DZZ
470             DOUBLE PRECISION, DIMENSION(3) :: ZETAA
471     
472             DXX = GSTEN(2,1,1,1) - GSTEN(1,1,1,1)
473             DYY = GSTEN(1,2,1,2) - GSTEN(1,1,1,2)
474     
475             ZETAA(1:2) = DESPOS(1:2) - GSTEN(1,1,1,1:2)
476     
477             ZETAA(1) = ZETAA(1)/DXX
478             ZETAA(2) = ZETAA(2)/DYY
479     
480             XXVAL(1)=1-ZETAA(1)
481             YYVAL(1)=1-ZETAA(2)
482             XXVAL(2)=ZETAA(1)
483             YYVAL(2)=ZETAA(2)
484     
485             IF(NO_K) THEN
486                DO JJ=1,2
487                   DO II=1,2
488                      WEIGHTFACTOR(II,JJ,1) = XXVAL(II)*YYVAL(JJ)
489                   ENDDO
490                ENDDO
491             ELSE
492                DZZ = GSTEN(1,1,2,3) - GSTEN(1,1,1,3)
493                ZETAA(3) = DESPOS(3) - GSTEN(1,1,1,3)
494                ZETAA(3) = ZETAA(3)/DZZ
495                ZZVAL(1)=1-ZETAA(3)
496                ZZVAL(2)=ZETAA(3)
497                DO KK=1,2
498                   DO JJ=1,2
499                      DO II=1,2
500                         WEIGHTFACTOR(II,JJ,KK) = XXVAL(II)*YYVAL(JJ)*ZZVAL(KK)
501                      ENDDO
502                   ENDDO
503                ENDDO
504             ENDIF
505     
506           END SUBROUTINE DRAG_WEIGHTFACTOR
507