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

1           MODULE GENERATE_PARTICLES
2     
3             DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: PARTICLE_COUNT
4     
5           CONTAINS
6     
7     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
8     !                                                                      C
9     !  SUBROUTINE: GENERATE_PARTICLE_CONFIG                                C
10     !                                                                      C
11     !  Purpose: Generate particle configuration based on maximum particle  C
12     !           radius and filling from top to bottom within specified     C
13     !           bounds                                                     C
14     !                                                                      C
15     !                                                                      C
16     !  Authors: Rahul Garg                                Date: 19-Mar-14  C
17     !                                                                      C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20           SUBROUTINE GENERATE_PARTICLE_CONFIG
21     
22           use mfix_pic, only: MPPIC
23           use discretelement, only: PIP, PARTICLES
24     ! Flag indicating that the IC region is defined.
25           use ic, only: IC_DEFINED
26     ! Parameter for detecting unspecified values, zero, and one
27           use param1, only: ONE
28     ! Maximum number of initial conditions
29           use param, only: DIMENSION_IC
30     ! IC Region gas volume fraction.
31           use ic, only: IC_EP_G
32     
33           use mpi_utility
34     
35     ! Use the error manager for posting error messages.
36     !---------------------------------------------------------------------//
37           use error_manager
38     
39           IMPLICIT NONE
40     
41           INTEGER :: ICV
42     
43     ! Initialize the error manager.
44           CALL INIT_ERR_MSG("Generate_Particle_Config")
45     
46           DO ICV = 1, DIMENSION_IC
47     
48              IF(.NOT.IC_DEFINED(ICV)) CYCLE
49              IF(IC_EP_G(ICV) == ONE) CYCLE
50     
51              IF(MPPIC) THEN
52                 CALL GENERATE_PARTICLE_CONFIG_MPPIC(ICV)
53              ELSE
54                 CALL GENERATE_PARTICLE_CONFIG_DEM(ICV)
55              ENDIF
56     
57           ENDDO
58     
59           CALL GLOBAL_SUM(PIP,PARTICLES)
60     
61           WRITE(ERR_MSG, 1004) PARTICLES
62      1004 FORMAT(/,'Total number of particles in the system: ',I15)
63     
64           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
65     
66           CALL FINL_ERR_MSG
67     
68           RETURN
69           END SUBROUTINE GENERATE_PARTICLE_CONFIG
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     ! Global Variables:
86     !---------------------------------------------------------------------//
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     ! solid phase diameters and densities.
104           use physprop, only: D_p0, RO_s0, MMAX
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
125     
126           use desgrid, only: dg_xstart, dg_ystart, dg_zstart
127           use desgrid, only: dg_xend, dg_yend, dg_zend
128     
129     ! direction wise spans of the domain and grid spacing in each direction
130           use geometry, only: xlength, ylength, zlength
131     
132           use cutcell, only : CARTESIAN_GRID
133           use stl_functions_des, only: CHECK_IF_PARTICLE_OVERLAPS_STL
134           use run, only: solids_model
135           use des_allocate, only: PARTICLE_GROW
136     
137           use desgrid, only: IofPOS, JofPOS, KofPOS
138           use desgrid, only: dg_is_ON_myPE_OWNs
139           use toleranc, only: compare
140     
141           use discretelement, only: max_pip, max_radius, xe, yn, zt
142           use error_manager
143           use functions
144           use param, only: dim_m
145           use param, only: dimension_i, dimension_j, dimension_k
146     
147           IMPLICIT NONE
148     
149     ! Dummy arguments
150     !---------------------------------------------------------------------//
151           INTEGER, INTENT(IN) :: ICV
152     
153     ! Local variables
154     !---------------------------------------------------------------------//
155     ! Starting positions in the axial directions
156           DOUBLE PRECISION :: xINIT, yINIT, zINIT
157     ! Fractor used to scale particle diameter
158           DOUBLE PRECISION :: lFAC
159     ! Particle position and velocity
160           DOUBLE PRECISION :: POS(3), VEL(3)
161     ! Number of particles in the lattice
162           INTEGER :: SEED_X, SEED_Y, SEED_Z
163     ! Loop indices phase/fluid cell
164           INTEGER :: M, MM, I, J, K, IJK, LB, UB
165     ! Loop indicies for seeding
166           INTEGER :: II, JJ, KK
167     ! Start and end bound for IC region.
168           DOUBLE PRECISION :: IC_START(3), IC_END(3)
169     ! Volume and lengths of the IC Region
170           DOUBLE PRECISION :: DOM_VOL, DOML(3)
171     ! Flag to skip the particle
172           LOGICAL :: SKIP
173     ! Diameter adjusted for space padding
174           DOUBLE PRECISION :: ADJ_DIA
175     ! Number of particles calculated from volume fracton
176           INTEGER :: rPARTS(DIM_M), tPARTS
177     ! Spacing between particles.
178           DOUBLE PRECISION :: lDEL, lDX, lDY, lDZ
179     ! Flag that the setup failed to fit the particles to the IC region
180           LOGICAL :: FIT_FAILED
181     ! Number of seeded particles
182           INTEGER :: pCOUNT(DIM_M), tCOUNT
183     
184           DOUBLE PRECISION :: SOLIDS_DATA(0:DIM_M)
185     
186           LOGICAL :: VEL_FLUCT
187           DOUBLE PRECISION :: VEL_SIG
188           DOUBLE PRECISION, ALLOCATABLE :: randVEL(:,:)
189     
190           logical :: report = .true.
191           logical :: found
192     
193     !......................................................................!
194     
195           CALL INIT_ERR_MSG("GENERATE_PARTICLE_CONFIG_DEM")
196     
197           WRITE(ERR_MSG,"(2/,'Generating initial particle configuration:')")
198           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
199     
200           SOLIDS_DATA = ZERO
201           CALL GET_IC_VOLUME(ICV, SOLIDS_DATA(0))
202     
203     ! setting particle seed spacing grid to be slightly greater than
204     ! the maximum particle diameter. seed at ~particle radii
205           lFAC = 1.05D0
206     
207     ! Setup local arrays with IC region bounds.
208           IC_START(1)=IC_X_W(ICV);   IC_END(1)=IC_X_E(ICV)
209           IC_START(2)=IC_Y_S(ICV);   IC_END(2)=IC_Y_N(ICV)
210           IC_START(3)=IC_Z_B(ICV);   IC_END(3)=IC_Z_T(ICV)
211     
212           DOML = IC_END-IC_START
213           IF(NO_K) DOML(3)=DZ(1)
214     
215     ! Volume of the IC region
216           DOM_VOL = DOML(1)*DOML(2)*DOML(3)
217     
218           rPARTS=0
219           DO M=MMAX+1,MMAX+DES_MMAX
220              IF(SOLIDS_MODEL(M) == 'DEM') THEN
221     ! Number of particles for phase M
222                 rPARTS(M) = &
223                    floor((6.0d0*IC_EP_S(ICV,M)*DOM_VOL)/(PI*(D_P0(M)**3)))
224              ENDIF
225           ENDDO
226     
227     ! Total number of particles in this IC region.
228           tPARTS = sum(rPARTS)
229           IF(tPARTS == 0) RETURN
230     
231           ADJ_DIA = 2.0d0*MAX_RADIUS*lFAC
232     
233     ! Attempt to seed particle throughout the IC region
234           FIT_FAILED=.FALSE.
235           IF(IC_DES_FIT_TO_REGION(ICV)) THEN
236              IF(NO_K) THEN
237                 lDEL = (DOML(1)-ADJ_DIA)*(DOML(2)-ADJ_DIA)
238                 lDEL = (lDEL/dble(tPARTS))**(1.0/2.0)
239                 SEED_X = max(1,ceiling((DOML(1)-ADJ_DIA)/lDEL))
240                 SEED_Y = max(1,ceiling((DOML(2)-ADJ_DIA)/lDEL))
241                 SEED_Z = 1
242              ELSE
243                 lDEL = (DOML(1)-ADJ_DIA)*(DOML(2)-ADJ_DIA)*(DOML(3)-ADJ_DIA)
244                 lDEL = (lDEL/dble(tPARTS))**(1.0/3.0)
245                 SEED_X = max(1,ceiling((DOML(1)-ADJ_DIA)/lDEL))
246                 SEED_Y = max(1,ceiling((DOML(2)-ADJ_DIA)/lDEL))
247                 SEED_Z = max(1,ceiling((DOML(3)-ADJ_DIA)/lDEL))
248              ENDIF
249              FIT_FAILED=(dble(SEED_X*SEED_Y*SEED_Z) < tPARTS)
250           ENDIF
251     
252     ! Generic filling
253           IF(.NOT.IC_DES_FIT_TO_REGION(ICV) .OR. FIT_FAILED) THEN
254              SEED_X = max(1,floor((IC_END(1)-IC_START(1)-ADJ_DIA)/ADJ_DIA))
255              SEED_Y = max(1,floor((IC_END(2)-IC_START(2)-ADJ_DIA)/ADJ_DIA))
256              SEED_Z = max(1,floor((IC_END(3)-IC_START(3)-ADJ_DIA)/ADJ_DIA))
257           ENDIF
258     
259           lDX = DOML(1)/dble(SEED_X)
260           lDY = DOML(2)/dble(SEED_Y)
261           IF(DO_K) THEN
262              lDZ = DOML(3)/dble(SEED_Z)
263           ELSE
264              lDZ = 0.0d0
265           ENDIF
266     
267           xINIT = IC_START(1)+HALF*lDX
268           yINIT = IC_START(2)+HALF*lDY
269           zINIT = IC_START(3)+HALF*lDZ
270     
271           M=1
272           pCOUNT = 0
273           tCOUNT = 0
274     
275           VEL_FLUCT = SET_VEL_FLUCT(ICV,M)
276     
277           JJ_LP: DO JJ=1, SEED_Y
278              POS(2) = YINIT + (JJ-1)*lDY
279              IF(compare(POS(2),dg_ystart) .OR. compare(POS(2),dg_yend))    &
280                 POS(2) = POS(2) + SMALL_NUMBER
281     
282           KK_LP: DO KK=1, SEED_Z
283              POS(3) = ZINIT + (KK-1)*lDZ
284              IF(DO_K) THEN
285                 IF(compare(POS(3),dg_zstart) .OR. compare(POS(3),dg_zend)) &
286                    POS(3) = POS(3) + SMALL_NUMBER
287              ENDIF
288     
289           II_LP: DO II=1, SEED_X
290              POS(1) = xINIT + (II-1)*lDX
291              IF(compare(POS(1),dg_xstart) .OR. compare(POS(1),dg_xend))    &
292                 POS(1) = POS(1) + SMALL_NUMBER
293     
294     ! Exit if all particles were seeded.
295              IF(tCOUNT > int(tPARTS)) THEN
296                 EXIT JJ_LP
297     ! Find the next phase that needs to be seeded
298              ELSEIF(pCOUNT(M) > int(rPARTS(M))) THEN
299                 MM_LP: DO MM=M+1,MMAX+DES_MMAX
300                    IF(rPARTS(MM) > 0.0) THEN
301                       M=MM
302                       EXIT MM_LP
303                    ENDIF
304                 ENDDO MM_LP
305                 IF(M > MMAX+DES_MMAX) EXIT JJ_LP
306                 VEL_FLUCT = SET_VEL_FLUCT(ICV,M)
307              ENDIF
308     
309              pCOUNT(M) = pCOUNT(M) + 1
310              tCOUNT = tCOUNT + 1
311     
312     ! Keep only particles that belong to this process.
313              IF(.NOT.dg_is_ON_myPE_OWNs(IofPOS(POS(1)), &
314                 JofPOS(POS(2)),KofPOS(POS(3)))) CYCLE
315     
316     ! Bin the parcel to the fuild grid.
317              K=1
318              IF(DO_K) CALL PIC_SEARCH(K, POS(3), ZT, DIMENSION_K, KMIN2, KMAX2)
319              CALL PIC_SEARCH(J, POS(2), YN, DIMENSION_J, JMIN2, JMAX2)
320              CALL PIC_SEARCH(I, POS(1), XE, DIMENSION_I, IMIN2, IMAX2)
321     
322     ! Skip cells that return invalid IJKs.
323              IF(DEAD_CELL_AT(I,J,K)) CYCLE
324     
325     ! Skip cells that are not part of the local fuild domain.
326              IJK = FUNIJK(I,J,K)
327              IF(.NOT.FLUID_AT(IJK)) CYCLE
328     
329              IF(CARTESIAN_GRID) THEN
330                 CALL CHECK_IF_PARTICLE_OVERLAPS_STL(POS, I, J, K, SKIP)
331                 IF(SKIP) CYCLE
332              ENDIF
333     
334              PIP = PIP + 1
335              CALL PARTICLE_GROW(PIP)
336              MAX_PIP = max(PIP,MAX_PIP)
337     
338              CALL SET_NORMAL(PIP)
339     
340              IF(VEL_FLUCT) THEN
341                 VEL(1) = randVEL(pCOUNT(M),1)
342                 VEL(2) = randVEL(pCOUNT(M),2)
343                 VEL(3) = randVEL(pCOUNT(M),3)
344              ELSE
345                 VEL(1) = IC_U_s(ICV,M)
346                 VEL(2) = IC_V_s(ICV,M)
347                 VEL(3) = IC_W_s(ICV,M)
348              ENDIF
349              IF(NO_K) VEL(3) = 0.0d0
350     
351     
352              DES_POS_NEW(PIP,:) = POS(:)
353              DES_VEL_NEW(PIP,:) = VEL(:)
354              OMEGA_NEW(PIP,:) = 0.0d0
355     
356              DES_RADIUS(PIP) = D_P0(M)*HALF
357              RO_SOL(PIP) =  RO_S0(M)
358     
359              PIJK(PIP,1) = I
360              PIJK(PIP,2) = J
361              PIJK(PIP,3) = K
362              PIJK(PIP,4) = IJK
363              PIJK(PIP,5) = M
364     
365              IF(DO_OLD) THEN
366                 DES_VEL_OLD(PIP,:) = DES_VEL_NEW(PIP,:)
367                 DES_POS_OLD(PIP,:) = DES_POS_NEW(PIP,:)
368                 OMEGA_OLD(PIP,:) = ZERO
369              ENDIF
370     
371              SOLIDS_DATA(M) = SOLIDS_DATA(M) + 1.0
372     
373           ENDDO II_LP
374           ENDDO KK_LP
375           ENDDO JJ_LP
376     
377     ! Collect the data
378           CALL GLOBAL_ALL_SUM(SOLIDS_DATA)
379     
380     ! Verify that the IC region volume is not zero.
381           IF(SOLIDS_DATA(0) <= 0.0d0) THEN
382              WRITE(ERR_MSG,1000) ICV, SOLIDS_DATA(0)
383              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
384           ENDIF
385     
386     1000 FORMAT('Error 1000: Invalid IC region volume: IC=',I3,' VOL=',&
387              ES15.4,/'Please correct the mfix.dat file.')
388     
389           WRITE(ERR_MSG,2000) ICV
390           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
391     
392           DO M=MMAX+1, MMAX+DES_MMAX
393              IF(SOLIDS_DATA(M) < SMALL_NUMBER) CYCLE
394              WRITE(ERR_MSG,2010) M, int(SOLIDS_DATA(M)), IC_EP_S(ICV,M),   &
395                 (dble(SOLIDS_DATA(M))*(Pi/6.0d0)*D_P0(M)**3)/SOLIDS_DATA(0)
396              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
397           ENDDO
398     
399           IF(allocated(randVEL)) deallocate(randVEL)
400     
401           CALL FINL_ERR_MSG
402     
403           RETURN
404     
405      2000 FORMAT(/2x,'|',43('-'),'|',/2x,'| IC Region: ',I3,28x,'|',/2x,   &
406              '|',43('-'),'|',/2x,'| Phase | Number of |    EPs    |    EP',&
407              's    |',/2x,'|   ID  | Particles | Specified |   Actual  |', &
408              /2x,'|-------|',3(11('-'),'|'))
409     
410      2010 FORMAT(2x,'|  ',I3,'  |',1x,I9,1x,'|',2(1x,ES9.2,1x,'|'),/2x,    &
411              '|-------|',3(11('-'),'|'))
412     
413     
414           CONTAINS
415     
416     !......................................................................!
417     ! Function: SET_VEL_FLUCT                                              !
418     ! Purpose: Set the flag for random velocity fluctuations. If needed    !
419     ! the random velocities are calculated.                                !
420     !......................................................................!
421           LOGICAL FUNCTION SET_VEL_FLUCT(lICV, lM)
422           INTEGER, INTENT(IN) :: lICV, lM
423           DOUBLE PRECISION :: VEL_SIG
424           SET_VEL_FLUCT=(IC_Theta_M(lICV,lM) /= 0.0d0)
425           IF(SET_VEL_FLUCT) THEN
426              if(allocated(randVEL)) deallocate(randVEL)
427              allocate(randVEL(100+int(rPARTS(lM)),3))
428              VEL_SIG = sqrt(IC_Theta_M(lICV,lM))
429              CALL NOR_RNO(randVEL(:,1), IC_U_s(lICV,lM),VEL_SIG)
430              CALL NOR_RNO(randVEL(:,2), IC_V_s(lICV,lM),VEL_SIG)
431              IF(DO_K) CALL NOR_RNO(randVEL(:,3),IC_W_s(lICV,lM),VEL_SIG)
432           ENDIF
433     
434           END FUNCTION SET_VEL_FLUCT
435     
436           END SUBROUTINE GENERATE_PARTICLE_CONFIG_DEM
437     
438     
439     
440     
441     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
442     !                                                                      !
443     !  Subroutine: GENERATE_PARTICLE_CONFIG_MMPPIC                         !
444     !  Author: Rahul Garg                                 Date: 3-May-2011 !
445     !                                                                      !
446     !  Purpose: Generates particle position distribution for MPPIC.        !
447     !                                                                      !
448     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
449           SUBROUTINE GENERATE_PARTICLE_CONFIG_MPPIC(ICV)
450     
451     ! Global variables
452     !---------------------------------------------------------------------//
453     ! Number of DES solids phases.
454           use discretelement, only: DES_MMAX
455     ! Flag indicating that the IC region is defined.
456           use ic, only: IC_DEFINED
457     ! IC Region bulk density (RO_s * EP_s)
458           use ic, only: IC_ROP_s
459     ! IC Region gas volume fraction.
460           use ic, only: IC_EP_G
461     ! MPPIC specific IC region specification.
462           use ic, only: IC_PIC_CONST_NPC, IC_PIC_CONST_STATWT
463     
464           use param1, only: UNDEFINED, UNDEFINED_I
465           use param1, only: ZERO, ONE, HALF
466     ! Maximum number of IC regions and solids phases
467           use param, only: DIMENSION_IC
468           use param, only: DIM_M
469           use physprop, only: mmax
470     
471     ! The accumulated number of particles in each IJK.
472           use mpi_utility, only: GLOBAL_ALL_SUM
473     
474           use error_manager
475     
476           use run, only: solids_model
477           IMPLICIT NONE
478     
479     ! Dummy arguments
480     !---------------------------------------------------------------------//
481           INTEGER, INTENT(IN) :: ICV
482     
483     ! Local variables
484     !---------------------------------------------------------------------//
485     ! Generic loop counters
486           INTEGER :: M
487     ! Actual volume of IC region
488           DOUBLE PRECISION :: IC_VOL
489     ! Solids data in IC Region by phase:
490           DOUBLE PRECISION :: SOLIDS_DATA(0:4*DIM_M)
491     !......................................................................!
492     
493           CALL INIT_ERR_MSG("GENERATE_PARTICLE_CONFIG_MPPIC")
494     
495           WRITE(ERR_MSG,"(2/,'Generating initial parcel configuration:')")
496           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
497     
498     
499           SOLIDS_DATA = ZERO
500           CALL GET_IC_VOLUME(ICV, SOLIDS_DATA(0))
501     
502           WRITE(ERR_MSG,2000) ICV
503           CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
504     
505     ! Set up the individual solids phases.
506           DO M=MMAX+1, DES_MMAX+MMAX
507              IF(SOLIDS_MODEL(M) == 'PIC') THEN
508                 IF(IC_ROP_S(ICV,M) == ZERO) CYCLE
509     ! Seed parcels with constant stastical weight.
510                 CALL GPC_MPPIC_CONST_NPC(ICV, M, SOLIDS_DATA(0), &
511                 SOLIDS_DATA((4*M-3):(4*M)))
512              ENDIF
513           ENDDO
514     
515     ! Collect the data
516           CALL GLOBAL_ALL_SUM(SOLIDS_DATA)
517     
518     ! Verify that the IC region volume is not zero.
519           IF(SOLIDS_DATA(0) <= 0.0d0) THEN
520              WRITE(ERR_MSG,1000) ICV, SOLIDS_DATA(0)
521              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
522           ENDIF
523     
524      1000 FORMAT('Error 1000: Invalid IC region volume: IC=',I3,' VOL=',&
525              ES15.4,/'Please correct the mfix.dat file.')
526     
527     ! Report solids information for the IC region.
528           DO M=MMAX+1, DES_MMAX+MMAX
529              IF(SOLIDS_MODEL(M) == 'PIC') THEN
530                 WRITE(ERR_MSG,2010) M, int(SOLIDS_DATA(4*M-3)),&
531                    int(SOLIDS_DATA(4*M-3)*SOLIDS_DATA(4*M-2)), &
532                    SOLIDS_DATA(4*M-2), SOLIDS_DATA(4*M-1),     &
533                    SOLIDS_DATA(4*M)/SOLIDS_DATA(0)
534                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
535              ENDIF
536           ENDDO
537     
538           CALL FINL_ERR_MSG
539     
540           RETURN
541     
542      2000 FORMAT(/2x,'|',67('-'),'|',/2x,'| IC Region: ',I3,52x,'|',/2x,   &
543              '|',67('-'),'|',/2x,'| Phase | Num. Comp | Num. Real ',       &
544              '| Stastical |    EPs    |    EPs    |',/2x,'|   ID  |  ',    &
545              'Parcels  | Particles |   Weight  | Specified |   Actual  |', &
546              /2x,'|-------|',5(11('-'),'|'))
547     
548      2010 FORMAT(2x,'|  ',I3,'  |',2(1x,I9,1x,'|'),3(1x,ES9.2,1x,'|'),/2x,&
549              '|-------|',5(11('-'),'|'))
550     
551           END SUBROUTINE GENERATE_PARTICLE_CONFIG_MPPIC
552     
553     
554     
555     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
556     !                                                                      !
557     !  Subroutine: GENERATE_PARTICLE_CONFIG_MMPPIC                         !
558     !  Author: Rahul Garg                                 Date: 3-May-2011 !
559     !                                                                      !
560     !  Purpose: generates particle position distribution for MPPIC         !
561     !                                                                      !
562     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
563           SUBROUTINE GPC_MPPIC_CONST_NPC(ICV, M, IC_VOL, sDATA)
564     
565     ! Global variables
566     !---------------------------------------------------------------------//
567     ! Constant: 3.14159...
568           use constant, only: PI
569     ! Cut_cell identifier array
570           use cutcell, only: cut_cell_at
571           use cutcell, only : CARTESIAN_GRID
572           use discretelement, only: XE, YN, ZT
573     
574           use des_allocate, only: PARTICLE_GROW
575           use discretelement
576     
577     ! Flag indicating that the IC region is defined.
578           use ic, only: IC_DEFINED
579     ! IC Region bulk density (RO_s * EP_s)
580           use ic, only: IC_EP_s
581           use ic, only: IC_I_w, IC_I_e, IC_J_s, IC_J_n, IC_K_b, IC_K_t
582     
583     ! initally specified velocity field and granular temperature
584           use ic, only: IC_U_s, IC_V_s, IC_W_s, IC_Theta_M
585           use ic, only: IC_PIC_CONST_NPC
586           use ic, only: IC_PIC_CONST_STATWT
587     
588           use mfix_pic, only: des_stat_wt
589           use mpi_utility
590     
591     ! Maximum number of IC regions and solids phases
592           use param, only: DIMENSION_IC
593           use param, only: DIM_M
594           use param1, only: UNDEFINED, UNDEFINED_I, ZERO, ONE, HALF
595     
596     ! solid phase diameters and densities.
597           use physprop, only: D_p0, RO_s0
598           use run, only: solids_model
599     
600           use randomno
601           use error_manager
602           use functions
603           use run, only: solids_model
604           use des_allocate, only: PARTICLE_GROW
605     
606           IMPLICIT NONE
607     
608     ! Dummy Arguments
609     !----------------------------------------------------------------------//
610     ! Index of IC region and solids phase
611           INTEGER, INTENT(IN) :: ICV, M
612     ! Specific volume of IC region (accounts for blocked cells)
613           DOUBLE PRECISION, INTENT(IN) :: IC_VOL
614     ! Data about solids in the IC region.
615           DOUBLE PRECISION, INTENT(OUT) :: sDATA(4)
616     
617     ! Local variables
618     !----------------------------------------------------------------------//
619           DOUBLE PRECISION :: EP_SM
620     
621     ! Number of real and comp. particles in a cell.
622           DOUBLE PRECISION ::  rPARTS
623           INTEGER :: maxPARTS
624           DOUBLE PRECISION :: DOML(3), IC_START(3)
625     ! Parcel position with random
626           DOUBLE PRECISION :: POS(3)
627     ! Average velocity and standard deivation
628           DOUBLE PRECISION :: IC_VEL(3), VEL_SIG
629     ! Flag to not keep parcel.
630           LOGICAL :: SKIP
631     ! Arrasy for assigning random position and velocities
632           DOUBLE PRECISION, ALLOCATABLE :: randVEL(:,:)
633           DOUBLE PRECISION :: RAND(3)
634     ! Statistical weights
635           DOUBLE PRECISION :: STAT_WT, SUM_STAT_WT
636     ! Volume of a parcel and total solids volume
637           DOUBLE PRECISION :: sVOL, sVOL_TOT
638     ! Counter for seeded parcles.
639           INTEGER :: SEEDED
640     ! Generic loop indices and loop counters
641           INTEGER :: I, J, K, IJK, LC, LC_MAX
642     !......................................................................!
643     
644           CALL INIT_ERR_MSG("GPC_MPPIC_CONST_NPC")
645     
646           maxPARTS=25
647           allocate(randVEL(maxPARTS,3))
648     
649           IC_VEL(1) = IC_U_s(ICV,M)
650           IC_VEL(2) = IC_V_s(ICV,M)
651           IC_VEL(3) = merge(IC_W_s(ICV,M),0.0d0,DO_K)
652     
653           VEL_SIG = sqrt(IC_Theta_M(ICV,M))
654     
655     ! Volume occupied by one particle
656           sVOL = (Pi/6.0d0)*(D_P0(M)**3.d0)
657     
658           SEEDED = 0
659           sVOL_TOT = 0.0d0
660           SUM_STAT_WT = 0.0d0
661     
662           DO K = IC_K_B(ICV), IC_K_T(ICV)
663           DO J = IC_J_S(ICV), IC_J_N(ICV)
664           DO I = IC_I_W(ICV), IC_I_E(ICV)
665     
666              IF(.not.IS_ON_myPE_wobnd(I,J,K)) cycle
667              IF(DEAD_CELL_AT(I,J,K)) cycle
668     
669              IJK = FUNIJK(I,J,K)
670              IF(.not.FLUID_AT(IJK)) cycle
671     
672              rPARTS = IC_EP_s(ICV,M)*VOL(IJK)/sVOL
673     
674     ! Seed parcels with a constant stastical weight
675              IF(IC_PIC_CONST_STATWT(ICV,M) /= ZERO) THEN
676                 STAT_WT = IC_PIC_CONST_STATWT(ICV,M)
677                 LC_MAX = int(rPARTS/STAT_WT)
678     
679     ! Seed parcels with a constant number per cell
680              ELSEIF(IC_PIC_CONST_NPC(ICV,M) /= 0) THEN
681                 LC_MAX = IC_PIC_CONST_NPC(ICV,M)
682                 STAT_WT = rPARTS/dble(LC_MAX)
683                 IF(CUT_CELL_AT(IJK)) LC_MAX = max(1,int(VOL(IJK)*dble(&
684                    IC_PIC_CONST_NPC(ICV,M))/(DX(I)*DY(J)*DZ(K))))
685              ENDIF
686     
687     ! Increase particle buffer
688              IF(LC_MAX > maxPARTS) THEN
689                 maxPARTS = 2*LC_MAX
690                 if(allocated(randVEL)) deallocate(randVEL)
691                 allocate(randVEL(maxPARTS,3))
692              ENDIF
693     
694              DO LC=1, merge(2,3,NO_K)
695                 IF(VEL_SIG > ZERO) THEN
696                    CALL NOR_RNO(randVEL(1:LC_MAX,LC), IC_VEL(LC), VEL_SIG)
697                 ELSE
698                    randVEL(1:LC_MAX,LC) = IC_VEL(LC)
699                 ENDIF
700              ENDDO
701              IF(NO_K) randVEL(1:LC_MAX,3) = 0.0d0
702     
703              IC_START(1) = XE(I-1)
704              IC_START(2) = YN(J-1)
705              IC_START(3) = ZERO;  IF(DO_K) IC_START(3) = ZT(K-1)
706     
707              DOML(1) = DX(I)
708              DOML(2) = DY(J)
709              DOML(3) = ZERO;  IF(DO_K) DOML(3) = DZ(K)
710     
711              DO LC=1,LC_MAX
712     
713                 CALL RANDOM_NUMBER(RAND)
714                 POS(:) = IC_START + DOML*RAND
715     
716                 IF(CARTESIAN_GRID) THEN
717                    CALL CHECK_IF_PARCEL_OVERLAPS_STL(POS, SKIP)
718                    DO WHILE(SKIP)
719                       CALL RANDOM_NUMBER(RAND)
720                       POS(:) = IC_START + DOML*RAND
721                       CALL CHECK_IF_PARCEL_OVERLAPS_STL(POS, SKIP)
722                    ENDDO
723                 ENDIF
724     
725                 PIP = PIP + 1
726                 CALL PARTICLE_GROW(PIP)
727                 MAX_PIP = max(PIP,MAX_PIP)
728     
729                 DES_POS_NEW(PIP,:) = POS(:)
730                 DES_VEL_NEW(PIP,:) = randVEL(LC,:)
731     
732                 DES_RADIUS(PIP) = D_P0(M)*HALF
733                 RO_SOL(PIP) =  RO_S0(M)
734     
735                 PIJK(PIP,1) = I
736                 PIJK(PIP,2) = J
737                 PIJK(PIP,3) = K
738                 PIJK(PIP,4) = IJK
739                 PIJK(PIP,5) = M
740     
741                 SUM_STAT_WT = SUM_STAT_WT + STAT_WT
742     
743                 DES_STAT_WT(PIP) = STAT_WT
744                 sVOL_TOT = sVOL_TOT + sVOL*STAT_WT
745     
746                 CALL SET_NORMAL(PIP)
747     
748                 SEEDED = SEEDED + 1
749     
750              ENDDO
751     
752           ENDDO
753           ENDDO
754           ENDDO
755     
756           IF(allocated(randVEL)) deallocate(randVEL)
757     
758           sDATA(1) = dble(SEEDED)
759           sDATA(2) = SUM_STAT_WT/dble(SEEDED)
760           sDATA(3) = IC_EP_S(ICV,M)
761           sDATA(4) = sVOL_TOT
762     
763           CALL FINL_ERR_MSG
764     
765           RETURN
766           END SUBROUTINE GPC_MPPIC_CONST_NPC
767     
768     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
769     !                                                                      !
770     !  Subroutine: GET_IC_VOLUME                                           !
771     !  Author: J.Musser                                 Date: 26-Aug-2015  !
772     !                                                                      !
773     !  Purpose: Calculate the actual volume of the IC region.              !
774     !                                                                      !
775     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
776           SUBROUTINE GET_IC_VOLUME(ICV, IC_VOL)
777     
778     ! IC region index bounds
779           use ic, only: IC_I_w, IC_I_e
780           use ic, only: IC_J_s, IC_J_n
781           use ic, only: IC_K_b, IC_K_t
782     
783     ! Volume of computational cells.
784           use geometry, only: VOL
785     
786           use functions
787           use compar, only: dead_cell_at
788     
789           IMPLICIT NONE
790     
791     ! Dummy Arguments
792     !---------------------------------------------------------------------//
793     ! Index of IC region
794           INTEGER, INTENT(IN) :: ICV
795     ! Total calculated volume of IC region
796           DOUBLE PRECISION, INTENT(OUT) :: IC_VOL
797     
798     ! Local variables
799     !---------------------------------------------------------------------//
800           INTEGER :: I, J, K, IJK
801     !......................................................................!
802     
803     
804           IC_VOL = 0.0d0
805           DO K = IC_K_B(ICV), IC_K_T(ICV)
806           DO J = IC_J_S(ICV), IC_J_N(ICV)
807           DO I = IC_I_W(ICV), IC_I_E(ICV)
808     
809              IF(.NOT.IS_ON_MYPE_WOBND(I,J,K)) CYCLE
810              IF(DEAD_CELL_AT(I,J,K)) CYCLE
811     
812              IJK = FUNIJK(I,J,K)
813              IF(FLUID_AT(IJK)) IC_VOL = IC_VOL + VOL(IJK)
814     
815           ENDDO
816           ENDDO
817           ENDDO
818     
819           RETURN
820           END SUBROUTINE GET_IC_VOLUME
821     
822           END MODULE GENERATE_PARTICLES
823