File: /nfs/home/0/users/jenkins/mfix.git/model/des/generate_particle_config.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15 SUBROUTINE GENERATE_PARTICLE_CONFIG
16
17
18
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
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
45
46 = .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)
121
122
123
124 = 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
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
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
225
226
227
228
229
230
231
232
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
252
253
254
255
256
257
258
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
267
268
269 IF(part%indomain) THEN
270
271
272
273
274
275 IF(FLUID_AT(IJK).and.(.not.CUT_CELL_AT(IJK))) THEN
276
277
278
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
293
294
295
296
297
298
299
300
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
321
322
323
324
325
326
327
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
336
337
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
353
354
355
356
357
358
359
360
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
379 CALL INIT_ERR_MSG("BIN_PARTICLES_TO_CELL")
380
381
382 => 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
387
388
389 %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.
400 ENDIF
401
402 part => part%next
403 ENDDO
404
405 CALL FINL_ERR_MSG
406 END SUBROUTINE BIN_PARTICLES_TO_CELL
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423 SUBROUTINE GENERATE_PARTICLE_CONFIG_DEM
424
425
426
427
428 USE discretelement, only: GENER_PART_CONFIG
429
430 USE discretelement, only: DES_RADIUS, RO_Sol
431
432 USE discretelement, only: des_pos_new, des_pos_old
433
434 USE discretelement, only: des_vel_new, des_vel_old
435
436 USE discretelement, only: DIMN
437
438 USE discretelement, only: XE, YN, ZT
439
440 USE discretelement, only: PIP
441
442 USE discretelement, only: DES_MMAX
443
444 USE discretelement, only: DES_D_p0, DES_RO_s
445
446 USE discretelement, only: VOL_IC_REGION
447
448 USE discretelement, only: PART_MPHASE_BYIC
449
450 USE discretelement, only: REALPART_MPHASE_BYIC
451
452
453 USE discretelement, only: PARTICLES
454
455 USE discretelement, only: PEA
456
457
458 USE discretelement, only:DES_GETINDEXFROMPOS
459 USE constant, only: PI
460
461 USE ic, only: IC_DEFINED
462
463 USE ic, only: IC_ROP_s
464
465 USE ic, only: IC_EP_S
466
467 USE ic, only: IC_EP_G
468
469 USE ic, only: IC_X_w, IC_X_e, IC_Y_s, IC_Y_n, IC_Z_b, IC_Z_t
470
471 USE ic, only: IC_U_s, IC_V_s, IC_W_s, IC_Theta_M
472
473 Use ic, only: IC_DES_FIT_TO_REGION
474
475 USE param1, only: UNDEFINED, UNDEFINED_I, ZERO, ONE, Half
476
477 USE param1, only: SMALL_NUMBER, LARGE_NUMBER
478
479 USE param, only: DIMENSION_IC
480
481 USe compar, only: istart1, iend1, jstart1, jend1, kstart1, kend1
482
483 USE randomno
484 USE mpi_utility
485
486 use error_manager
487
488
489
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
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
508 double precision :: max_ic_pt(3)
509
510 double precision :: refit_fac(3), volijk
511
512 double precision :: PART_CEN_MIN(DIMN), PART_CEN_MAX(DIMN)
513
514
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
526 = 0
527
528 PCOUNT_BYIC_BYPHASE(:,:) = 0
529
530
531 = 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
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
595 (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
609
610
611
612
613 (:) = 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
647
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
662
663
664
665 (ICV,M) = PCOUNT_BYIC_BYPHASE(ICV,M) + 1
666
667
668
669
670
671
672 = 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
712
713 = 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
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
733
734 => part_list_byic
735 DO WHILE (ASSOCIATED(part))
736
737 => part
738 part => part%next
739
740
741
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
750
751
752 = 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
772
773
774 %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
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
813
814
815
816
817
818
819 SUBROUTINE GENERATE_PARTICLE_CONFIG_MPPIC
820
821
822
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
832 USE discretelement, only: DES_MMAX
833
834
835 USE discretelement, only: DES_D_p0, DES_RO_s
836
837 USE discretelement, only: PART_MPHASE_BYIC
838
839 USE discretelement, only: REALPART_MPHASE_BYIC
840
841
842 USE discretelement, only: PART_MPHASE
843
844 USE discretelement, only: REALPART_MPHASE
845
846 USE mfix_pic, only : mppic
847
848
849 USE mfix_pic, only: rnp_pic
850
851 USE mfix_pic, only: cnp_pic
852
853 USE mfix_pic, only: cnp_array
854
855
856 USE ic, only: IC_DEFINED
857
858 USE ic, only: IC_ROP_s
859
860 USE ic, only: IC_EP_G
861
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
868 USE ic, only: IC_U_s, IC_V_s, IC_W_s, IC_Theta_M
869
870
871 USE cutcell, only: cut_cell_at
872
873 USE param1, only: UNDEFINED, UNDEFINED_I, ZERO, ONE
874
875
876 USE ic, only: IC_PIC_CONST_NPC, IC_PIC_CONST_STATWT
877
878
879 USE param, only: DIMENSION_IC
880
881
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
893
894 DOUBLE PRECISION :: EP_SM
895
896 LOGICAL :: CONST_NPC, CONST_STATWT
897
898 INTEGER :: I, J, K, IJK, ICV, M, COUNT_IC, IDIM
899
900
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
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
972
973
974
975
976 (:) = 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
996 (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
1024
1025
1026
1027 = 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
1111
1112
1113
1114
1115
1116 = 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
1145
1146
1147 = 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
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
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
1212
1213
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
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
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
1267
1268 use error_manager
1269
1270 Implicit none
1271 CHARACTER (LEN=*), intent(in) :: part_fname
1272 logical , intent(in) :: writeindomain
1273
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
1279 REAL POS_Z, VEL_W
1280
1281
1282 CHARACTER(LEN=12) :: S_TIME_CHAR = ''
1283
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),'" ',&
1319 'group="" part="0" ',&
1320 'file="',TRIM(VTP_FNAME),'"/>'
1321
1322 enddo
1323
1324
1325
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
1339 write(vtp_unit,"(a,es24.16,a)") '<!-- Time =',s_time,'s -->'
1340
1341
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
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
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
1370
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
1387 => 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
1402 write(vtp_unit,"(9x,a)") '<CellData></CellData>'
1403
1404
1405
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
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