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

1           MODULE GENERATE_PARTICLES
2     
3           CONTAINS
4     
5     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
6     !                                                                      C
7     !  SUBROUTINE: GENERATE_PARTICLE_CONFIG                                C
8     !                                                                      C
9     !  Purpose: Generate particle configuration based on maximum particle  C
10     !           radius and filling from top to bottom within specified     C
11     !           bounds                                                     C
12     !                                                                      C
13     !                                                                      C
14     !  Authors: Rahul Garg                                Date: 19-Mar-14  C
15     !                                                                      C
16     !                                                                      C
17     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
18           SUBROUTINE GENERATE_PARTICLE_CONFIG
19     
20           use mfix_pic, only: MPPIC
21           use discretelement, only: PIP, PARTICLES
22     ! Flag indicating that the IC region is defined.
23           use ic, only: IC_DEFINED
24     ! Parameter for detecting unspecified values, zero, and one
25           use param1, only: ONE
26     ! Maximum number of initial conditions
27           use param, only: DIMENSION_IC
28     ! IC Region gas volume fraction.
29           use ic, only: IC_EP_G
30     
31           use mpi_utility
32     
33     ! Use the error manager for posting error messages.
34     !---------------------------------------------------------------------//
35           use error_manager
36     
37           IMPLICIT NONE
38     
39           INTEGER :: ICV
40     
41     ! Initialize the error manager.
42           CALL INIT_ERR_MSG("Generate_Particle_Config")
43     
44           DO ICV = 1, DIMENSION_IC
45     
46              IF(.NOT.IC_DEFINED(ICV)) CYCLE
47              IF(IC_EP_G(ICV) == ONE) CYCLE
48     
49              IF(MPPIC) THEN
50                 CALL GENERATE_PARTICLE_CONFIG_MPPIC(ICV)
51              ELSE
52                 CALL GENERATE_PARTICLE_CONFIG_DEM(ICV)
53              ENDIF
54     
55           ENDDO
56     
57           CALL GLOBAL_SUM(PIP,PARTICLES)
58     
59           WRITE(ERR_MSG, 1004) PARTICLES
60      1004 FORMAT(/,'Total number of particles in the system: ',I15)
61     
62           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
63     
64           CALL FINL_ERR_MSG
65     
66           RETURN
67           END SUBROUTINE GENERATE_PARTICLE_CONFIG
68     
69     
70     
71     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
72     !                                                                      !
73     !  SUBROUTINE: GENERATE_PARTICLE_CONFIG                                !
74     !  Authors: Rahul Garg                               Date: 21-Mar-2014 !
75     !                                                                      !
76     !  Purpose: Generate particle configuration for DEM solids for each IC !
77     !           region. Now using the particle linked lists for initial    !
78     !           build                                                      !
79     !           This routine will ultimately supersede the older rouine    !
80     !           that has not been deleted yet                              !
81     !                                                                      !
82     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
83           SUBROUTINE GENERATE_PARTICLE_CONFIG_DEM(ICV)
84     
85     
86     ! Global Variables:
87     ! particle radius and density
88           use discretelement, only: DES_RADIUS, RO_Sol
89     ! particle position new and old
90           use discretelement, only: DES_POS_NEW, DES_POS_OLD
91     ! particle velocity new and old
92           use discretelement, only: DES_VEL_NEW, DES_VEL_OLD
93     ! Simulation dimension (2D/3D)
94           use discretelement, only: DIMN
95     ! Number of particles in the system (current)
96           use discretelement, only: PIP
97     ! Number of DEM solids phases.
98           use discretelement, only: DES_MMAX
99     ! Flag to use _OLD variables
100           use discretelement, only: DO_OLD
101     ! Angular velocity
102           use discretelement, only: OMEGA_OLD, OMEGA_NEW, PIJK
103     ! DEM solid phase diameters and densities.
104           use physprop, only: D_p0, RO_s0, SMAX
105     ! IC Region solids volume fraction.
106           use ic, only: IC_EP_S
107     
108     ! Constant: 3.14159...
109           use constant, only: PI
110     ! min and max physical co-ordinates of IC regions in each direction
111           use ic, only: IC_X_w, IC_X_e, IC_Y_s, IC_Y_n, IC_Z_b, IC_Z_t
112     ! initally specified velocity field and granular temperature
113           use ic, only: IC_U_s, IC_V_s, IC_W_s, IC_Theta_M
114     ! Flag to extend the lattice distribution in a given IC to available area
115           use ic, only: IC_DES_FIT_TO_REGION
116     ! Parameter for detecting unspecified values, zero, and one
117           use param1, only: UNDEFINED, UNDEFINED_I, ZERO, ONE, Half
118     ! Parameter for small and large numbers
119           use param1, only: SMALL_NUMBER, LARGE_NUMBER
120     
121     ! to access random number generator subroutines
122           use randomno
123           use mpi_utility
124           use functions, only: SET_NORMAL, FUNIJK, FLUID_AT
125     
126           use error_manager
127     
128     ! direction wise spans of the domain and grid spacing in each direction
129           use geometry, only: xlength, ylength, zlength
130     
131           use cutcell, only : CARTESIAN_GRID, CUT_CELL_AT
132           use STL_PREPROC_DES, only: CHECK_IF_PARTICLE_OVERLAPS_STL
133           use run, only: solids_model
134           use des_allocate, only: PARTICLE_GROW
135     
136           use discretelement, only: MAX_RADIUS
137     
138           use discretelement, only: XE, YN, ZT
139     
140           use param, only: DIM_M, DIMENSION_I, DIMENSION_J, DIMENSION_K
141           use functions, only: IS_ON_MYPE_WOBND
142     
143           IMPLICIT NONE
144     
145           INTEGER, INTENT(IN) :: ICV
146     
147     ! Local variables
148     !---------------------------------------------------------------------//
149     ! Starting positions in the axial directions
150           DOUBLE PRECISION :: xINIT, yINIT, zINIT
151     ! Fractor used to scale particle diameter
152           DOUBLE PRECISION :: lFAC
153     ! Particle position and velocity
154           DOUBLE PRECISION :: POS(3), VEL(3)
155     ! Number of particles in the lattice
156           INTEGER :: SEED_X, SEED_Y, SEED_Z
157     ! Loop indices phase/fluid cell
158           INTEGER :: M, MM, I, J, K, IJK, LB, UB
159     ! Loop indicies for seeding
160           INTEGER :: II, JJ, KK
161     ! Start and end bound for IC region.
162           DOUBLE PRECISION :: IC_START(3), IC_END(3)
163     ! Volume and lengths of the IC Region
164           DOUBLE PRECISION :: DOM_VOL, DOML(3)
165     ! Flag to skip the particle
166           LOGICAL :: SKIP
167     ! Diameter adjusted for space padding
168           DOUBLE PRECISION :: ADJ_DIA
169     ! Number of particles calculated from volume fracton
170           INTEGER :: rPARTS(DIM_M), tPARTS
171     ! Spacing between particles.
172           DOUBLE PRECISION :: lDEL, lDX, lDY, lDZ
173     ! Flag that the setup failed to fit the particles to the IC region
174           LOGICAL :: FIT_FAILED
175     ! Number of seeded particles
176           INTEGER :: pCOUNT(DIM_M), tCOUNT
177     
178           DOUBLE PRECISION :: SOLIDS_DATA(0:DIM_M)
179     
180           LOGICAL :: VEL_FLUCT
181           DOUBLE PRECISION :: VEL_SIG
182           DOUBLE PRECISION, ALLOCATABLE :: randVEL(:,:)
183     
184     !......................................................................!
185     
186           CALL INIT_ERR_MSG("GENERATE_PARTICLE_CONFIG_DEM")
187     
188           WRITE(ERR_MSG,"(2/,'Generating initial particle configuration:')")
189           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
190     
191           SOLIDS_DATA = ZERO
192           CALL GET_IC_VOLUME(ICV, SOLIDS_DATA(0))
193     
194     ! setting particle seed spacing grid to be slightly greater than
195     ! the maximum particle diameter. seed at ~particle radii
196           lFAC = 1.05D0
197     
198     ! Setup local arrays with IC region bounds.
199           IC_START(1)=IC_X_W(ICV);   IC_END(1)=IC_X_E(ICV)
200           IC_START(2)=IC_Y_S(ICV);   IC_END(2)=IC_Y_N(ICV)
201           IC_START(3)=IC_Z_B(ICV);   IC_END(3)=IC_Z_T(ICV)
202     
203           DOML = IC_END-IC_START
204           IF(NO_K) DOML(3)=DZ(1)
205     
206     ! Volume of the IC region
207           DOM_VOL = DOML(1)*DOML(2)*DOML(3)
208     
209           rPARTS=0
210           VEL_FLUCT = .FALSE.
211           DO M=1,SMAX+DES_MMAX
212              IF(SOLIDS_MODEL(M) == 'DEM') THEN
213     ! Number of particles for phase M
214                 rPARTS(M) = &
215                    floor((6.0d0*IC_EP_S(ICV,M)*DOM_VOL)/(PI*(D_P0(M)**3)))
216     ! Set flags for random velocity fluctuations
217                 VEL_FLUCT = VEL_FLUCT .OR. (IC_Theta_M(ICV,M) /= 0.0d0)
218              ENDIF
219           ENDDO
220     
221     ! Total number of particles in this IC region.
222           tPARTS = sum(rPARTS)
223           IF(tPARTS == 0) RETURN
224     
225           ADJ_DIA = 2.0d0*MAX_RADIUS*lFAC
226     
227     ! Calculate velocity fluctuations if specified.
228           IF(VEL_FLUCT) THEN
229              allocate(randVEL(tPARTS,3))
230              LB=1
231              DO M=1,SMAX+DES_MMAX
232                 IF(SOLIDS_MODEL(M) == 'DEM' .AND. rPARTS(M) > 0.0d0) THEN
233                    UB=LB+int(rPARTS(M))-1
234                    IF(IC_Theta_M(ICV,1) /= 0.0d0)THEN
235                       VEL_SIG = sqrt(IC_Theta_M(ICV,1))
236                       CALL NOR_RNO(randVEL(LB:UB,1), IC_U_s(ICV,1), VEL_SIG)
237                       CALL NOR_RNO(randVEL(LB:UB,2), IC_V_s(ICV,2), VEL_SIG)
238                       IF(DO_K) CALL NOR_RNO(randVEL(LB:UB,3), &
239                          IC_W_s(ICV,3), VEL_SIG)
240                    ELSE
241                       randVEL(LB:UB,1) = ZERO
242                       randVEL(LB:UB,2) = ZERO
243                       IF(DO_K) randVEL(LB:UB,3) = ZERO
244                    ENDIF
245                    LB=UB+1
246                 ENDIF
247              ENDDO
248           ENDIF
249     
250     
251     ! Attempt to seed particle throughout the IC region
252           FIT_FAILED=.FALSE.
253           IF(IC_DES_FIT_TO_REGION(ICV)) THEN
254              IF(NO_K) THEN
255                 lDEL = (DOML(1)-ADJ_DIA)*(DOML(2)-ADJ_DIA)
256                 lDEL = (lDEL/dble(tPARTS))**(1.0/2.0)
257                 SEED_X = max(1,ceiling((DOML(1)-ADJ_DIA)/lDEL))
258                 SEED_Y = max(1,ceiling((DOML(2)-ADJ_DIA)/lDEL))
259                 SEED_Z = 1
260              ELSE
261                 lDEL = (DOML(1)-ADJ_DIA)*(DOML(2)-ADJ_DIA)*(DOML(3)-ADJ_DIA)
262                 lDEL = (lDEL/dble(tPARTS))**(1.0/3.0)
263                 SEED_X = max(1,ceiling((DOML(1)-ADJ_DIA)/lDEL))
264                 SEED_Y = max(1,ceiling((DOML(2)-ADJ_DIA)/lDEL))
265                 SEED_Z = max(1,ceiling((DOML(3)-ADJ_DIA)/lDEL))
266              ENDIF
267              FIT_FAILED=(dble(SEED_X*SEED_Y*SEED_Z) < tPARTS)
268           ENDIF
269     
270     ! Generic filling
271           IF(.NOT.IC_DES_FIT_TO_REGION(ICV) .OR. FIT_FAILED) THEN
272              SEED_X = max(1,floor((IC_END(1)-IC_START(1)-ADJ_DIA)/ADJ_DIA))
273              SEED_Y = max(1,floor((IC_END(2)-IC_START(2)-ADJ_DIA)/ADJ_DIA))
274              SEED_Z = max(1,floor((IC_END(3)-IC_START(3)-ADJ_DIA)/ADJ_DIA))
275           ENDIF
276     
277           lDX = DOML(1)/dble(SEED_X)
278           lDY = DOML(2)/dble(SEED_Y)
279           IF(DO_K) THEN
280              lDZ = DOML(3)/dble(SEED_Z)
281           ELSE
282              lDZ = 0.0d0
283           ENDIF
284     
285           xINIT = IC_START(1)+HALF*lDX
286           yINIT = IC_START(2)+HALF*lDY
287           zINIT = IC_START(3)+HALF*lDZ
288     
289     
290           M=1
291           pCOUNT = 0
292           tCOUNT = 0
293           JJ_LP: DO  JJ=1, SEED_Y
294           KK_LP: DO  KK=1, SEED_Z
295           II_LP: DO  II=1, SEED_X
296     
297     ! Exit if all particles were seeded.
298              IF(tCOUNT > int(tPARTS)) THEN
299                 EXIT JJ_LP
300     ! Find the next phase that needs to be seeded
301              ELSEIF(pCOUNT(M) > int(rPARTS(M))) THEN
302                 MM_LP: DO MM=M+1,SMAX+DES_MMAX
303                    IF(rPARTS(MM) > 0.0) THEN
304                       M=MM
305                       EXIT MM_LP
306                    ENDIF
307                 ENDDO MM_LP
308                 IF(MM > SMAX+DES_MMAX) THEN
309                    EXIT JJ_LP
310                 ELSEIF(IC_Theta_M(ICV,MM) /= 0.0d0) THEN
311                    VEL_SIG = sqrt(IC_Theta_M(ICV,MM))
312                    CALL NOR_RNO(randVEL(:,1), 0.0d0, VEL_SIG)
313                    CALL NOR_RNO(randVEL(:,2), 0.0d0, VEL_SIG)
314                    IF(DO_K) CALL NOR_RNO(randVEL(:,3), 0.0d0, VEL_SIG)
315                 ENDIF
316              ENDIF
317     
318              pCOUNT(M) = pCOUNT(M) + 1
319              tCOUNT = tCOUNT + 1
320     
321     ! Position the particle
322              POS(1) = xINIT + (II-1)*lDX
323              POS(2) = YINIT + (JJ-1)*lDY
324              POS(3) = ZINIT + (KK-1)*lDZ
325     
326     ! Bin the parcel to the fuild grid.
327              K=1
328              IF(DO_K) CALL PIC_SEARCH(K, POS(3), ZT, DIMENSION_K, KMIN2, KMAX2)
329              CALL PIC_SEARCH(J, POS(2), YN, DIMENSION_J, JMIN2, JMAX2)
330              CALL PIC_SEARCH(I, POS(1), XE, DIMENSION_I, IMIN2, IMAX2)
331     
332     ! Skip cells that are not part of the local fuild domain.
333              IF(.NOT.IS_ON_MYPE_WOBND(I,J,K)) CYCLE
334              IF(DEAD_CELL_AT(I,J,K)) CYCLE
335     
336              IJK = FUNIJK(I,J,K)
337              IF(.NOT.FLUID_AT(IJK)) CYCLE
338     
339              IF(CARTESIAN_GRID) THEN
340                 CALL CHECK_IF_PARTICLE_OVERLAPS_STL(POS, I, J, K, SKIP)
341                 IF(SKIP) CYCLE
342              ENDIF
343     
344              PIP = PIP + 1
345              CALL PARTICLE_GROW(PIP)
346     
347              CALL SET_NORMAL(PIP)
348     
349              IF(VEL_FLUCT) THEN
350                 VEL(1) = randVEL(tCOUNT,1)
351                 VEL(2) = randVEL(tCOUNT,2)
352                 VEL(3) = randVEL(tCOUNT,3)
353              ELSE
354                 VEL(1) = IC_U_s(ICV,M)
355                 VEL(2) = IC_V_s(ICV,M)
356                 VEL(3) = IC_W_s(ICV,M)
357              ENDIF
358              IF(NO_K) VEL(3) = 0.0d0
359     
360     
361              DES_POS_NEW(:,PIP) = POS(:)
362              DES_VEL_NEW(:,PIP) = VEL(:)
363              OMEGA_NEW(:,PIP) = 0.0d0
364     
365              DES_RADIUS(PIP) = D_P0(M)*HALF
366              RO_SOL(PIP) =  RO_S0(M)
367     
368              PIJK(PIP,1) = I
369              PIJK(PIP,2) = J
370              PIJK(PIP,3) = K
371              PIJK(PIP,4) = IJK
372              PIJK(PIP,5) = M
373     
374              IF(DO_OLD) THEN
375                 DES_VEL_OLD(:,PIP) = DES_VEL_NEW(:,PIP)
376                 DES_POS_OLD(:,PIP) = DES_POS_NEW(:,PIP)
377                 OMEGA_OLD(:,PIP) = ZERO
378              ENDIF
379     
380              SOLIDS_DATA(M) = SOLIDS_DATA(M) + 1.0
381     
382           ENDDO II_LP
383           ENDDO KK_LP
384           ENDDO JJ_LP
385     
386     ! Collect the data
387           CALL GLOBAL_ALL_SUM(SOLIDS_DATA)
388     
389     ! Verify that the IC region volume is not zero.
390           IF(SOLIDS_DATA(0) <= 0.0d0) THEN
391              WRITE(ERR_MSG,1000) ICV, SOLIDS_DATA(0)
392              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
393           ENDIF
394      
395     1000 FORMAT('Error 1000: Invalid IC region volume: IC=',I3,' VOL=',&
396              ES15.4,/'Please correct the mfix.dat file.')
397     
398           WRITE(ERR_MSG,2000) ICV
399           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
400     
401           DO M=1, SMAX+DES_MMAX
402              IF(SOLIDS_DATA(M) < SMALL_NUMBER) CYCLE
403              WRITE(ERR_MSG,2010) M, int(SOLIDS_DATA(M)), IC_EP_S(ICV,M),   &
404                 (dble(SOLIDS_DATA(M))*(Pi/6.0d0)*D_P0(M)**3)/SOLIDS_DATA(0)
405              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
406           ENDDO
407     
408           IF(allocated(randVEL)) deallocate(randVEL)
409     
410           CALL FINL_ERR_MSG
411     
412           RETURN
413     
414      2000 FORMAT(/2x,'|',43('-'),'|',/2x,'| IC Region: ',I3,28x,'|',/2x,   &
415              '|',43('-'),'|',/2x,'| Phase | Number of |    EPs    |    EP',&
416              's    |',/2x,'|   ID  | Particles | Specified |   Actual  |', &
417              /2x,'|-------|',3(11('-'),'|'))
418     
419      2010 FORMAT(2x,'|  ',I3,'  |',1x,I9,1x,'|',2(1x,ES9.2,1x,'|'),/2x,    &
420              '|-------|',3(11('-'),'|'))
421     
422           END SUBROUTINE GENERATE_PARTICLE_CONFIG_DEM
423     
424     
425     
426     
427     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
428     !                                                                      !
429     !  Subroutine: GENERATE_PARTICLE_CONFIG_MMPPIC                         !
430     !  Author: Rahul Garg                                 Date: 3-May-2011 !
431     !                                                                      !
432     !  Purpose: Generates particle position distribution for MPPIC.        !
433     !                                                                      !
434     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
435           SUBROUTINE GENERATE_PARTICLE_CONFIG_MPPIC(ICV)
436     
437     ! Number of DES solids phases.
438           use discretelement, only: DES_MMAX
439     ! Flag indicating that the IC region is defined.
440           use ic, only: IC_DEFINED
441     ! IC Region bulk density (RO_s * EP_s)
442           use ic, only: IC_ROP_s
443     ! IC Region gas volume fraction.
444           use ic, only: IC_EP_G
445     ! MPPIC specific IC region specification.
446           use ic, only: IC_PIC_CONST_NPC, IC_PIC_CONST_STATWT
447     
448           use param1, only: UNDEFINED, UNDEFINED_I
449           use param1, only: ZERO, ONE, HALF
450     ! Maximum number of IC regions and solids phases
451           use param, only: DIMENSION_IC
452           use param, only: DIM_M
453     
454           use mpi_utility, only: GLOBAL_ALL_SUM
455     
456           use error_manager
457     
458     
459           IMPLICIT NONE
460     
461     
462           INTEGER, INTENT(IN) :: ICV
463     
464     ! Local variables
465     !---------------------------------------------------------------------//
466     ! Generic loop counters
467           INTEGER :: M
468     ! Actual volume of IC region
469           DOUBLE PRECISION :: IC_VOL
470     ! Solids data in IC Region by phase:
471           DOUBLE PRECISION :: SOLIDS_DATA(0:4*DIM_M)
472     !......................................................................!
473     
474           CALL INIT_ERR_MSG("GENERATE_PARTICLE_CONFIG_MPPIC")
475     
476           WRITE(ERR_MSG,"(2/,'Generating initial parcel configuration:')")
477           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
478     
479     
480           SOLIDS_DATA = ZERO
481           CALL GET_IC_VOLUME(ICV, SOLIDS_DATA(0))
482     
483           WRITE(ERR_MSG,2000) ICV
484           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
485     
486     ! Set up the individual solids phases.
487           DO M=1, DES_MMAX
488              IF(IC_ROP_S(ICV,M) == ZERO) CYCLE
489     ! Seed parcels with constant stastical weight.
490              IF(IC_PIC_CONST_STATWT(ICV,M) /= ZERO) THEN
491                 CALL GPC_MPPIC_CONST_STATWT(ICV, M, SOLIDS_DATA(0), &
492                    SOLIDS_DATA((4*M-3):(4*M)))
493     ! Seed parcels with a constant number per cell
494              ELSEIF(IC_PIC_CONST_NPC(ICV,M) /= 0) THEN
495                 CALL GPC_MPPIC_CONST_NPC(ICV, M, SOLIDS_DATA(0), &
496                    SOLIDS_DATA((4*M-3):(4*M)))
497              ENDIF
498           ENDDO
499     
500     ! Collect the data
501           CALL GLOBAL_ALL_SUM(SOLIDS_DATA)
502     
503     ! Verify that the IC region volume is not zero.
504           IF(SOLIDS_DATA(0) <= 0.0d0) THEN
505              WRITE(ERR_MSG,1000) ICV, SOLIDS_DATA(0)
506              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
507           ENDIF
508     
509      1000 FORMAT('Error 1000: Invalid IC region volume: IC=',I3,' VOL=',&
510              ES15.4,/'Please correct the mfix.dat file.')
511     
512     ! Report solids information for the IC region.
513           DO M=1, DES_MMAX
514              WRITE(ERR_MSG,2010) M, int(SOLIDS_DATA(4*M-3)),&
515                 int(SOLIDS_DATA(4*M-3)*SOLIDS_DATA(4*M-2)), &
516                 SOLIDS_DATA(4*M-2), SOLIDS_DATA(4*M-1),     &
517                 SOLIDS_DATA(4*M)/SOLIDS_DATA(0)
518              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
519           ENDDO
520     
521           CALL FINL_ERR_MSG
522     
523           RETURN
524     
525      2000 FORMAT(/2x,'|',67('-'),'|',/2x,'| IC Region: ',I3,52x,'|',/2x,   &
526              '|',67('-'),'|',/2x,'| Phase | Num. Comp | Num. Real ',       &
527              '| Stastical |    EPs    |    EPs    |',/2x,'|   ID  |  ',    &
528              'Parcels  | Particles |   Weight  | Specified |   Actual  |', &
529              /2x,'|-------|',5(11('-'),'|'))
530     
531      2010 FORMAT(2x,'|  ',I3,'  |',2(1x,I9,1x,'|'),3(1x,ES9.2,1x,'|'),/2x,&
532              '|-------|',5(11('-'),'|'))
533     
534           END SUBROUTINE GENERATE_PARTICLE_CONFIG_MPPIC
535     
536     
537     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
538     !                                                                      !
539     !  Subroutine: GPC_MPPIC_CONST_STATWT                                  !
540     !  Author: J.Musser                                 Date:  26-Aug-2015 !
541     !                                                                      !
542     !  Purpose: generates particle position distribution for MPPIC         !
543     !                                                                      !
544     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
545           SUBROUTINE GPC_MPPIC_CONST_STATWT(ICV, M, IC_VOL, sDATA)
546     
547     !-----------------------------------------------
548     ! Modules
549     !-----------------------------------------------
550           use cutcell, only : CARTESIAN_GRID
551     
552     ! IC Region Bounds
553           use ic, only: IC_X_w, IC_X_e
554           use ic, only: IC_Y_s, IC_Y_n
555           use ic, only: IC_Z_b, IC_Z_t
556     ! IC Region bulk density (RO_s * EP_s)
557           use ic, only: IC_ROP_s
558     ! IC Region velocity field and granular temperature
559           use ic, only: IC_U_s, IC_V_s, IC_W_s, IC_Theta_M
560     ! MPPIC specific IC region specification.
561           use ic, only: IC_PIC_CONST_STATWT
562     ! Flag for Cartesian cut-cell 
563           use cutcell, only : CARTESIAN_GRID
564     ! DES solid phase diameters and densities.
565           use discretelement, only: DES_D_p0, DES_RO_s
566     ! Total number of active parcles
567           use discretelement, only: PIP
568     ! Parcel position and velocity
569           use discretelement, only: DES_POS_NEW, DES_VEL_NEW
570     ! Parcel radius and density
571           use discretelement, only: DES_RADIUS, RO_SOL
572     ! Stastical weight of each parcel
573           use mfix_pic, only: DES_STAT_WT
574     ! Paricle fluid I/J/K/IJK/Phase information
575           use discretelement, only: PIJK
576     ! The east/north/top location of grid partition
577           use discretelement, only: XE, YN, ZT
578     ! Fluid grid cell dimensions and mesh size
579           USE geometry, only: IMIN2, IMAX2
580           USE geometry, only: JMIN2, JMAX2
581           USE geometry, only: KMIN2, KMAX2
582           USE geometry, only: DO_K, NO_K, DZ
583     
584     ! Global parameters
585     !----------------------------------------------------------------------//
586           use param1, only: ZERO, HALF, ONE
587     ! Maximum cell partitions in the axial directions
588           use param, only: DIMENSION_I, DIMENSION_J, DIMENSION_K
589     ! Maximum number of IC regions and solids phases
590           use param, only: DIMENSION_IC
591     ! Constant: 3.14159...
592           use constant, only: PI
593     
594     
595     ! Module procedures
596     !----------------------------------------------------------------------//
597           use functions, only: FUNIJK, FLUID_AT, DEAD_CELL_AT
598           use stl_preproc_des, only: CHECK_IF_PARTICLE_OVERLAPS_STL
599           use cutcell, only: cut_cell_at
600           use randomno, only: UNI_RNO, NOR_RNO
601           use functions, only: SET_NORMAL
602           use functions, only: IS_ON_MYPE_WOBND
603           use des_allocate, only: PARTICLE_GROW
604     
605           use error_manager
606     
607     
608           IMPLICIT NONE
609     
610     
611     ! Dummy Arguments
612     !----------------------------------------------------------------------//
613     ! Index of IC region and solids phase
614           INTEGER, INTENT(IN) :: ICV, M
615     ! Specific volume of IC region (accounts for blocked cells)
616           DOUBLE PRECISION, INTENT(IN) :: IC_VOL
617     ! Data about solids in the IC region.
618           DOUBLE PRECISION, INTENT(OUT) :: sDATA(4)
619     
620     ! Local variables
621     !----------------------------------------------------------------------//
622     ! Solids phase M volume fraction in IC region
623           DOUBLE PRECISION :: EP_SM
624     ! Start/End bounds of IC region
625           DOUBLE PRECISION :: IC_START(3), IC_END(3)
626     ! IC region lenghts and volume
627           DOUBLE PRECISION :: DOML(3), DOM_VOL
628     ! Number of real particles calculated from volume fracton
629           DOUBLE PRECISION ::  rPARTS
630     ! Number of computational parcels
631           INTEGER :: cPARTS
632     ! Calculated statistical weight of parcels
633           DOUBLE PRECISION :: STAT_WT
634     ! Volume of a parcel and total solids volume
635           DOUBLE PRECISION :: sVOL, sVOL_TOT
636     ! Parcel position and velocity
637           DOUBLE PRECISION :: POS(3), VEL(3)
638     ! Arrays for assigning random position and velocities
639           DOUBLE PRECISION, ALLOCATABLE :: randPOS(:,:), randVEL(:,:)
640     ! Standard deivation of initial velocities
641           DOUBLE PRECISION :: VEL_SIG
642     ! Flag to skip parcel
643           LOGICAL :: SKIP
644     ! Counter for seeded parcles.
645           INTEGER :: SEEDED
646     ! Generic fluid cell indices and loop counters
647           INTEGER :: I, J, K, IJK, LC
648     !......................................................................!
649     
650           CALL INIT_ERR_MSG("GENERATE_PARTICLE_CONFIG_MPPIC")
651     
652     ! Setup local arrays with IC region bounds.
653           IC_START(1)=IC_X_W(ICV);   IC_END(1)=IC_X_E(ICV)
654           IC_START(2)=IC_Y_S(ICV);   IC_END(2)=IC_Y_N(ICV)
655           IC_START(3)=IC_Z_B(ICV);   IC_END(3)=IC_Z_T(ICV)
656     
657           DOML = IC_END-IC_START
658           IF(NO_K) DOML(3)=DZ(1)
659     
660     ! Volume of the IC region
661           DOM_VOL = DOML(1)*DOML(2)*DOML(3)
662     ! Solids volume fraction in IC region
663           EP_SM = IC_ROP_S(ICV,M)/DES_RO_s(M)
664     
665     ! Total number of real and computational particles
666           rPARTS = 6.d0*EP_SM*DOM_VOL/(PI*(Des_D_P0(M)**3.d0))
667           cPARTS = max(1,int(rPARTS/real(IC_PIC_CONST_STATWT(ICV,M))))
668     
669     ! The 'actual' statistical weight
670           STAT_WT = rPARTS/real(cPARTS)
671     
672     ! Obtain random variations for parcel position 
673           ALLOCATE(randPOS(cPARTS,3))
674     
675           DO LC = 1, merge(2,3,NO_K)
676              CALL UNI_RNO(RANDPOS(1:cPARTS,LC))
677           ENDDO
678           IF(NO_K) randPOS(:,3) = 0.0d0
679     
680     ! Obtain initial velocities with random variations
681           ALLOCATE(randVEL(cPARTS,3))
682     
683           randVEL(:,1) = IC_U_s(ICV,M)
684           randVEL(:,2) = IC_V_s(ICV,M)
685           randVEL(:,3) = merge(IC_W_s(ICV,M), 0.0d0, DO_K)
686     
687           VEL_SIG = sqrt(IC_Theta_M(ICV,M))
688           IF(VEL_SIG /= ZERO) THEN
689              CALL NOR_RNO(randVEL(:,1), IC_U_S(ICV,M), VEL_SIG)
690              CALL NOR_RNO(randVEL(:,2), IC_V_S(ICV,M), VEL_SIG)
691              IF(DO_K) CALL NOR_RNO(randVEL(:,3), IC_W_S(ICV,M), VEL_SIG)
692           ENDIF
693     
694     ! Volume occupied by one parcel
695           sVOL = (Pi/6.0d0)*(DES_D_P0(M)**3.d0)*STAT_WT
696           sVOL_TOT = IC_VOL*EP_SM
697     
698           SEEDED = 0
699           DO LC = 1, cPARTS
700     
701     ! Generate the initial parcel position.
702              POS = IC_START + DOML*RANDPOS(LC,:)
703              VEL = 0.0d0
704     
705     ! Bin the parcel to the fuild grid.
706              K=1
707              IF(DO_K) CALL PIC_SEARCH(K, POS(3), ZT, DIMENSION_K, KMIN2, KMAX2)
708              CALL PIC_SEARCH(J, POS(2), YN, DIMENSION_J, JMIN2, JMAX2)
709              CALL PIC_SEARCH(I, POS(1), XE, DIMENSION_I, IMIN2, IMAX2)
710     
711     ! Skip cells that are not part of the local fuild domain.
712              IF(.NOT.IS_ON_MYPE_WOBND(I,J,K)) CYCLE
713              IF(DEAD_CELL_AT(I,J,K)) CYCLE
714     
715              IJK = FUNIJK(I, J, K)
716              IF(.NOT.FLUID_AT(IJK)) CYCLE
717     
718              IF(CARTESIAN_GRID) THEN
719                 CALL CHECK_IF_PARCEL_OVERLAPS_STL(POS, SKIP)
720                 IF(SKIP) CYCLE
721              ENDIF
722     
723              PIP = PIP + 1
724              CALL PARTICLE_GROW(PIP)
725     
726              DES_POS_NEW(:,PIP) = POS(:)
727              DES_VEL_NEW(:,PIP) = VEL(:)
728     
729              DES_RADIUS(PIP) = DES_D_P0(M)*HALF
730              RO_SOL(PIP) =  DES_RO_S(M)
731     
732              PIJK(PIP,1) = I
733              PIJK(PIP,2) = J
734              PIJK(PIP,3) = K
735              PIJK(PIP,4) = IJK
736              PIJK(PIP,5) = M
737     
738              DES_STAT_WT(PIP) = STAT_WT
739     
740              CALL SET_NORMAL(PIP)
741     
742              SEEDED = SEEDED + 1
743     
744              IF(sVOL_TOT <= sVOL*dble(SEEDED)) EXIT
745           ENDDO
746     
747           sDATA(1) = dble(SEEDED)
748           sDATA(2) = STAT_WT
749           sDATA(3) = EP_SM
750           sDATA(4) = sVOL*dble(SEEDED)
751     
752           IF(allocated(randPOS)) deallocate(randPOS)
753           IF(allocated(randPOS)) deallocate(randVEL)
754     
755           CALL FINL_ERR_MSG
756     
757           RETURN
758           END SUBROUTINE GPC_MPPIC_CONST_STATWT
759     
760     
761     
762     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
763     !                                                                      !
764     !  Subroutine: GENERATE_PARTICLE_CONFIG_MMPPIC                         !
765     !  Author: Rahul Garg                                 Date: 3-May-2011 !
766     !                                                                      !
767     !  Purpose: generates particle position distribution for MPPIC         !
768     !                                                                      !
769     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
770           SUBROUTINE GPC_MPPIC_CONST_NPC(ICV, M, IC_VOL, sDATA)
771     
772     ! Global variables
773     !---------------------------------------------------------------------//
774           use cutcell, only : CARTESIAN_GRID
775           use discretelement, only: XE, YN, ZT
776     
777     ! Number of DES solids phases.
778           use discretelement, only: DES_MMAX
779     
780     ! DEM solid phase diameters and densities.
781           use discretelement, only: DES_D_p0, DES_RO_s
782     
783     ! Implied total number of physical particles
784           use mfix_pic, only: rnp_pic
785     ! Total number of computational particles/parcels
786           use mfix_pic, only: cnp_pic
787     
788           use mfix_pic, only: des_stat_wt
789     
790     ! Flag indicating that the IC region is defined.
791           use ic, only: IC_DEFINED
792     ! IC Region bulk density (RO_s * EP_s)
793           use ic, only: IC_EP_s
794     
795           use ic, only: IC_I_w, IC_I_e, IC_J_s, IC_J_n, IC_K_b, IC_K_t
796     
797     ! initally specified velocity field and granular temperature
798           use ic, only: IC_U_s, IC_V_s, IC_W_s, IC_Theta_M
799           use ic, only: IC_PIC_CONST_NPC
800     
801     ! Cut_cell identifier array
802           use cutcell, only: cut_cell_at
803     
804     ! Maximum number of IC regions and solids phases
805           use param, only: DIMENSION_IC
806           use param, only: DIM_M
807           use param1, only: UNDEFINED, UNDEFINED_I, ZERO, ONE, HALF
808     
809     ! Constant: 3.14159...
810           use constant, only: PI
811     
812           use mpi_utility
813     
814           use randomno
815           use error_manager
816           use functions
817     !      use STL_PREPROC_DES, only: CHECK_IF_PARTICLE_OVERLAPS_STL
818           use run, only: solids_model
819           use des_allocate, only: PARTICLE_GROW
820     
821           IMPLICIT NONE
822     
823     ! Dummy Arguments
824     !----------------------------------------------------------------------//
825     ! Index of IC region and solids phase
826           INTEGER, INTENT(IN) :: ICV, M
827     ! Specific volume of IC region (accounts for blocked cells)
828           DOUBLE PRECISION, INTENT(IN) :: IC_VOL
829     ! Data about solids in the IC region.
830           DOUBLE PRECISION, INTENT(OUT) :: sDATA(4)
831     
832     ! Local variables
833     !----------------------------------------------------------------------//
834           DOUBLE PRECISION :: EP_SM
835     
836     
837     ! Number of real and comp. particles in a cell.
838           DOUBLE PRECISION ::  rPARTS
839           INTEGER :: cPARTS
840           DOUBLE PRECISION :: DOML(3), IC_START(3)
841     ! Parcel position with random
842           DOUBLE PRECISION :: POS(3)
843     ! Average velocity and standard deivation
844           DOUBLE PRECISION :: IC_VEL(3), VEL_SIG
845     ! Flag to not keep parcel.
846           LOGICAL :: SKIP
847     ! Arrasy for assigning random position and velocities
848           DOUBLE PRECISION, ALLOCATABLE :: randPOS(:,:), randVEL(:,:)
849     ! Statistical weights
850           DOUBLE PRECISION :: STAT_WT, SUM_STAT_WT
851     ! Volume of a parcel and total solids volume
852           DOUBLE PRECISION :: sVOL, sVOL_TOT
853     ! Counter for seeded parcles.
854           INTEGER :: SEEDED
855     ! Generic loop indices and loop counters
856           INTEGER :: I, J, K, IJK, LC, LC_MAX
857     !......................................................................!
858     
859           CALL INIT_ERR_MSG("GPC_MPPIC_CONST_NPC")
860     
861           cPARTS = IC_PIC_CONST_NPC(ICV,M)
862     
863           allocate(randPOS(cPARTS,3))
864           allocate(randVEL(cPARTS,3))
865     
866           IC_VEL(1) = IC_U_s(ICV,M)
867           IC_VEL(2) = IC_V_s(ICV,M)
868           IC_VEL(3) = merge(IC_W_s(ICV,M),0.0d0,DO_K)
869     
870           VEL_SIG = sqrt(IC_Theta_M(ICV,M))
871     
872     ! Volume occupied by one particle
873           sVOL = (Pi/6.0d0)*(DES_D_P0(M)**3.d0)
874     
875           SEEDED = 0
876           sVOL_TOT = 0.0d0
877           SUM_STAT_WT = 0.0d0
878     
879           DO K = IC_K_B(ICV), IC_K_T(ICV)
880           DO J = IC_J_S(ICV), IC_J_N(ICV)
881           DO I = IC_I_W(ICV), IC_I_E(ICV)
882     
883              IF(.not.IS_ON_myPE_wobnd(I,J,K)) cycle
884              IF(DEAD_CELL_AT(I,J,K)) cycle
885     
886              IJK = FUNIJK(I,J,K)
887              IF(.not.FLUID_AT(IJK)) cycle
888     
889              rPARTS = IC_EP_s(ICV,M)*VOL(IJK)/sVOL
890     
891              LC_MAX = cPARTS
892              IF(CUT_CELL_AT(IJK)) LC_MAX = &
893                 max(1, int(VOL(IJK)*dble(cPARTS)/(DX(I)*DY(J)*DZ(K))))
894     
895              DO LC=1, merge(2,3,NO_K)
896                 CALL UNI_RNO(RANDPOS(1:LC_MAX,LC))
897                 IF(VEL_SIG > ZERO) THEN
898                    CALL NOR_RNO(randVEL(1:LC_MAX,LC), IC_VEL(LC), VEL_SIG)
899                 ELSE
900                    randVEL(1:LC_MAX,LC) = IC_VEL(LC)
901                 ENDIF
902              ENDDO
903              IF(NO_K) THEN
904                 randPOS(1:LC_MAX,3) = 0.0d0
905                 randVEL(1:LC_MAX,3) = 0.0d0
906              ENDIF
907     
908              IC_START(1) = XE(I-1)
909              IC_START(2) = YN(J-1)
910              IC_START(3) = ZERO;  IF(DO_K) IC_START(3) = ZT(K-1)
911     
912              DOML(1) = DX(I)
913              DOML(2) = DY(J)
914              DOML(3) = ZERO;  IF(DO_K) DOML(3) = DZ(K)
915     
916              DO LC=1,LC_MAX
917     
918                 POS(:) = IC_START + DOML*randPOS(LC,:)
919     
920                 IF(CARTESIAN_GRID) THEN
921                    CALL CHECK_IF_PARCEL_OVERLAPS_STL(POS, SKIP)
922                    IF(SKIP) CYCLE
923                 ENDIF
924     
925                 PIP = PIP + 1
926                 CALL PARTICLE_GROW(PIP)
927      
928                 DES_POS_NEW(:,PIP) = POS(:)
929                 DES_VEL_NEW(:,PIP) = randVEL(LC,:)
930      
931                 DES_RADIUS(PIP) = DES_D_P0(M)*HALF
932                 RO_SOL(PIP) =  DES_RO_S(M)
933      
934                 PIJK(PIP,1) = I
935                 PIJK(PIP,2) = J
936                 PIJK(PIP,3) = K
937                 PIJK(PIP,4) = IJK
938                 PIJK(PIP,5) = M
939      
940                 STAT_WT = rPARTS/dble(LC_MAX)
941                 SUM_STAT_WT = SUM_STAT_WT + STAT_WT
942     
943                 DES_STAT_WT(PIP) = STAT_WT
944                 sVOL_TOT = sVOL_TOT + sVOL*STAT_WT
945      
946                 CALL SET_NORMAL(PIP)
947      
948                 SEEDED = SEEDED + 1
949     
950              ENDDO
951     
952           ENDDO
953           ENDDO
954           ENDDO
955     
956           IF(allocated(randPOS)) deallocate(randPOS)
957           IF(allocated(randPOS)) deallocate(randVEL)
958     
959           sDATA(1) = dble(SEEDED)
960           sDATA(2) = SUM_STAT_WT/dble(SEEDED)
961           sDATA(3) = IC_EP_S(ICV,M)
962           sDATA(4) = sVOL_TOT
963     
964           CALL FINL_ERR_MSG
965     
966           RETURN
967           END SUBROUTINE GPC_MPPIC_CONST_NPC
968     
969     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
970     !                                                                      !
971     !  Subroutine: GET_IC_VOLUME                                           !
972     !  Author: J.Musser                                 Date: 26-Aug-2015  !
973     !                                                                      !
974     !  Purpose: Calculate the actual volume of the IC region.              !
975     !                                                                      !
976     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
977           SUBROUTINE GET_IC_VOLUME(ICV, IC_VOL)
978     
979     ! IC region index bounds
980           use ic, only: IC_I_w, IC_I_e
981           use ic, only: IC_J_s, IC_J_n
982           use ic, only: IC_K_b, IC_K_t
983     
984     ! Volume of computational cells.
985           use geometry, only: VOL
986     
987           use functions
988     
989           IMPLICIT NONE
990     
991     ! Dummy Arguments
992     !---------------------------------------------------------------------//
993     ! Index of IC region
994           INTEGER, INTENT(IN) :: ICV
995     ! Total calculated volume of IC region
996           DOUBLE PRECISION, INTENT(OUT) :: IC_VOL
997     
998     ! Local variables
999     !---------------------------------------------------------------------//
1000           INTEGER :: I, J, K, IJK
1001     !......................................................................!
1002     
1003     
1004           IC_VOL = 0.0d0
1005           DO K = IC_K_B(ICV), IC_K_T(ICV)
1006           DO J = IC_J_S(ICV), IC_J_N(ICV)
1007           DO I = IC_I_W(ICV), IC_I_E(ICV)
1008     
1009              IF(.NOT.IS_ON_MYPE_WOBND(I,J,K)) CYCLE
1010              IF(DEAD_CELL_AT(I,J,K)) CYCLE
1011     
1012              IJK = FUNIJK(I,J,K)
1013              IF(FLUID_AT(IJK)) IC_VOL = IC_VOL + VOL(IJK)
1014     
1015           ENDDO
1016           ENDDO
1017           ENDDO
1018     
1019           RETURN
1020           END SUBROUTINE GET_IC_VOLUME
1021     
1022           END MODULE GENERATE_PARTICLES
1023