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

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