File: /nfs/home/0/users/jenkins/mfix.git/model/des/generate_particle_config.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  SUBROUTINE: GENERATE_PARTICLE_CONFIG                                C
4     !                                                                      C
5     !  Purpose: Generate particle configuration based on maximum particle  C
6     !           radius and filling from top to bottom within specified     C
7     !           bounds                                                     C
8     !                                                                      C
9     !                                                                      C
10     !  Authors: Rahul Garg                                Date: 19-Mar-14  C
11     !                                                                      C
12     !                                                                      C
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
14     
15           SUBROUTINE GENERATE_PARTICLE_CONFIG
16     
17     !-----------------------------------------------
18     ! Modules
19     !-----------------------------------------------
20           USE mfix_pic, only: MPPIC
21           USE DES_LINKED_LIST_DATA, only : orig_part_list, particle
22           USE des_linked_list_funcs
23           USE cutcell, only : CARTESIAN_GRID
24           USE discretelement, only: dimn, xe, yn, zt
25           USE discretelement, only: particles, pip, DEBUG_DES
26     
27           USE mpi_utility
28     
29           USE mfix_pic, only : mppic
30     
31     ! Use the error manager for posting error messages.
32     !---------------------------------------------------------------------//
33           use error_manager
34     
35     
36           IMPLICIT NONE
37           INTEGER :: IPE
38           type(particle), pointer :: part => null()
39     
40           INTEGER :: LORIG_ALL(0:Numpes-1)
41           INTEGER :: LDEL_ALL(0:Numpes-1), LREM_ALL(0:Numpes-1)
42     
43           Logical :: write_vtp_files, indomain
44     ! Initialize the error manager.
45     
46           WRITE_VTP_FILES = .true.
47     
48           CALL INIT_ERR_MSG("Generate_Particle_Config")
49     
50           IF(MPPIC) THEN
51              CALL GENERATE_PARTICLE_CONFIG_MPPIC
52           ELSE
53              CALL GENERATE_PARTICLE_CONFIG_DEM
54           ENDIF
55     
56           LORIG_ALL = 0
57     
58     
59           IF(CARTESIAN_GRID) then
60              write(ERR_MSG, 3000)
61              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
62      3000 FORMAT('Cartesian grid detected with particle configuration,'/&
63              'deleting particles located outsidte the domain.')
64     
65              indomain = .true.
66              IF(WRITE_VTP_FILES .AND. DEBUG_DES) &
67                 CALL DBG_WRITE_PART_VTP_FILE_FROM_LIST("BEFORE_DELETION", indomain)
68     
69              IF(MPPIC) THEN
70                 write(err_msg, '(A)') 'Now Marking particles to be deleted'
71                 call flush_err_msg(header = .false., footer = .false.)
72     
73                 CALL MARK_PARTS_TOBE_DEL_PIC_STL
74              ELSE
75     
76                 write(err_msg, '(A)') 'Now Marking particles to be deleted'
77                 call flush_err_msg(header = .false., footer = .false.)
78     
79                 CALL MARK_PARTS_TOBE_DEL_DEM_STL
80              END IF
81     
82              write(ERR_MSG,"('Finished marking particles to delete.')")
83              CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
84     
85              IF(WRITE_VTP_FILES .AND. DEBUG_DES) &
86                 CALL DBG_WRITE_PART_VTP_FILE_FROM_LIST("AFTER_DELETION", indomain)
87     
88              indomain = .false.
89              IF(WRITE_VTP_FILES .AND. DEBUG_DES) &
90                 CALL DBG_WRITE_PART_VTP_FILE_FROM_LIST("DELETED", indomain)
91     
92           ENDIF
93     
94     
95           write(err_msg, '(/,70("-"),/,A)') 'Particle statistics reporting'
96     
97           call flush_err_msg(header = .false., footer = .false.)
98     
99           LDEL_ALL = 0
100           LREM_ALL = 0
101     
102           part => orig_part_list
103           DO WHILE (ASSOCIATED(part))
104              LORIG_ALL(mype) = LORIG_ALL(mype) + 1
105     
106              IF(part%INDOMAIN) then
107                 LREM_ALL(mype) = LREM_ALL(mype) + 1
108              else
109                 LDEL_ALL(mype) = LDEL_ALL(mype) + 1
110              ENDIF
111     
112              part => part%next
113           ENDDO
114     
115     
116           CALL global_all_sum(LORIG_ALL)
117           CALL global_all_sum(LREM_ALL)
118           CALL global_all_sum(LDEL_ALL)
119     
120           PIP = LREM_ALL(mype) !Setting pip on each processor equal to remaining paticles
121     
122           !total number of particles used to allocate the arrays is set equal to sum
123           !of remaining particles on all processors
124           PARTICLES  = SUM(LREM_ALL(0:numPEs-1))
125     
126           Do ipe = 0, numPEs - 1
127              write(err_msg, 1002) ipe, Lorig_all(ipe), Ldel_all(ipe), Lrem_all(ipe)
128      1002    FORMAT(/, &
129              'For proc', 2x, i5, / &
130              'Total number of particles originally seeded       : ', 2x, I15, / &
131              'Total number of particles found outside the domain: ', 2x, I15, / &
132              'Total number of particles remaining               : ', 2x, I15, /)
133     
134              CALL FLUSH_ERR_MSG (header = .false., Footer = .false.)
135     
136              if( Lrem_all(ipe) + Ldel_all(ipe) .ne. Lorig_all(ipe)) then
137     
138                 write(err_msg, 1003) ipe, Lrem_all(ipe) + Ldel_all(ipe), Lorig_all(ipe)
139                 CALL FLUSH_ERR_MSG (ABORT = .true.)
140              endif
141     
142      1003    Format('Sanity check failed for particle seeding on proc:', 2x, i5, / &
143              '# of particles deleted + # of particles remaining ',  2x, I15, / &
144              'not equal to original number of particles',  2x, I15, / &
145              'MFIX will now exit')
146     
147           ENDDO
148     
149           write(err_msg, 1004) PARTICLES
150      1004 FORMAT(/, &
151           'Total number of particles in the system: ', 2x, i15)
152     
153           CALL FLUSH_ERR_MSG (header = .false.)
154     
155           CALL FINL_ERR_MSG
156           END SUBROUTINE GENERATE_PARTICLE_CONFIG
157     
158           SUBROUTINE COPY_PARTICLE_CONFIG_FROMLISTS
159           USE DES_LINKED_LIST_Data, only : orig_part_list, particle
160           USE des_linked_list_funcs
161           use error_manager
162           USE mfix_pic
163     ! Flag: solids_model array
164           use run, only: solids_model
165     
166           USE discretelement
167           use des_allocate, only: PARTICLE_GROW
168     
169           IMPLICIT NONE
170           type(particle), pointer :: part => null()
171           integer :: count_part, phase
172     
173     ! Initialize the error manager.
174           CALL INIT_ERR_MSG("MARK_PARTS_TOBE_DEL_DEM_STL")
175     
176           CALL PARTICLE_GROW(PIP)
177     
178           part => orig_part_list
179           count_part = 0
180           DO WHILE (ASSOCIATED(part))
181              IF(part%INDOMAIN) then
182                 count_part = count_part + 1
183     
184                 des_pos_new(1:dimn, count_part) = part%position(1:dimn)
185                 des_vel_new(1:dimn, count_part) = part%velocity(1:dimn)
186     
187                 DES_RADIUS(count_part) = part%rad
188                 RO_SOL(count_part)  = part%dens
189     
190                 PIJK(count_part,1:4) = part%cell(1:4)
191                 PIJK(count_part,5) = part%M
192                 phase  = part%M
193                 IF (DO_OLD) THEN
194                    des_vel_old(:, count_part) = des_vel_new(:, count_part)
195                    des_pos_old(:, count_part) = des_pos_new(:, count_part)
196                 ENDIF
197     
198                 PEA(count_part,1) = .TRUE.
199     
200                 IF(SOLIDS_MODEL(phase).eq.'PIC') then
201                    DES_STAT_WT(count_part) = part%STATWT
202                 ELSEIF(SOLIDS_MODEL(phase).eq.'DEM') then
203                    IF (DO_OLD) OMEGA_OLD(:, count_part) = zero
204                    OMEGA_NEW(:, count_part) = zero
205                 ENDIF
206     
207              ENDIF
208              part => part%next
209     
210           ENDDO
211           if(count_part.ne.pip) then
212              write(err_msg, 1000) count_part, pip
213              call flush_err_msg(abort = .true.)
214           endif
215     
216      1000 format('Particles copied from linked lists:', 2x, i15, /, &
217           'not equal to particles earlier calculated to be inside domain', 2x, i15, / &
218           'This should not have happened, exitting')
219           CALL FINL_ERR_MSG
220     
221           CALL DEALLOC_PART_LIST(orig_part_list)
222           END SUBROUTINE COPY_PARTICLE_CONFIG_FROMLISTS
223     
224     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
225     !                                                                      !
226     ! Subroutine: MARK_PARTS_TOBE_DEL_DEM_ST                               !
227     ! Author: R. Garg                                     Date: 21-Mar-14  !
228     !                                                                      !
229     ! Purpose: Mark particles that are outside the domain using the factes !
230     !          information                                                 !
231     !                                                                      !
232     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
233           SUBROUTINE MARK_PARTS_TOBE_DEL_DEM_STL
234     
235           USE DES_LINKED_LIST_Data, only : orig_part_list, particle
236           USE calc_collision_wall
237           USE compar
238           USE cutcell, only : cut_cell_at
239           USE des_linked_list_funcs
240           USE discretelement, only: dimn
241           USE error_manager
242           USE functions
243           USE geometry
244           USE indices
245     
246           IMPLICIT NONE
247           INTEGER :: IJK
248           LOGICAL :: DELETE_PART
249           type(particle), pointer :: part => null()
250     
251     ! This call will delete the particles outside the domain. It will then
252     ! re-arrange the arrays such that the active particles are in a block.
253     ! This works only for MPPIC as the cell information is added to
254     ! particles in generate_particle_config_mppic itself. For DEM particles,
255     ! the cell information is not added in generate_particle_config but is
256     ! done only during first call to particles_in_cell.
257     
258     ! Initialize the error manager.
259           CALL INIT_ERR_MSG("MARK_PARTS_TOBE_DEL_DEM_STL")
260     
261           part => orig_part_list
262           DO WHILE (ASSOCIATED(part))
263              DELETE_PART = .false.
264              IJK = part%cell(4)
265     
266     ! Only look for particles that are in domain. Particles in dead cells
267     ! have already been flagged as outside of the domain in
268     ! BIN_PARTICLES_TO_CELL.
269              IF(part%indomain) THEN
270     
271     ! Since checking if a particle is on other side of a triangle is tricky,
272     ! for safety, initially remove all the particles in the cut-cell.
273     ! now cut-cell and actualy geometry could be off, so this might not work
274     ! very robustly.
275                 IF(FLUID_AT(IJK).and.(.not.CUT_CELL_AT(IJK))) THEN
276     
277     ! Check if any of the particles are overlapping with the walls. This
278     ! will include the normal ghost cell walls and also the cut cells.
279                    CALL CHECK_IF_PARTICLE_OVELAPS_STL(PART%POSITION(:), &
280                       PART%RAD, DELETE_PART)
281                 ELSE
282                    DELETE_PART = .TRUE.
283                 ENDIF
284              ENDIF
285              IF(DELETE_PART) PART%INDOMAIN = .FALSE.
286              PART => PART%NEXT
287            ENDDO
288     
289            RETURN
290            END SUBROUTINE MARK_PARTS_TOBE_DEL_DEM_STL
291     
292     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
293     !                                                                      !
294     ! Subroutine: MARK_PARTS_TOBE_DEL_PIC_ST                               !
295     ! Author: R. Garg                                     Date: 25-Mar-14  !
296     !                                                                      !
297     ! Purpose: Mark particles that are outside the domain using the factes !
298     !          information                                                 !
299     !                                                                      !
300     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
301           SUBROUTINE MARK_PARTS_TOBE_DEL_PIC_STL
302     
303           USE DES_LINKED_LIST_Data, only : orig_part_list, particle
304           USE des_linked_list_funcs
305           USE discretelement, only: dimn
306     
307           USE cutcell, only : cut_cell_at
308           USE indices
309           USE geometry
310           USE compar
311           USE error_manager
312           USE functions
313     
314           IMPLICIT NONE
315           INTEGER :: IJK
316           LOGICAL :: DELETE_PART
317     
318           type(particle), pointer :: part => null()
319     
320     ! This call will delete the particles outside the domain. It will then
321     ! re-arrange the arrays such that the active particles are in a block.
322     ! This works only for MPPIC as the cell information is added to
323     ! particles in generate_particle_config_mppic itself. For DEM particles,
324     ! the cell information is not added in generate_particle_config but is
325     ! done only during first call to particles_in_cell.
326     
327     ! Initialize the error manager.
328           CALL INIT_ERR_MSG("MARK_PARTS_TOBE_DEL_PIC_STL")
329           part => orig_part_list
330           DO WHILE (ASSOCIATED(part))
331              DELETE_PART = .false.
332              IJK = part%cell(4)
333              IF(part%indomain) THEN
334     
335     ! Only look for particles that are in domain. Particles in dead cells
336     ! have already been flagged as outside of the domain in
337     !  BIN_PARTICLES_TO_CELL.
338     
339                 IF(FLUID_AT(IJK)) then
340     
341                    CALL CHECK_IF_PARCEL_OVELAPS_STL(part%position(1:dimn), &
342                         DELETE_PART)
343                 ELSE
344                    DELETE_PART = .true.
345                 ENDIF
346              ENDIF
347              if(delete_part) part%indomain = .false.
348              part => part%next
349            ENDDO
350            END SUBROUTINE MARK_PARTS_TOBE_DEL_PIC_STL
351     
352     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
353     !                                                                      !
354     ! Subroutine: BIN_PARTICLES_TO_CELL                                    !
355     ! Author: R. Garg                                     Date: 21-Mar-14  !
356     !                                                                      !
357     ! Purpose: Routine to bin particles in particle lists using brute force!
358     !          technique needed to determine out of domain particles       !
359     !                                                                      !
360     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
361     
362           SUBROUTINE BIN_PARTICLES_TO_CELL
363     
364           USE DES_LINKED_LIST_Data, only : orig_part_list, particle
365           USE des_linked_list_funcs
366           USE discretelement, only: dimn, xe, yn, zt
367           USE indices
368           USE geometry
369           USE compar
370           USE mpi_utility
371           USE discretelement, only:DES_GETINDEXFROMPOS
372           USE error_manager
373           USE functions
374     
375           IMPLICIT NONE
376           type(particle), pointer :: part => null()
377     
378     ! Initialize the error manager.
379           CALL INIT_ERR_MSG("BIN_PARTICLES_TO_CELL")
380     
381           !Compute the particle cell coordinates by brute forces
382           part => orig_part_list
383           DO WHILE (ASSOCIATED(part))
384              part%cell(1) = DES_GETINDEXFROMPOS(ISTART1, IEND1, &
385              part%position(1), XE(1:size(XE,1)-1),'X','I')
386              !-1 above since xe ranges from 0:IMAX3, so size is imax3 + 1.
387              !therefore, (1:size(xe,1)) will give 1:imax3 + 1, resulting in a seg error.
388     
389              part%cell(2) = DES_GETINDEXFROMPOS(JSTART1,JEND1,part%position(2), &
390              YN(1:size(YN,1)-1),'Y','J')
391     
392              part%cell(3) = 1
393              IF(DO_K) part%cell(3) = DES_GETINDEXFROMPOS(KSTART1, &
394              KEND1,part%position(3),ZT(1:size(ZT,1)-1),'Z','K')
395     
396              IF(.NOT.DEAD_CELL_AT(part%cell(1), part%cell(2), part%cell(3))) THEN
397                 part%cell(4) = FUNIJK(part%cell(1), part%cell(2), part%cell(3))
398              ELSE
399                 part%indomain = .false.  ! flag particle as outside domain
400              ENDIF
401     
402              part => part%next !step
403           ENDDO
404     
405           CALL FINL_ERR_MSG
406           END SUBROUTINE BIN_PARTICLES_TO_CELL
407     
408     
409     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
410     !                                                                      C
411     !  SUBROUTINE: GENERATE_PARTICLE_CONFIG                                C
412     !                                                                      C
413     !  Purpose: Generate particle configuration for DEM solids for each IC C
414     !           region. Now using the particle linked lists for initial    C
415     !           build                                                      C
416     !           This routine will ultimately supersede the older rouine    C
417     !           that has not been deleted yet                              C
418     !                                                                      C
419     !  Authors: Rahul Garg                               Date: 21-Mar-2014 C
420     !                                                                      C
421     !                                                                      C
422     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
423           SUBROUTINE GENERATE_PARTICLE_CONFIG_DEM
424     
425     
426     ! Global Variables:
427     ! Runtime Flag: Generate initial particle configuation.
428           USE discretelement, only: GENER_PART_CONFIG
429     ! particle radius and density
430           USE discretelement, only: DES_RADIUS, RO_Sol
431     ! particle position new and old
432           USE discretelement, only: des_pos_new, des_pos_old
433     ! particle velocity new and old
434           USE discretelement, only: des_vel_new, des_vel_old
435     ! Simulation dimension (2D/3D)
436           USE discretelement, only: DIMN
437     ! X , Y, and Z position of cell faces of computational fluid grid
438           USE discretelement, only: XE, YN, ZT
439     ! Number of particles in the system (current)
440           USE discretelement, only:  PIP
441     ! Number of DEM solids phases.
442           USE discretelement, only: DES_MMAX
443     ! DEM solid phase diameters and densities.
444           USE discretelement, only: DES_D_p0, DES_RO_s
445     ! Computed volume of IC region for seeding
446           USE discretelement, only: VOL_IC_REGION
447     ! Number of particles estimated to be seeded, per phase in each IC region
448           USE discretelement, only: PART_MPHASE_BYIC
449     ! Number of particles actually seeded, per phase in each IC region
450           USE discretelement, only: REALPART_MPHASE_BYIC
451     
452     ! Number of particles to read from input file.
453           USE discretelement, only: PARTICLES
454     
455           USE discretelement, only: PEA
456     ! Constant: 3.14159...
457     
458           USE discretelement, only:DES_GETINDEXFROMPOS
459           USE constant, only: PI
460     ! Flag indicating that the IC region is defined.
461           USE ic, only: IC_DEFINED
462     ! IC Region bulk density (RO_s * EP_s)
463           USE ic, only: IC_ROP_s
464     ! IC Region solids volume fraction.
465           USE ic, only: IC_EP_S
466     ! IC Region gas volume fraction.
467           USE ic, only: IC_EP_G
468     ! min and max physical co-ordinates of IC regions in each direction
469           USE ic, only: IC_X_w, IC_X_e, IC_Y_s, IC_Y_n, IC_Z_b, IC_Z_t
470     ! initally specified velocity field and granular temperature
471           USE ic, only: IC_U_s, IC_V_s, IC_W_s, IC_Theta_M
472     ! Flag to extend the lattice distribution in a given IC to available area
473           Use ic, only: IC_DES_FIT_TO_REGION
474     ! Parameter for detecting unspecified values, zero, and one
475           USE param1, only: UNDEFINED, UNDEFINED_I, ZERO, ONE, Half
476     ! Parameter for small and large numbers
477           USE param1, only: SMALL_NUMBER, LARGE_NUMBER
478     ! Maximum number of initial conditions
479           USE param, only: DIMENSION_IC
480     ! first and last index of the physical cells in regular MFIX grid
481           USe compar, only: istart1, iend1, jstart1, jend1, kstart1, kend1
482     ! to access random number generator subroutines
483           USE randomno
484           USE mpi_utility
485     
486           use error_manager
487     
488     
489     ! direction wise spans of the domain and grid spacing in each direction
490           Use geometry, only: xlength, ylength, zlength
491     
492           USE DES_LINKED_LIST_Data, only : orig_part_list, particle
493           USE des_linked_list_funcs
494           USE functions
495           IMPLICIT NONE
496     !-----------------------------------------------
497     ! Local variables
498     !-----------------------------------------------
499           INTEGER :: M, ICV, I,J, K, idim, IJK
500           INTEGER :: lproc_parcount, pcount_byic_byphase(dimension_ic, DES_MMAX)
501           INTEGER :: seed_x, seed_y, seed_z
502           INTEGER :: TOTAL_PARTS_IC, TMP_PART_COUNT_INTHIS_IC
503           double precision lmax_dia,lfac,xp,yp,zp
504           double precision :: XSTART_IC, YSTART_IC, ZSTART_IC, adj_dia, ep_sm
505           double precision :: XEND_IC, YEND_IC, ZEND_IC
506           double precision :: xinit, yinit, zinit, ymax , doml(3)
507           !max point where the particle center cud be placed
508           double precision ::  max_ic_pt(3)
509           !factor to re-fit the configuration to the IC region size
510           double precision :: refit_fac(3), volijk
511           !particles min and max positions
512           double precision :: PART_CEN_MIN(DIMN), PART_CEN_MAX(DIMN)
513     
514           !Particle mean velocity and standard deviation
515           double precision :: vel_mean(3), vel_sig(3)
516     
517           double precision, dimension(:,:), allocatable :: pvel_temp
518     
519           double precision :: rad, dens, position(dimn), velocity(dimn), wtp
520     
521           type(particle), pointer :: part_list_byic => NULL(), part => NULL(), part_old => NULL()
522     
523           CALL INIT_ERR_MSG("GENERATE_PARTICLE_CONFIG_DEM")
524     
525     ! initializing particle count
526           lproc_parcount = 0
527     
528           PCOUNT_BYIC_BYPHASE(:,:)  = 0
529     ! setting particle seed spacing grid to be slightly greater than
530     ! the maximum particle diameter. seed at ~particle radii
531           lfac = 1.05d0
532     
533     
534           if(.not.allocated(part_mphase_byic)) &
535                ALLOCATE(PART_MPHASE_BYIC(DIMENSION_IC, DES_MMAX))
536           IF(.not.allocated(REALPART_MPHASE_BYIC)) &
537                ALLOCATE(REALPART_MPHASE_BYIC(DIMENSION_IC, DES_MMAX))
538     
539           write(err_msg, 2015)
540           CALL FLUSH_ERR_MSG(header = .false., FOOTER = .false.)
541     
542      2015 FORMAT('IC region wise report on particle initialization')
543     
544           WTP = 1.d0
545     
546           IC_LOOP : DO ICV = 1, DIMENSION_IC
547           IF(associated(part_list_byic)) Nullify(part_list_byic)
548     
549              PART_MPHASE_BYIC(ICV,:) = 0
550              REALPART_MPHASE_BYIC(ICV,:) = 0
551              IF (.not.IC_DEFINED(ICV)) cycle
552     
553              PART_CEN_MIN(:) = LARGE_NUMBER
554              PART_CEN_MAX(:) = SMALL_NUMBER
555     
556              TOTAL_PARTS_IC = 0
557     
558              LMAX_DIA = SMALL_NUMBER
559              DO M = 1, DES_MMAX
560                 DOML(1) = IC_X_E(ICV) - IC_X_W(ICV)
561                 DOML(2) = IC_Y_N(ICV) - IC_Y_S(ICV)
562                 DOML(3) = IC_Z_T(ICV) - IC_Z_B(ICV)
563     
564                 IF(NO_K) DOML(3) = ZLENGTH
565     
566                 VOLIJK = DOML(1)*DOML(2)*DOML(3)
567     
568                 PART_MPHASE_BYIC(ICV, M) = FLOOR((6.D0*IC_EP_S(ICV, M)*VOLIJK)/&
569                      (PI*(DES_D_P0(M)**3)))
570                 ! setting a local value of maximum diameter for each IC region
571                 if(PART_MPHASE_BYIC(ICV,M).gt.0) then
572                    total_parts_ic = total_parts_ic + PART_MPHASE_BYIC(ICV,M)
573                    LMAX_DIA = MAX(LMAX_DIA, DES_D_P0(M))
574                 endif
575              ENDDO
576     
577              IF(total_parts_ic.eq.0) cycle IC_LOOP
578     
579              WRITE(ERR_MSG,2016) ICV
580     
581     2016     FORMAT(/1X,70('-')/, 5x, &
582                   'IC number         = ', I4)
583     
584              CALL FLUSH_ERR_MSG(HEADER = .false., Footer = .false.)
585     
586     
587              XSTART_IC = IC_X_W(ICV)
588              YSTART_IC = IC_Y_S(ICV)
589              ZSTART_IC = IC_Z_B(ICV)
590              XEND_IC   = IC_X_E(ICV)
591              YEND_IC   = IC_Y_N(ICV)
592              ZEND_IC   = IC_Z_T(ICV)
593     
594              !Max location of any particle in this IC
595              MAX_IC_PT(1) = XEND_IC - 0.5d0*LMAX_DIA*LFAC
596              MAX_IC_PT(2) = YEND_IC - 0.5d0*LMAX_DIA*LFAC
597              MAX_IC_PT(3) = ZEND_IC - 0.5d0*LMAX_DIA*LFAC
598              XINIT = XSTART_IC
599              YINIT = YSTART_IC
600              ZINIT = ZSTART_IC
601     
602              DO M = 1, DES_MMAX
603                 if(PART_MPHASE_BYIC(ICV,M).eq.0) cycle
604     
605                 VEL_MEAN(1) = IC_U_S(ICV, M)
606                 VEL_MEAN(2) = IC_V_S(ICV, M)
607                 IF(DO_K) VEL_MEAN(3) = IC_W_S(ICV, M)
608                 !granular temp is defined as (Variance uprime + Variance vprime + Variance wprime)/3
609                 !assuming equal energy in each direction
610                 !Variance uprime  = IC_Theta
611                 !Stdev (or sigma) = sqrt(Variance)
612     
613                 VEL_SIG(:) = sqrt(IC_Theta_M(ICV, M))
614                 write(ERR_MSG,2022) M,  &
615                      vel_mean(:), IC_theta_m(ICV, M), vel_sig(:)
616                 CALL FLUSH_ERR_MSG(HEADER = .false., Footer = .false.)
617     
618     
619      2022       FORMAT(1X,70('.'),/5x, &
620                 'PHASE INDEX, M                              =  ', I5,2X, /5x, &
621                 'INITIALIZING SOLIDS VELOCITY FIELD', /5x, &
622                 'Mean velocity direction wise                =  ', 3(G15.8,2X), /5x, &
623                 'Use specified initial granular temperature  =  ', (G15.8,2X), /5x, &
624                 'Velocity standard deviation direction wise  =  ', 3(G15.8,2X))
625     
626                 allocate(pvel_temp( PART_MPHASE_BYIC(ICV,M) , DIMN))
627     
628                 do IDIM = 1, DIMN
629                    pvel_temp(:, idim) = vel_mean(idim)
630     
631                    IF(vel_sig(idim).gt.zero) then
632                       CALL nor_rno(pvel_temp(1:PART_MPHASE_BYIC(ICV,M),IDIM), &
633                            vel_mean(idim),vel_sig(idim))
634                    ENDIF
635                 ENDDO
636     
637                 SEED_X = 1
638                 SEED_Y = 1
639                 SEED_Z = 1
640     
641                 ADJ_DIA = LFAC*DES_D_P0(M)
642     
643                 SEED_X = FLOOR((XEND_IC - XINIT)/ADJ_DIA)
644                 SEED_Y = FLOOR((YEND_IC - YINIT)/ADJ_DIA)
645                 SEED_Z = FLOOR((ZEND_IC - ZINIT)/ADJ_DIA)
646                 !write(*,*) 'adj_dia = ', adj_dia, lfac, lmax_dia
647                 !write(*,*) 'seedx  = ', seed_x, seed_y, seed_z
648                 if(NO_K) seed_z = 1
649     
650                 outerloop :  DO  J = 1, SEED_Y
651                    DO  K = 1, SEED_Z
652                       DO  I = 1, SEED_X
653                          XP = XINIT + I*ADJ_DIA - DES_D_P0(M)*0.5D0
654                          YP = YINIT + J*ADJ_DIA - DES_D_P0(M)*0.5D0
655                          ZP = ZINIT + K*ADJ_DIA - DES_D_P0(M)*0.5D0
656     
657                          TMP_PART_COUNT_INTHIS_IC = (PCOUNT_BYIC_BYPHASE(ICV,M)) + 1
658     
659                          IF(TMP_PART_COUNT_INTHIS_IC.GT.PART_MPHASE_BYIC(ICV,M)) EXIT outerloop
660     
661                          !Associate this new particle to the solid phase id based on the map_to_proc defined
662                          !earlier
663     
664     
665                          PCOUNT_BYIC_BYPHASE(ICV,M)  = PCOUNT_BYIC_BYPHASE(ICV,M) + 1
666     
667                          !computational FLUID grid (xe,yn, zt)
668                          !IF ( XP.GE.XE(ISTART1-1) .AND.XP.LT.XE(IEND1) &
669                          !     .AND.YP.GE.YN(JSTART1-1) .AND.YP.LT.YN(JEND1) &
670                          !     .AND.ZP.GE.ZT(KSTART1-1) .AND.ZP.LT.ZT(KEND1)) THEN
671     
672                          RAD = DES_D_P0(M)*HALF
673                          DENS  =  DES_RO_S(M)
674                          POSITION(1) = XP
675                          POSITION(2) = YP
676                          IF(DO_K) POSITION(3) = ZP
677                          VELOCITY(1:DIMN) = PVEL_TEMP(PCOUNT_BYIC_BYPHASE(ICV,M),1:DIMN)
678     
679                          CALL GEN_AND_ADD_TO_PART_LIST(PART_LIST_BYIC, M, POSITION(1:DIMN), &
680                          VELOCITY(1:DIMN), RAD, DENS, WTP)
681     
682     
683                          PART_CEN_MIN(1)  = MIN(XP , PART_CEN_MIN(1))
684                          PART_CEN_MIN(2)  = MIN(YP , PART_CEN_MIN(2))
685     
686                          PART_CEN_MAX(1)  = MAX(XP , PART_CEN_MAX(1))
687                          PART_CEN_MAX(2)  = MAX(YP , PART_CEN_MAX(2))
688                          IF(DO_K) THEN
689                             PART_CEN_MIN(3)  = MIN(ZP, PART_CEN_MIN(3))
690                             PART_CEN_MAX(3)  = MAX(ZP, PART_CEN_MAX(3))
691                          ENDIF
692     
693                          YMAX = YP + DES_D_P0(M)*0.5D0
694     
695                       end DO
696                    end DO
697                 end DO outerloop
698     
699                 YINIT = YMAX
700     
701                 DEALLOCATE(pvel_temp)
702     
703                 EP_SM = IC_EP_S(ICV,M)
704                 WRITE(ERR_MSG,2017) EP_SM
705     
706                 CALL FLUSH_ERR_MSG(HEADER = .false., Footer = .false.)
707     
708      2017       FORMAT(5x, 'Solids vol fraction for M   =  ', (G15.8,2X))
709     
710     
711              end DO                 ! DO M = 1, DES_MMAX
712     
713              refit_fac = 1.d0
714              IF(IC_DES_FIT_TO_REGION(ICV)) THEN
715     
716                 DO IDIM = 1, DIMN
717                    IF((PART_CEN_MAX(IDIM)-PART_CEN_MIN(IDIM).GT.LMAX_DIA).AND. &
718                         (MAX_IC_PT(IDIM) - PART_CEN_MAX(IDIM).GT.LMAX_DIA)) THEN
719     
720                       REFIT_FAC(IDIM)  = MAX_IC_PT(IDIM)/PART_CEN_MAX(IDIM)
721                       !write(*,*) ' REFI, IDIM =', IDIM, REFIT_FAC(IDIM)
722                    END IF
723                 END DO
724     
725                 write(err_msg, 2020) ICV, refit_fac(1:dimn)
726                 CALL FLUSH_ERR_MSG(HEADER = .false., Footer = .false.)
727     
728      2020       Format(/5x, 'Refitting to box for IC', I4, /5x,   &
729                 'Refitting factors (1-DIMN): ', 3(2x,g17.8))
730     
731     
732              END IF              !DES_IC_FIT_TO_REGION
733     
734              part => part_list_byic
735              DO WHILE (ASSOCIATED(part))
736                 !
737                 part_old => part
738                 part => part%next
739                 !Increment part here itself since the later
740                 !part of the logic under this loop will be executed
741                 !processor wise only
742                 do idim = 1, dimn
743                    part_old%position(idim) = &
744                         part_old%position(idim)*REFIT_FAC(IDIM)
745                 ENDDO
746     
747                 I = DES_GETINDEXFROMPOS(ISTART1, IEND1, &
748                 part_old%position(1), XE(1:size(XE,1)-1),'X','I')
749                 !-1 above since xe ranges from 0:IMAX3, so size is imax3 + 1.
750                 !therefore, (1:size(xe,1)) will give 1:imax3 + 1, resulting in a seg error.
751     
752                 J = DES_GETINDEXFROMPOS(JSTART1,JEND1, &
753                 part_old%position(2),YN(1:size(YN,1)-1),'Y','J')
754     
755                 K = 1
756                 IF(DO_K) K = DES_GETINDEXFROMPOS(KSTART1, &
757                 KEND1,part_old%position(3),ZT(1:size(ZT,1)-1),'Z','K')
758     
759                 IF(.not.IS_ON_myPE_wobnd(I,J,K)) cycle
760     
761                 IF(DEAD_CELL_AT(I,J,K)) cycle
762     
763                 IJK = FUNIJK(I, J, K)
764                 IF(.not.FLUID_AT(IJK)) cycle
765     
766                 CALL GEN_AND_ADD_TO_PART_LIST(ORIG_PART_LIST, part_old%M, &
767                 part_old%POSITION(1:DIMN),part_old%VELOCITY(1:DIMN), &
768                 part_old%RAD, &
769                 part_old%DENS, part_old%statwt)
770     
771                 !find its linearized cell ID which is used in
772                 !flagging out of domain particles
773     
774                 orig_part_list%cell(1) = I
775                 orig_part_list%cell(2) = J
776                 orig_part_list%cell(3) = K
777                 orig_part_list%cell(4) = FUNIJK(I, J, K)
778     
779                 REALPART_MPHASE_BYIC(ICV, part_old%M) =  &
780                      REALPART_MPHASE_BYIC(ICV, part_old%M) + 1
781              ENDDO
782     
783              !delete this list from the memory
784              CALL DEALLOC_PART_LIST(part_list_byic)
785     
786     
787              CALL global_all_sum(REALPART_MPHASE_BYIC(ICV, 1:DES_MMAX))
788     
789              DO M = 1, DES_MMAX
790     
791                 WRITE(ERR_MSG,'(//, 70("-"),/5x, A, I5,/5x,A, I15, /5x, A, I15)') &
792                      'For Phase M:              ', M,  &
793                      '# of particles estimated      :',  PART_MPHASE_BYIC(ICV, M), &
794                      '# of particles acutally seeded:', INT(REALPART_MPHASE_BYIC(ICV, M))
795     
796     
797                 CALL FLUSH_ERR_MSG(header = .false. ,footer = .false.)
798     
799              END DO
800     
801     
802           end DO IC_LOOP
803     
804     
805           CALL FINL_ERR_MSG
806     
807           RETURN
808         END SUBROUTINE GENERATE_PARTICLE_CONFIG_DEM
809     
810     
811     
812     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
813     !                                                                      C
814     !  Subroutine: GENERATE_PARTICLE_CONFIG_MMPPIC                         C
815     !  Purpose: generates particle position distribution for MPPIC         C
816     !  Author: Rahul Garg                                 Date: 3-May-2011 C
817     !                                                                      C
818     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
819           SUBROUTINE GENERATE_PARTICLE_CONFIG_MPPIC
820     
821     !-----------------------------------------------
822     ! Modules
823     !-----------------------------------------------
824           USE DES_LINKED_LIST_DATA, only : orig_part_list, particle
825           USE des_linked_list_funcs
826           USE cutcell, only : CARTESIAN_GRID
827           USE discretelement, only: dimn, xe, yn, zt
828           USE discretelement, only: particles, pip
829     
830           USE discretelement, only:DES_GETINDEXFROMPOS
831     ! Number of DES solids phases.
832           USE discretelement, only: DES_MMAX
833     
834     ! DEM solid phase diameters and densities.
835           USE discretelement, only: DES_D_p0, DES_RO_s
836     ! Number of particles seeded, per phase in each IC region
837           USE discretelement, only: PART_MPHASE_BYIC
838     ! Number of implied real particles by phase by ic region
839           USE discretelement, only: REALPART_MPHASE_BYIC
840     
841     ! Global Number of particles seeded per phase
842           USE discretelement, only: PART_MPHASE
843     ! Global Number of implied real particles
844           USE discretelement, only: REALPART_MPHASE
845     
846           USE mfix_pic, only : mppic
847     
848     ! Implied total number of physical particles
849           USE mfix_pic, only: rnp_pic
850     ! Total number of computational particles/parcels
851           USE mfix_pic, only: cnp_pic
852     
853           USE mfix_pic, only: cnp_array
854     
855     ! Flag indicating that the IC region is defined.
856           USE ic, only: IC_DEFINED
857     ! IC Region bulk density (RO_s * EP_s)
858           USE ic, only: IC_ROP_s
859     ! IC Region gas volume fraction.
860           USE ic, only: IC_EP_G
861     ! IC Region solid volume fraction.
862           USE ic, only: IC_EP_S
863     
864           USE ic, only: IC_X_w, IC_X_e, IC_Y_s, IC_Y_n, IC_Z_b, IC_Z_t
865           USE ic, only: IC_I_w, IC_I_e, IC_J_s, IC_J_n, IC_K_b, IC_K_t
866     
867     ! initally specified velocity field and granular temperature
868           USE ic, only: IC_U_s, IC_V_s, IC_W_s, IC_Theta_M
869     
870     ! Cut_cell identifier array
871           USE cutcell, only: cut_cell_at
872     
873           USE param1, only: UNDEFINED, UNDEFINED_I, ZERO, ONE
874     
875     ! MPPIC specific IC region specification.
876           USE ic, only: IC_PIC_CONST_NPC, IC_PIC_CONST_STATWT
877     
878     ! Maximum number of IC regions and solids phases
879           USE param, only: DIMENSION_IC
880     
881     ! Constant: 3.14159...
882           USE constant, only: PI
883     
884           USE mpi_utility
885     
886           USE randomno
887           USE error_manager
888           USE functions
889     
890           IMPLICIT NONE
891     !-----------------------------------------------
892     ! Local variables
893     !-----------------------------------------------
894           DOUBLE PRECISION :: EP_SM
895     ! Temp logical variables for checking constant npc and statwt specification
896           LOGICAL :: CONST_NPC, CONST_STATWT
897     
898           INTEGER :: I, J, K, IJK, ICV, M, COUNT_IC, IDIM
899     
900     ! Number of real and comp. particles in a cell.
901           DOUBLE PRECISION ::  REAL_PARTS(DIM_M)
902           INTEGER :: COMP_PARTS(DIM_M), IPCOUNT
903           DOUBLE PRECISION :: DOML(3), CORD_START(3) , CORD_END(3)
904           DOUBLE PRECISION :: STAT_WT
905     
906           DOUBLE PRECISION :: RAD, DENS, POSITION(DIMN), VELOCITY(DIMN)
907           DOUBLE PRECISION :: VOLIJK, VOLIJK_UNCUT
908           DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: RANDPOS
909           double precision :: vel_mean(3), vel_sig(3)
910     
911     
912           double precision, dimension(:,:), allocatable :: pvel_temp
913     
914           type(particle), pointer :: part_list_byic
915     
916           CALL INIT_ERR_MSG("GENERATE_PARTICLE_CONFIG_MPPIC")
917           !ALLOCATE(CNP_ARRAY(DIM3_TEMP, DES_MMAX))
918     
919           ALLOCATE(RNP_PIC(DES_MMAX))
920           ALLOCATE(CNP_PIC(DES_MMAX))
921     
922           if(.not.allocated(part_mphase_byic)) &
923                ALLOCATE(PART_MPHASE_BYIC(DIMENSION_IC, DES_MMAX))
924           IF(.not.allocated(REALPART_MPHASE_BYIC)) &
925                ALLOCATE(REALPART_MPHASE_BYIC(DIMENSION_IC, DES_MMAX))
926     
927     
928           RNP_PIC(:) = ZERO
929           CNP_PIC(:) = ZERO
930     
931           COUNT_IC = 0
932     
933           write(ERR_MSG,'(/,A,/,A,/)') 'Generating particle configuration for MPPIC model', &
934                'Per IC, per Solid phase reporting follows'
935           CALL FLUSH_ERR_MSG(header = .false. ,footer = .false.)
936     
937           DO ICV = 1, DIMENSION_IC
938              PART_MPHASE_BYIC(ICV,:) = 0
939              REALPART_MPHASE_BYIC(ICV,:) = 0
940     
941              IF(associated(part_list_byic)) Nullify(part_list_byic)
942     
943              IF(.not.ic_defined(icv)) cycle
944     
945              IF (IC_EP_G(ICV).eq.ONE) cycle
946     
947              COUNT_IC = COUNT_IC + 1
948     
949     
950              MLOOP: DO M = 1, DES_MMAX
951                 COMP_PARTS(M) = 0
952                 REAL_PARTS(M) = zero
953     
954                 EP_SM = IC_ROP_S(ICV,M)/DES_RO_s(M)
955                 IF(EP_SM.eq.zero) cycle
956                 write(ERR_MSG,'(/,70("-"), /, 5x, A, i4, 2x, A, i2)') 'Generating for IC#', ICV, ' and phase', M
957                 CALL FLUSH_ERR_MSG(header = .false., footer = .false.)
958     
959                 WRITE(ERR_MSG,2017) DES_D_P0(M), EP_SM
960                 CALL FLUSH_ERR_MSG(HEADER = .false., FOOTER = .false.)
961     
962     2017        FORMAT(/, 5x, &
963                      'Diameter                    =  ', (ES15.7,2X), /5x, &
964                      'Solids vol fraction for M   =  ', (ES15.7,2X), /5x)
965     
966                 VEL_MEAN = ZERO
967                 VEL_SIG = ZERO
968                 VEL_MEAN(1) = IC_U_S(ICV, M)
969                 VEL_MEAN(2) = IC_V_S(ICV, M)
970                 IF(DO_K) VEL_MEAN(3) = IC_W_S(ICV, M)
971                 !granular temp is defined as (Variance uprime + Variance vprime + Variance wprime)/3
972                 !assuming equal energy in each direction
973                 !Variance uprime  = IC_Theta
974                 !Stdev (or sigma) = sqrt(Variance)
975     
976                 VEL_SIG(:) = sqrt(IC_Theta_M(ICV, M))
977     
978                 write(ERR_MSG,2022) &
979                      vel_mean(:), IC_theta_m(ICV, M), vel_sig(:)
980                 CALL FLUSH_ERR_MSG(HEADER = .false., Footer = .false.)
981     
982     
983     2022        FORMAT(/,10x, 'Initializing solids velocity field', /10x, &
984                      'Mean velocity direction wise                =  ', 3(G15.8,2X), /10x, &
985                      'Use specified initial granular temperature  =  ', (G15.8,2X), /10x, &
986                      'Velocity standard deviation direction wise  =  ', 3(G15.8,2X))
987     
988                 CONST_NPC    = (IC_PIC_CONST_NPC   (ICV, M) .ne. 0) &
989                      .AND. (EP_SM.gt.Zero)
990     
991                 CONST_STATWT = (IC_PIC_CONST_STATWT(ICV, M) .ne. ZERO  ) &
992                      .AND. (EP_SM.gt.Zero)
993     
994                 IF(CONST_NPC) then
995                    !Estimate the number of parcels
996                    COMP_PARTS(M) = IC_PIC_CONST_NPC(ICV,M)*(IC_K_T(ICV) - IC_K_B(ICV) +1)* &
997                         (IC_J_N(ICV) - IC_J_S(ICV) + 1)* (IC_I_E(ICV) - IC_I_W(ICV) + 1)
998     
999                    WRITE(ERR_MSG, '(/, 10x, A, 2x, i5,/10x, A, I20)')  &
1000                         'Constant Number of particles per cell specified         = ', IC_PIC_CONST_NPC(ICV, M), &
1001                         'Number of estimated parcels  =                          = ', comp_parts(m)
1002     
1003                    CALL FLUSH_ERR_MSG(header = .false., FOOTER = .false.)
1004                    comp_parts(m) = 0
1005     
1006                    DO K = IC_K_B(ICV), IC_K_T(ICV)
1007                       DO J = IC_J_S(ICV), IC_J_N(ICV)
1008                          DO I = IC_I_W(ICV), IC_I_E(ICV)
1009     
1010                             IF(.not.IS_ON_myPE_wobnd(I,J,K)) cycle
1011                             IF(DEAD_CELL_AT(I,J,K)) cycle
1012                             IJK = FUNIJK(I,J,K)
1013     
1014                             IF(.not.FLUID_AT(IJK)) cycle
1015     
1016                             CORD_START(1) = XE(I-1)
1017                             CORD_START(2) = YN(J-1)
1018                             CORD_START(3) = ZT(K-1)
1019                             DOML(1) = DX(I)
1020                             DOML(2) = DY(J)
1021                             IF(DO_K) DOML(3) = DZ(K)
1022     
1023                             !VOL(IJK) is calculated by set_geometry1, which is
1024                             !called in mfix.f after calling get_data. At this
1025                             !stage, VOL(IJK) = zero
1026                             !VOLIJK = VOL(IJK)
1027                             VOLIJK       = DX(I)*DY(J)*DZ(K)
1028                             VOLIJK_UNCUT = DX(I)*DY(J)*DZ(K)
1029     
1030                             REAL_PARTS(M) = 0.
1031                             COMP_PARTS(M) = 0
1032     
1033                             REAL_PARTS(M) = 6.d0*EP_SM*VOLIJK/ &
1034                                  (PI*(Des_D_P0(M)**3.d0))
1035     
1036     
1037                             COMP_PARTS(M) = IC_PIC_CONST_NPC(ICV,M)
1038                             IF(CUT_CELL_AT(IJK)) COMP_PARTS(M) = &
1039                                  MAX(1,INT(VOLIJK*real(COMP_PARTS(M))/VOLIJK_UNCUT))
1040     
1041                             STAT_WT = REAL_PARTS(M)/REAL(COMP_PARTS(M))
1042     
1043     
1044                             ALLOCATE(RANDPOS    (COMP_PARTS(M), DIMN))
1045                             ALLOCATE(PVEL_TEMP  (COMP_PARTS(M), DIMN))
1046     
1047     
1048                             DO IDIM = 1, DIMN
1049                                CALL UNI_RNO(RANDPOS(1:COMP_PARTS(M), IDIM))
1050                                PVEL_TEMP(:, IDIM) = VEL_MEAN(IDIM)
1051                                IF(VEL_SIG(IDIM).GT.ZERO) THEN
1052                                   CALL NOR_RNO(PVEL_TEMP(1:COMP_PARTS(M),IDIM), &
1053                                        VEL_MEAN(IDIM),VEL_SIG(IDIM))
1054                                ENDIF
1055                             ENDDO
1056     
1057                             DO IPCOUNT = 1, COMP_PARTS(M)
1058     
1059                                RNP_PIC(M) = RNP_PIC(M) + STAT_WT
1060                                CNP_PIC(M) = CNP_PIC(M) + 1
1061                                PART_MPHASE_BYIC(ICV, M) =  PART_MPHASE_BYIC(ICV, M) &
1062                                     + 1
1063                                REALPART_MPHASE_BYIC(ICV, M) =  REALPART_MPHASE_BYIC(ICV, M) &
1064                                     + STAT_WT
1065     
1066                                RAD = DES_D_P0(M)*HALF
1067                                DENS  =  DES_RO_S(M)
1068                                POSITION(:) = CORD_START(:) + RANDPOS(IPCOUNT, :)*DOML(:)
1069                                VELOCITY(1:DIMN) = PVEL_TEMP(IPCOUNT,1:DIMN)
1070     
1071                                CALL GEN_AND_ADD_TO_PART_LIST(PART_LIST_BYIC, M, POSITION(1:DIMN), &
1072                                     VELOCITY(1:DIMN), RAD, DENS, STAT_WT)
1073                                PART_LIST_BYIC%cell(1) = I
1074                                PART_LIST_BYIC%cell(2) = J
1075                                PART_LIST_BYIC%cell(3) = K
1076                                PART_LIST_BYIC%cell(4) = IJK
1077                             ENDDO
1078                             deallocate(pvel_temp)
1079                             deallocate(randpos)
1080                          ENDDO
1081                       ENDDO
1082                    ENDDO
1083     
1084     
1085                 ELSEIF(CONST_STATWT) THEN
1086     
1087     
1088                    CORD_START(1) = IC_X_W(ICV)
1089                    CORD_START(2) = IC_Y_S(ICV)
1090                    CORD_START(3) = IC_Z_B(ICV)
1091                    CORD_END(1)   = IC_X_E(ICV)
1092                    CORD_END(2)   = IC_Y_N(ICV)
1093                    CORD_END(3)   = IC_Z_T(ICV)
1094     
1095                    DOML(:) = CORD_END(:) - CORD_START(:)
1096     
1097                    IF(NO_K) DOML(3) = DZ(1)
1098     
1099                    VOLIJK = DOML(1)*DOML(2)*DOML(3)
1100     
1101                    REAL_PARTS(M) = 0.
1102                    COMP_PARTS(M) = 0
1103     
1104                    REAL_PARTS(M) = 6.d0*EP_SM*VOLIJK/ &
1105                         (PI*(Des_D_P0(M)**3.d0))
1106     
1107                    COMP_PARTS(M) = MAX(1, &
1108                         INT(REAL_PARTS(M)/REAL(IC_PIC_CONST_STATWT(ICV,M))))
1109     
1110                    ! although the weight was specified in the input file, the integer
1111                    ! number of CP's requires recalculating statistical weight. This
1112                    ! slightly different statistical weight will ensure that the initial
1113                    ! volume fraction is as inputted. If the input statwt_pic is used,
1114                    ! then the initial volume fraction might be slightly less than the
1115                    ! input initial volume fraction
1116                    STAT_WT = REAL_PARTS(M)/REAL(COMP_PARTS(M))
1117     
1118                    WRITE(ERR_MSG, '(/10x, A, 2x, ES15.7, /10x, A, ES15.7, /10x, A, i20)') &
1119                         'Constant statistical weight specified                  = ', IC_PIC_CONST_STATWT(ICV, M), &
1120                         'Statistical weight computed (could be a little different from specified value)  = ', STAT_WT, &
1121                         'Number of estimated parcels                             = ', comp_parts(m)
1122     
1123                    CALL FLUSH_ERR_MSG(header = .false., FOOTER = .false.)
1124     
1125                    ALLOCATE(RANDPOS    (COMP_PARTS(M), DIMN))
1126                    ALLOCATE(PVEL_TEMP  (COMP_PARTS(M), DIMN))
1127     
1128     
1129                    DO IDIM = 1, DIMN
1130                       CALL UNI_RNO(RANDPOS(1:COMP_PARTS(M), IDIM))
1131                       PVEL_TEMP(:, IDIM) = VEL_MEAN(IDIM)
1132                       IF(VEL_SIG(IDIM).GT.ZERO) THEN
1133                          CALL NOR_RNO(PVEL_TEMP(1:COMP_PARTS(M),IDIM), &
1134                               VEL_MEAN(IDIM),VEL_SIG(IDIM))
1135                       ENDIF
1136                    ENDDO
1137     
1138                    DO IPCOUNT = 1, COMP_PARTS(M)
1139     
1140                       POSITION(:) = CORD_START(:) + RANDPOS(IPCOUNT, :)*DOML(:)
1141     
1142                       I = DES_GETINDEXFROMPOS(ISTART1, IEND1, &
1143                            position(1), XE(1:size(XE,1)-1),'X','I')
1144                       !-1 above since xe ranges from 0:IMAX3, so size is imax3 + 1.
1145                       !therefore, (1:size(xe,1)) will give 1:imax3 + 1, resulting in a seg error.
1146     
1147                       J = DES_GETINDEXFROMPOS(JSTART1,JEND1,position(2), &
1148                            YN(1:size(YN,1)-1),'Y','J')
1149     
1150                       K = 1
1151                       IF(DO_K) K = DES_GETINDEXFROMPOS(KSTART1, &
1152                            KEND1,position(3),ZT(1:size(ZT,1)-1),'Z','K')
1153     
1154                       IF(.not.IS_ON_myPE_wobnd(I,J,K)) cycle
1155     
1156                       IF(DEAD_CELL_AT(I,J,K)) cycle
1157     
1158                       IJK = FUNIJK(I, J, K)
1159                       IF(.not.FLUID_AT(IJK)) cycle
1160     
1161                       RNP_PIC(M) = RNP_PIC(M) + STAT_WT
1162                       CNP_PIC(M) = CNP_PIC(M) + 1
1163                       PART_MPHASE_BYIC(ICV, M) =  PART_MPHASE_BYIC(ICV, M) &
1164                            + 1
1165     
1166                       REALPART_MPHASE_BYIC(ICV, M) =  REALPART_MPHASE_BYIC(ICV, M) &
1167                            + STAT_WT
1168     
1169                       RAD = DES_D_P0(M)*HALF
1170                       DENS  =  DES_RO_S(M)
1171                       VELOCITY(1:DIMN) = PVEL_TEMP(IPCOUNT,1:DIMN)
1172     
1173                       CALL GEN_AND_ADD_TO_PART_LIST(PART_LIST_BYIC, M, POSITION(1:DIMN), &
1174                            VELOCITY(1:DIMN), RAD, DENS, STAT_WT)
1175                       PART_LIST_BYIC%cell(1) = I
1176                       PART_LIST_BYIC%cell(2) = J
1177                       PART_LIST_BYIC%cell(3) = K
1178                       PART_LIST_BYIC%cell(4) = IJK
1179                    ENDDO
1180                    deallocate(pvel_temp)
1181                    deallocate(randpos)
1182     
1183                 ENDIF
1184     
1185                 CALL global_all_sum(PART_MPHASE_BYIC(ICV, M))
1186                 CALL global_all_sum(REALPART_MPHASE_BYIC(ICV, M))
1187     
1188     
1189                 WRITE(ERR_MSG, '(10x, A, I15, /, 10x, A, ES15.7)') &
1190                      '# of computational particles or parcels acutally seeded = ', PART_MPHASE_BYIC(ICV, M), &
1191                      '# of implied real particles for above parcel count          = ', REALPART_MPHASE_BYIC(ICV, M)
1192     
1193                 CALL FLUSH_ERR_MSG(header = .false., FOOTER = .false.)
1194     
1195              ENDDO MLOOP
1196              ! Now merge this IC specific list to the master list
1197              if(associated(part_list_byic)) then
1198                 CALL merge_part_lists(part_list_byic, orig_part_list)
1199                 nullify(part_list_byic)
1200              endif
1201     
1202     
1203           ENDDO                     !ICV = 1, DIMENSION_IC
1204     
1205     
1206           CALL global_all_sum(RNP_PIC)
1207           CALL global_all_sum(CNP_PIC)
1208     
1209           PART_MPHASE(1:DES_MMAX)     = CNP_PIC(1:DES_MMAX)
1210           REALPART_MPHASE(1:DES_MMAX) = RNP_PIC(1:DES_MMAX)
1211           !particles will be set in gener_part_config routine
1212     
1213     ! Local reporting for the user
1214           WRITE(ERR_MSG, 2015)  COUNT_IC, DES_MMAX
1215           CALL FLUSH_ERR_MSG(HEADER = .false., FOOTER = .false.)
1216     
1217     2015  FORMAT(/, 'Done generating the initial configuration for PIC model', /, &
1218           'Total IC region count with non zero solids = ', I5,2X, /, &
1219           'Total discrete solid phases                = ', I5,2X, /, &
1220           'Global reporting per solid phase follows')
1221     
1222           DO M = 1, DES_MMAX
1223              WRITE(err_msg, 2023) M, PART_MPHASE(M), REALPART_MPHASE(M)
1224              CALL FLUSH_ERR_MSG(HEADER = .false., FOOTER = .false.)
1225           enddo
1226     
1227     
1228     2023  FORMAT(/5x, &
1229                'For solid phase M:                          ',  I5, /5x, &
1230                'Total number of parcels                    = ', I15, /5x, &
1231                'Total number of implied physical particles = ', ES15.7)
1232     
1233           WRITE(ERR_MSG,*) ''
1234           CALL FLUSH_ERR_MSG(HEADER = .false.)
1235     
1236           CALL FINL_ERR_MSG
1237           END SUBROUTINE GENERATE_PARTICLE_CONFIG_MPPIC
1238     
1239     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1240     !                                                                      C
1241     !  SUBROUTINE: DBG_WRITE_PART_VTP_FILE_FROM_LIST                       C
1242     !                                                                      C
1243     !  Purpose: Subroutine to write initial particle lists to vtp files    C
1244     !           for debugging purporse.                                    C
1245     !           This only writes process wise files since constructor does C
1246     !           not exist to gather linked-lists on root processor.        C
1247     !                                                                      C
1248     !  Authors: Rahul Garg                               Date: 21-Mar-2014 C
1249     !                                                                      C
1250     !                                                                      C
1251     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1252           SUBROUTINE DBG_WRITE_PART_VTP_FILE_FROM_LIST(part_fname, writeindomain)
1253           USE run
1254           USE param1
1255           USE discretelement, only : dimn, s_time
1256           USE compar
1257           USE constant
1258           USE cutcell
1259           USE funits
1260           USE indices
1261           USE physprop
1262           USE parallel
1263           USE DES_LINKED_LIST_DATA, only : orig_part_list, particle
1264           USE des_linked_list_funcs
1265           use geometry, only: NO_K
1266     ! Use the error manager for posting error messages.
1267     !---------------------------------------------------------------------//
1268           use error_manager
1269     
1270           Implicit none
1271           CHARACTER (LEN=*), intent(in) :: part_fname
1272           logical , intent(in) :: writeindomain
1273           !facet id and particle id
1274           Integer ::  vtp_unit , k, nparts, pvd_unit, ipe
1275           CHARACTER(LEN=100) :: vtp_fname, pvd_fname
1276     
1277           type(particle), pointer :: part => null()
1278           ! dummy values to maintain format for dimn=2
1279           REAL POS_Z, VEL_W
1280     
1281     ! formatted solids time
1282           CHARACTER(LEN=12) :: S_TIME_CHAR = ''
1283     ! Initialize the error manager.
1284           CALL INIT_ERR_MSG("WRITE_PARTICLE_VTP_FILE")
1285           vtp_unit = 1002
1286     
1287           POS_Z = 0.0
1288           vel_w = 0.0
1289     
1290           nparts = 0
1291           part => orig_part_list
1292           DO WHILE (ASSOCIATED(part))
1293              if(part%INDOMAIN.and.writeindomain) then
1294                 nparts = nparts + 1
1295              elseif(.not.part%INDOMAIN.and..not.writeindomain) then
1296                 nparts = nparts + 1
1297              endif
1298              part => part%next
1299           ENDDO
1300     
1301           if(mype.eq.pe_io) then
1302              write(pvd_fname,'(A,"_", A, ".pvd")') &
1303              trim(run_name),trim(part_fname)
1304              pvd_unit = 1002
1305              open(pvd_unit, file = pvd_fname, form='formatted')
1306     
1307     
1308              WRITE(PVD_UNIT,"(A)")'<?xml version="1.0"?>'
1309              WRITE(PVD_UNIT,"(A,A)")'<VTKFile type="Collection" ',&
1310              'version="0.1" byte_order="LittleEndian">'
1311              WRITE(PVD_UNIT,"(3X,A)")'<Collection>'
1312     
1313              S_TIME_CHAR = '0.0'
1314              do ipe = 0, numpes-1
1315                 write(vtp_fname,'(A,"_", A, "_",I5.5,".vtp")') &
1316                 trim(run_name),trim(part_fname), ipe
1317                 WRITE(PVD_UNIT,"(6X,A,A,A,A,A,A,A)")&
1318                 '<DataSet timestep="',TRIM(S_TIME_CHAR),'" ',& ! simulation time
1319                 'group="" part="0" ',& ! necessary file data
1320                 'file="',TRIM(VTP_FNAME),'"/>' ! file name of vtp
1321     
1322              enddo
1323     
1324     
1325     ! Write the closing tags
1326              WRITE(PVD_UNIT,"(3X,A)")'</Collection>'
1327              WRITE(PVD_UNIT,"(A)")'</VTKFile>'
1328     
1329              CLOSE(PVD_UNIT)
1330     
1331           endif
1332     
1333           write(vtp_fname,'(A,"_", A, "_",I5.5,".vtp")') &
1334           trim(run_name),trim(part_fname), mype
1335     
1336           open(vtp_unit, file = vtp_fname, form='formatted')
1337           write(vtp_unit,"(a)") '<?xml version="1.0"?>'
1338     ! Write the S_TIME as a comment for reference
1339           write(vtp_unit,"(a,es24.16,a)") '<!-- Time =',s_time,'s -->'
1340     
1341     ! Write the necessary header information for a PolyData file type
1342           write(vtp_unit,"(a,a)") '<VTKFile type="PolyData"',&
1343           ' version="0.1" byte_order="LittleEndian">'
1344           write(vtp_unit,"(3x,a)") '<PolyData>'
1345     
1346     ! Write Piece tag and identify the number of particles in the system.
1347           write(vtp_unit,"(6x,a,i10.10,a,a)")&
1348           '<Piece NumberOfPoints="',nparts,&
1349           '" NumberOfVerts="0" ',&
1350           'NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0">'
1351           write(vtp_unit,"(9x,a)")&
1352           '<PointData Scalars="Diameter" Vectors="Velocity">'
1353     
1354     ! Write the diameter data.
1355           write(vtp_unit,"(12x,a)")&
1356           '<DataArray type="Float32" Name="Diameter" format="ascii">'
1357           part => orig_part_list
1358           DO WHILE (ASSOCIATED(part))
1359              if(part%INDOMAIN.and.writeindomain) then
1360                 write (vtp_unit,"(15x,es13.6)") (real(2.d0*part%rad))
1361              elseif(.not.part%INDOMAIN.and..not.writeindomain) then
1362                 write (vtp_unit,"(15x,es13.6)") (real(2.d0*part%rad))
1363              endif
1364              part => part%next
1365           ENDDO
1366     
1367           write(vtp_unit,"(12x,a)") '</DataArray>'
1368     
1369           ! Write velocity data. Force to three dimensions. So for 2D runs, a
1370     ! dummy value of zero is supplied as the 3rd point.
1371           write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
1372           'Name="Velocity" NumberOfComponents="3" format="ascii">'
1373           if(NO_K) then
1374              part => orig_part_list
1375              DO WHILE (ASSOCIATED(part))
1376                 if(part%INDOMAIN.and.writeindomain) then
1377                    write (vtp_unit,"(15x,3(es13.6,3x))")&
1378                    (real(part%velocity(k)),k=1,dimn),vel_w
1379                 elseif(.not.part%INDOMAIN.and..not.writeindomain) then
1380                    write (vtp_unit,"(15x,3(es13.6,3x))")&
1381                    (real(part%velocity(k)),k=1,dimn),vel_w
1382                 endif
1383                 part => part%next
1384              ENDDO
1385     
1386           else                      ! 3d
1387              part => orig_part_list
1388              DO WHILE (ASSOCIATED(part))
1389                 if(part%INDOMAIN.and.writeindomain) then
1390                    write (vtp_unit,"(15x,3(es13.6,3x))")&
1391                         (real(part%velocity(k)),k=1,dimn)
1392                 elseif(.not.part%INDOMAIN.and..not.writeindomain) then
1393                    write (vtp_unit,"(15x,3(es13.6,3x))")&
1394                         (real(part%velocity(k)),k=1,dimn)
1395                 endif
1396                 part => part%next
1397              ENDDO
1398           endif
1399           write(vtp_unit,"(12x,a,/9x,a)") '</DataArray>','</PointData>'
1400     
1401     ! Skip CellData tag, no data. (vtp format style)
1402           write(vtp_unit,"(9x,a)") '<CellData></CellData>'
1403     
1404     ! Write position data. Point data must be supplied in 3 dimensions. So
1405     ! for 2D runs, a dummy value of zero is supplied as the 3rd point.
1406           write(vtp_unit,"(9x,a)") '<Points>'
1407           write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
1408           'Name="Position" NumberOfComponents="3" format="ascii">'
1409           if(NO_K) then
1410     
1411              part => orig_part_list
1412              DO WHILE (ASSOCIATED(part))
1413                 if(part%INDOMAIN.and.writeindomain) then
1414                    write (vtp_unit,"(15x,3(es13.6,3x))")&
1415                    (real(part%position(k)),k=1,dimn),pos_z
1416                 elseif(.not.part%INDOMAIN.and..not.writeindomain) then
1417                    write (vtp_unit,"(15x,3(es13.6,3x))")&
1418                         (real(part%position(k)),k=1,dimn),pos_z
1419                 endif
1420                 part => part%next
1421              ENDDO
1422           else
1423     
1424              part => orig_part_list
1425              DO WHILE (ASSOCIATED(part))
1426                 if(part%INDOMAIN.and.writeindomain) then
1427                    write (vtp_unit,"(15x,3(es13.6,3x))")&
1428                    (real(part%position(k)),k=1,dimn)
1429                 elseif(.not.part%INDOMAIN.and..not.writeindomain) then
1430                    write (vtp_unit,"(15x,3(es13.6,3x))")&
1431                         (real(part%position(k)),k=1,dimn)
1432                 endif
1433     
1434                 part => part%next
1435              ENDDO
1436           endif
1437           write(vtp_unit,"(12x,a,/9x,a)")'</DataArray>','</Points>'
1438     
1439     ! Write tags for data not included (vtp format style)
1440           write(vtp_unit,"(9x,a,/9x,a,/9x,a,/9x,a)")'<Verts></Verts>',&
1441           '<Lines></Lines>','<Strips></Strips>','<Polys></Polys>'
1442           write(vtp_unit,"(6x,a,/3x,a,/a)")&
1443           '</Piece>','</PolyData>','</VTKFile>'
1444     
1445           close(vtp_unit, status = 'keep')
1446     
1447     
1448           write(ERR_MSG, 1000) part_fname
1449     
1450      1000 FORMAT( 'For debugging purposes, processor wise particle configurations ', &
1451           'written to files containing ', /, A, /)
1452     
1453           CALL FLUSH_ERR_MSG(header=.false., footer = .false.)
1454     
1455           CALL FINL_ERR_MSG
1456           END SUBROUTINE  DBG_WRITE_PART_VTP_FILE_FROM_LIST
1457     
1458     
1459