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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: PIC_APPLY_WALLBC_STL                                    !
4     !  Author: R. Garg                                                     !
5     !                                                                      !
6     !  Purpose: Detect collisions between PIC particles and STLs.          !
7     !                                                                      !
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9           SUBROUTINE PIC_APPLY_WALLBC_STL
10     
11     ! Global Variables:
12     !---------------------------------------------------------------------//
13     ! Particle postions: Current/Previous
14           use discretelement, only: DES_POS_NEW, DES_POS_OLD
15     ! Paricle velocities
16           use discretelement, only: DES_VEL_NEW
17     ! Particle radi
18           use discretelement, only: DES_RADIUS
19     ! Flag indicating state of particle
20           use discretelement, only: PEA
21     ! Max number of particles on this process
22           use discretelement, only: MAX_PIP
23     ! The number of neighbor facets for each DES grid cell
24           use discretelement, only: cellneighbor_facet_num
25     ! Facet informatin binned to the DES grid
26           use discretelement, only: cellneighbor_facet
27     ! Map from particle to DES grid cell.
28           use discretelement, only: DG_PIJK
29     ! Flag indicating that the index cell contains no STLs
30           use stl, only: NO_NEIGHBORING_FACET_DES
31     ! STL Vertices
32           use stl, only: VERTEX
33     ! STL Facet normals
34           use stl, only: NORM_FACE
35     ! The number of mass inflow/outflow BCs
36           use pic_bc, only: PIC_BCMI, PIC_BCMO
37     
38     ! Global Parameters:
39     !---------------------------------------------------------------------//
40           use param1, only: ZERO, ONE, UNDEFINED
41     ! DES array dimensionality (3)
42           use discretelement, only: DIMN
43     
44     ! Module Procedures:
45     !---------------------------------------------------------------------//
46           use des_stl_functions, only: intersectLnPlane
47           use des_stl_functions, only: checkPTonTriangle
48     
49           use error_manager
50     
51           IMPLICIT NONE
52     
53     ! Local Variables:
54     !---------------------------------------------------------------------//
55     ! Loop counters.
56           INTEGER NF, LL
57     ! Loop counter for SLTs in curring DES grid cell
58           INTEGER :: CELL_COUNT
59     ! reference point and direction of the line
60           DOUBLE PRECISION, DIMENSION(DIMN) :: REF_LINE,  DIR_LINE
61     ! reference point and normal of the plane
62           DOUBLE PRECISION, DIMENSION(DIMN) :: REF_PLANE, NORM_PLANE
63     ! line is parameterized as p = p_ref + t * dir_line, t is line_param
64           DOUBLE PRECISION :: LINE_T
65     ! logical for determining if a point is on the triangle or not
66           LOGICAL :: ONTRIANGLE
67     ! The line and plane intersection point
68           DOUBLE PRECISION, DIMENSION(DIMN) :: POINT_ONPLANE
69     ! Old and new position of the parcel
70     
71     ! distance parcel travelled out of domain along the facet normal
72           DOUBLE PRECISION :: dist_from_facet
73           DOUBLE PRECISION :: veldotnorm
74     
75           CALL INIT_ERR_MSG("PIC_APPLY_WALLBC_STL")
76     
77           DO LL = 1, MAX_PIP
78     
79     ! Skip non-existent particles
80              IF(.NOT.PEA(LL,1)) CYCLE
81     
82     ! If no neighboring facet in the surrounding 27 cells, then exit
83              IF (NO_NEIGHBORING_FACET_DES(DG_PIJK(LL))) CYCLE
84     
85     ! Loop through the STLs associated with the DES grid cell.
86              CELLLOOP: DO CELL_COUNT=1, CELLNEIGHBOR_FACET_NUM(DG_PIJK(LL))
87     
88     ! Get the index of the neighbor facet
89                 NF = CELLNEIGHBOR_FACET(DG_PIJK(LL))%P(CELL_COUNT)
90     
91     ! Set to -UNDEFINED, because non zero values imply intersection.
92     ! with the plane.
93                 LINE_T  = -UNDEFINED
94     
95                 ONTRIANGLE = .FALSE.
96     
97     ! Parametrize a line as p = p_0 + t (p_1 - p_0) and intersect with the
98     ! triangular plane. If 0<t<1, then this ray connect old and new
99     ! positions intersects with the facet.
100                 REF_LINE(1:DIMN) = DES_POS_OLD(:,LL)
101                 DIR_LINE(1:DIMN) = DES_POS_NEW(:,LL) - DES_POS_OLD(:,LL)
102     
103                 NORM_PLANE = NORM_FACE(:,NF)
104                 REF_PLANE  = VERTEX(1,:,NF)
105     
106                 VELDOTNORM = DOT_PRODUCT(NORM_PLANE, DES_VEL_NEW(:,LL))
107     
108     ! If the normal velocity of parcel is in the same direction as facet
109     ! normal, then the parcel could not have crossed this facet.
110                 IF(VELDOTNORM.GT.0.d0) CYCLE
111     
112                 CALL INTERSECTLNPLANE(REF_LINE, DIR_LINE, REF_PLANE,       &
113                      NORM_PLANE, LINE_T)
114     
115     ! line_t = one implies that the particle is sitting on the stl face.
116                 IF(LINE_T.GE.ZERO.AND.LINE_T.LT.ONE) THEN
117     
118                    POINT_ONPLANE = REF_LINE + LINE_T*DIR_LINE
119     
120     ! Check through barycentric coordinates to determine if the orthogonal
121     ! projection lies on the facet or not. If it does, then this point is
122     ! deemed to be on the non-fluid side of the facet.
123                    CALL CHECKPTONTRIANGLE(POINT_ONPLANE(:), VERTEX(:,:,NF),&
124                       ONTRIANGLE)
125     
126                    IF(ONTRIANGLE) THEN
127     
128                       DIST_FROM_FACET = ABS(DOT_PRODUCT(DES_POS_NEW(:,LL) -&
129                          POINT_ONPLANE(:), NORM_PLANE(:)))
130     
131                       DES_POS_NEW(:,LL) = POINT_ONPLANE +                  &
132                          DIST_FROM_FACET*NORM_PLANE
133     
134     ! In the reflect routine, it is assumed that the normal points into the
135     ! system which is consitent with the definition of normal for facets.
136                       CALL PIC_REFLECT_PART(LL, NORM_PLANE(:))
137     
138                       EXIT CELLLOOP
139                    ENDIF
140                 ENDIF
141     
142              END DO CELLLOOP
143     
144           END DO
145     
146     ! Seed new parcels entering the system.
147           IF(PIC_BCMI > 0) CALL PIC_MI_BC
148           IF(PIC_BCMO > 0) CALL PIC_MO_BC
149     
150           CALL FINL_ERR_MSG
151     
152           RETURN
153           END SUBROUTINE PIC_APPLY_WALLBC_STL
154     
155     
156     
157     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
158     !                                                                      !
159     !  Subroutine: Mass_OUTFLOW_PIC                                        !
160     !  Author: R. Garg                                   Date: 23-Jun-14   !
161     !                                                                      !
162     !  Purpose:  Routine to delete out of domain parcels for PIC           !
163     !  implementation                                                      !
164     !                                                                      !
165     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
166           SUBROUTINE PIC_MO_BC
167     
168           use discretelement
169           use pic_bc
170           use bc
171           USE error_manager
172           USE mpi_utility
173     
174     
175           implicit none
176     
177           INTEGER :: IJK
178           INTEGER :: LC, LP, NP
179           INTEGER :: BCV, BCV_I
180     
181           DOUBLE PRECISION :: DIST
182     
183           INTEGER :: PIP_DEL_COUNT
184           INTEGER, dimension(0:numpes-1) :: del_count_all
185           INTEGER :: IPROC
186     
187     
188           CALL INIT_ERR_MSG("PIC_MO_BC")
189     
190           PIP_DEL_COUNT = 0
191           DO BCV_I = 1, PIC_BCMO
192     
193              BCV = PIC_BCMO_MAP(BCV_I)
194     
195              SELECT CASE (BC_PLANE(BCV))
196              END SELECT
197     
198              DO LC=PIC_BCMO_IJKSTART(BCV_I), PIC_BCMO_IJKEND(BCV_I)
199                 IJK = PIC_BCMO_IJK(LC)
200                 DO LP= 1,PINC(IJK)
201     
202                    NP = PIC(IJK)%p(LP)
203     
204                    IF(PEA(NP,4)) CYCLE
205     
206                    SELECT CASE (BC_PLANE(BCV))
207                    CASE('S'); DIST = YN(BC_J_s(BCV)-1) - DES_POS_NEW(2,NP)
208                    CASE('N'); DIST = DES_POS_NEW(2,NP) - YN(BC_J_s(BCV))
209                    CASE('W'); DIST = XE(BC_I_w(BCV)-1) - DES_POS_NEW(1,NP)
210                    CASE('E'); DIST = DES_POS_NEW(1,NP) - XE(BC_I_w(BCV))
211                    CASE('B'); DIST = ZT(BC_K_b(BCV)-1) - DES_POS_NEW(3,NP)
212                    CASE('T'); DIST = DES_POS_NEW(3,NP) - ZT(BC_K_b(BCV))
213                    END SELECT
214     
215                    IF(DIST .lt. zero ) THEN
216     
217     ! The parcel has left this bc_plane
218     
219                       CALL DELETE_PARCEL(NP)
220                       PIP_DEL_COUNT = PIP_DEL_COUNT+1
221                    ENDIF
222     
223                 ENDDO
224              ENDDO
225           ENDDO
226     
227     
228           IF(PIC_REPORT_DELETION_STATS) then
229     
230              del_count_all(:) = 0
231              del_count_all(mype) = pip_del_count
232              call global_all_sum(del_count_all(0:numpes-1))
233     
234              IF(SUM(del_COUNT_ALL(:)).GT.0) THEN
235                 WRITE(err_msg,'(/,2x,A,2x,i10)')&
236                 'TOTAL NUMBER OF PARCELS DELETED GLOBALLY = ', SUM(DEL_COUNT_ALL(:))
237     
238                 call flush_err_msg(header = .false., footer = .false.)
239     
240                 DO IPROC = 0, NUMPES-1
241                    IF(DMP_LOG) WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
242                         'PARCELS ADDED ON PROC:', IPROC,' EQUAL TO', pip_del_count
243                 ENDDO
244              ENDIF
245     
246           ENDIF
247     
248     
249           CALL FINL_ERR_MSG
250           RETURN
251           END SUBROUTINE PIC_MO_BC
252     
253     
254     
255     
256     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
257     !                                                                      !
258     !  Subroutine: DELETE_PARCEL                                           !
259     !  Author: R. Garg                                    Date: 23-Jun-14  !
260     !                                                                      !
261     !  Purpose:  Routine to delete parcel                                  !
262     !                                                                      !
263     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
264           SUBROUTINE DELETE_PARCEL(NP)
265     
266           USE compar
267           USE constant
268           USE des_bc
269           USE discretelement
270           USE funits
271           USE geometry
272           USE indices
273           USE param1
274           USE physprop
275           USE mfix_pic
276           USE functions
277     
278           IMPLICIT NONE
279     
280     
281           INTEGER, INTENT(IN) :: NP
282     
283           PEA(NP,:) = .FALSE.
284     
285           DES_POS_OLD(:,NP) = ZERO
286           DES_POS_NEW(:,NP) = ZERO
287           DES_VEL_OLD(:,NP) = ZERO
288           DES_VEL_NEW(:,NP) = ZERO
289           OMEGA_OLD(:,NP) = ZERO
290           OMEGA_NEW(:,NP) = ZERO
291           DES_RADIUS(NP) = ZERO
292           PMASS(NP) = ZERO
293           PVOL(NP) = ZERO
294           RO_Sol(NP) = ZERO
295           OMOI(NP) = ZERO
296     
297           DES_STAT_WT(NP) = ZERO
298     
299           FC(:,NP) = ZERO
300     
301           PIP = PIP - 1
302     
303           RETURN
304           END SUBROUTINE DELETE_PARCEL
305     
306     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
307     !                                                                      C
308     !  Subroutine: MPPIC_MI_BC                                             C
309     !  Purpose:                                                            C
310     !                                                                      C
311     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
312     
313           SUBROUTINE  PIC_MI_BC
314     
315           USE run
316           USE discretelement
317           USE geometry
318           USE constant
319           USE cutcell
320           USE mfix_pic
321           USE pic_bc
322           USE bc
323           USE mpi_utility
324           USE randomno
325           use error_manager
326           USE functions
327     
328           IMPLICIT NONE
329     !-----------------------------------------------
330     ! Local variables
331     !-----------------------------------------------
332           INTEGER :: WDIR, IDIM, IPCOUNT, LAST_EMPTY_SPOT, NEW_SPOT
333           INTEGER :: BCV, BCV_I, L, LC, PIP_ADD_COUNT, IPROC
334           INTEGER ::  IFLUID, JFLUID, KFLUID, IJK, M
335           DOUBLE PRECISION :: CORD_START(3), DOML(3),WALL_NORM(3)
336           DOUBLE PRECISION :: AREA_INFLOW, VEL_INFLOW(DIMN), STAT_WT
337     
338           LOGICAL :: DELETE_PART
339           DOUBLE PRECISION :: EPS_INFLOW(DES_MMAX)
340           DOUBLE PRECISION REAL_PARTS(DES_MMAX), COMP_PARTS(DES_MMAX), VEL_NORM_MAG, VOL_INFLOW, VOL_IJK
341     
342     
343     ! Temp logical variables for checking constant npc and statwt specification
344           LOGICAL :: CONST_NPC, CONST_STATWT
345     
346           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: RANDPOS
347           INTEGER :: CNP_CELL_COUNT
348           DOUBLE PRECISION :: RNP_CELL_COUNT
349           integer :: lglobal_id
350           integer, dimension(0:numpes-1) :: add_count_all
351     
352           LOGICAL, parameter :: setDBG = .true.
353           LOGICAL :: dFlag
354     
355           type :: ty_spotlist
356              integer spot
357              type(ty_spotlist),pointer :: next => NULL()
358           end type ty_spotlist
359     
360           type(ty_spotlist),pointer :: &
361                cur_spotlist => NULL(), &
362                prev_spotlist => NULL(), &
363                temp_spotlist => NULL()
364     
365     !-----------------------------------------------
366     
367           CALL INIT_ERR_MSG("PIC_MI_BC")
368     
369           dFlag = (DMP_LOG .AND. setDBG)
370           PIP_ADD_COUNT = 0
371           LAST_EMPTY_SPOT = 0
372     
373           ALLOCATE(cur_spotlist); cur_spotlist%spot = -1
374     
375           DO BCV_I = 1, PIC_BCMI
376              BCV = PIC_BCMI_MAP(BCV_I)
377     
378              WALL_NORM(1:3) =  PIC_BCMI_NORMDIR(BCV_I,1:3)
379     
380              !Find the direction of the normal for this wall
381              WDIR = 0
382              DO IDIM = 1, DIMN
383                 WDIR = WDIR + ABS(WALL_NORM(IDIM))*IDIM
384              end DO
385     
386              DO LC=PIC_BCMI_IJKSTART(BCV_I), PIC_BCMI_IJKEND(BCV_I)
387                 IJK = PIC_BCMI_IJK(LC)
388     
389                 IF(.NOT.FLUID_AT(IJK)) CYCLE
390     
391                 IFLUID = I_OF(IJK)
392                 JFLUID = J_OF(IJK)
393                 KFLUID = K_OF(IJK)
394     
395                 CORD_START(1) = XE(IFLUID) - PIC_BCMI_OFFSET (BCV_I,1)*DX(IFLUID)
396     
397                 CORD_START(2) = YN(JFLUID) - PIC_BCMI_OFFSET (BCV_I,2)*DY(JFLUID)
398     
399     
400                 CORD_START(3) = merge(zero, ZT(KFLUID) - PIC_BCMI_OFFSET (BCV_I,3)*DZ(KFLUID), no_k)
401     
402                 DOML(1) = DX(IFLUID)
403                 DOML(2) = DY(JFLUID)
404                 DOML(3) = MERGE(DZ(1), DZ(KFLUID), NO_K)
405     
406                 AREA_INFLOW = DOML(1)*DOML(2)*DOML(3)/DOML(WDIR)
407     
408                 VOL_IJK = DOML(1)*DOML(2)*DOML(3)
409     
410                 DOML(WDIR) = ZERO
411                 !set this to zero as the particles will
412                 !be seeded only on the BC plane
413     
414                 DO M = 1, DES_MMAX
415     
416                    IF(SOLIDS_MODEL(M) /= 'PIC') CYCLE
417                    EPS_INFLOW(M) = BC_ROP_S(BCV, M)/DES_RO_S(M)
418                    VEL_INFLOW(1) = BC_U_S(BCV, M)
419                    VEL_INFLOW(2) = BC_V_S(BCV, M)
420                    VEL_INFLOW(3) = BC_W_S(BCV, M)
421     
422                    VEL_NORM_MAG = ABS(DOT_PRODUCT(VEL_INFLOW(1:DIMN), WALL_NORM(1:DIMN)))
423                    VOL_INFLOW   = AREA_INFLOW*VEL_NORM_MAG*DTSOLID
424     
425                    REAL_PARTS(M) = 6.d0*EPS_INFLOW(M)*VOL_INFLOW/(PI*(DES_D_p0(M)**3.d0))
426                    COMP_PARTS(M) = zero
427     
428                    CONST_NPC    = (BC_PIC_MI_CONST_NPC   (BCV, M) .ne. 0)
429                    CONST_STATWT = (BC_PIC_MI_CONST_STATWT(BCV, M) .ne. ZERO)
430                    IF(CONST_NPC) THEN
431                       IF(EPS_INFLOW(M).GT.ZERO) &
432                       COMP_PARTS(M) = REAL(BC_PIC_MI_CONST_NPC(BCV, M))* &
433                       VOL_INFLOW/VOL_IJK
434                    ELSEIF(CONST_STATWT) THEN
435                       COMP_PARTS(M) = REAL_PARTS(M)/ &
436                       BC_PIC_MI_CONST_STATWT(BCV, M)
437                    ENDIF
438     
439                    pic_bcmi_rnp(LC,M) = pic_bcmi_rnp(LC,M) + REAL_PARTS(M)
440     
441                    pic_bcmi_cnp(LC,M) = pic_bcmi_cnp(LC,M) + COMP_PARTS(M)
442     
443     
444                    IF(pic_bcmi_cnp(LC,M).GE.1.d0) THEN
445                       CNP_CELL_COUNT = INT(pic_bcmi_cnp(LC,M))
446                       pic_bcmi_cnp(LC,M) = pic_bcmi_cnp(LC,M) - CNP_CELL_COUNT
447     
448                       RNP_CELL_COUNT = pic_bcmi_rnp(LC,M)
449                       pic_bcmi_rnp(LC,M)  = zero
450     
451                       !set pic_bcmi_rnp to zero to reflect that all real particles have been seeded
452                       STAT_WT = RNP_CELL_COUNT/REAL(CNP_CELL_COUNT)
453                       ALLOCATE(RANDPOS(CNP_CELL_COUNT*DIMN))
454                       CALL UNI_RNO(RANDPOS(:))
455     
456                       DO IPCOUNT = 1, CNP_CELL_COUNT
457     
458     
459                          CALL PIC_FIND_EMPTY_SPOT(LAST_EMPTY_SPOT, NEW_SPOT)
460     
461                          DES_POS_OLD(1:DIMN, NEW_SPOT) =  CORD_START(1:DIMN) &
462                               & + RANDPOS((IPCOUNT-1)*DIMN+1: &
463                               & (IPCOUNT-1)*DIMN+DIMN)*DOML(1:DIMN)
464                          DES_POS_NEW(:, NEW_SPOT) = DES_POS_OLD(:,NEW_SPOT)
465                          DES_VEL_OLD(1:DIMN,NEW_SPOT) = VEL_INFLOW(1:DIMN)
466     
467                          DES_VEL_NEW(:, NEW_SPOT) = DES_VEL_OLD(:, NEW_SPOT)
468     
469                          DES_RADIUS(NEW_SPOT) = DES_D_p0(M)*HALF
470     
471                          RO_Sol(NEW_SPOT) =  DES_RO_S(M)
472     
473                          DES_STAT_WT(NEW_SPOT) = STAT_WT
474     
475                          PIJK(NEW_SPOT, 1) = IFLUID
476                          PIJK(NEW_SPOT, 2) = JFLUID
477                          PIJK(NEW_SPOT, 3) = KFLUID
478                          PIJK(NEW_SPOT, 4) = IJK
479                          PIJK(NEW_SPOT, 5) = M
480     
481                          PVOL(NEW_SPOT) = (4.0d0/3.0d0)*Pi*DES_RADIUS(NEW_SPOT)**3
482                          PMASS(NEW_SPOT) = PVOL(NEW_SPOT)*RO_SOL(NEW_SPOT)
483     
484                          FC(:, NEW_SPOT) = zero
485                          DELETE_PART = .false.
486                          IF(PIC_BCMI_INCL_CUTCELL(BCV_I)) &
487                               CALL CHECK_IF_PARCEL_OVELAPS_STL &
488                               (des_pos_new(1:dimn, NEW_SPOT), &
489                               DELETE_PART)
490     
491                          IF(.NOT.DELETE_PART) THEN
492     
493                             PIP = PIP+1
494                             PIP_ADD_COUNT = PIP_ADD_COUNT + 1
495                             PEA(NEW_SPOT, 1) = .true.
496                             PEA(NEW_SPOT, 2:4) = .false.
497                             ! add to the list
498                             ALLOCATE(temp_spotlist)
499                             temp_spotlist%spot = new_spot
500                             temp_spotlist%next => cur_spotlist
501                             cur_spotlist => temp_spotlist
502                             nullify(temp_spotlist)
503     
504                          ELSE
505                             PEA(NEW_SPOT, 1) = .false.
506                             LAST_EMPTY_SPOT = NEW_SPOT - 1
507                          ENDIF
508     
509     
510                          !WRITE(*,'(A,2(2x,i5), 2x, A,2x, 3(2x,i2),2x, A, 3(2x,g17.8))') &
511                          !   'NEW PART AT ', NEW_SPOT, MAX_PIP, 'I, J, K = ', IFLUID, JFLUID, KFLUID, 'POS =', DES_POS_NEW(:,NEW_SPOT)
512                          !IF(DMP_LOG) WRITE(UNIT_LOG,'(A,2x,i5, 2x, A,2x, 3(2x,i2),2x, A, 3(2x,g17.8))') &
513                          !    'NEW PART AT ', NEW_SPOT, 'I, J, K = ', IFLUID, JFLUID, KFLUID, 'POS =', DES_POS_NEW(:,NEW_SPOT)
514     
515                          !WRITE(*,*) 'WDIR, DOML = ', WDIR, DOML(:)
516                       END DO
517                       DEALLOCATE(RANDPOS)
518     
519                    end IF
520                 end DO
521              end DO
522           end DO
523     
524     !Now assign global id to new particles added
525           add_count_all(:) = 0
526           add_count_all(mype) = pip_add_count
527           call global_all_sum(add_count_all(0:numpes-1))
528           lglobal_id = imax_global_id + sum(add_count_all(0:mype-1))
529     
530           do l = 1,pip_add_count
531              lglobal_id = lglobal_id + 1
532              iglobal_id(cur_spotlist%spot)= lglobal_id
533              prev_spotlist=> cur_spotlist
534              cur_spotlist => cur_spotlist%next
535              deallocate(prev_spotlist)
536           end do
537           deallocate(cur_spotlist)
538           imax_global_id = imax_global_id + sum(add_count_all(0:numpes-1))
539     
540           IF(PIC_REPORT_SEEDING_STATS) then
541     
542              IF(SUM(ADD_COUNT_ALL(:)).GT.0) THEN
543                 WRITE(err_msg,'(/,2x,A,2x,i10)') &
544                    'TOTAL NUMBER OF PARCELS ADDED GLOBALLY = ',&
545                     SUM(ADD_COUNT_ALL(:))
546     
547                 call flush_err_msg(header = .false., footer = .false.)
548     
549                 DO IPROC = 0, NUMPES-1
550                    IF(DMP_LOG) WRITE(UNIT_LOG, '(/,A,i4,2x,A,2x,i5)') &
551                       'PARCELS ADDED ON PROC:', IPROC,&
552                       ' EQUAL TO', ADD_COUNT_ALL(IPROC)
553                 ENDDO
554              ENDIF
555     
556           ENDIF
557     
558           CALL FINL_ERR_MSG
559           RETURN
560           END SUBROUTINE PIC_MI_BC
561     
562     
563     
564     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
565     !                                                                      C
566     !  Subroutine: PIC_FIND_EMPTY_SPOT                                   C
567     !  Purpose:                                                            C
568     !                                                                      C
569     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
570     
571           SUBROUTINE PIC_FIND_EMPTY_SPOT(LAST_INDEX, EMPTY_SPOT)
572           USE funits
573           USE mpi_utility
574           USE error_manager
575           USE discretelement, only: max_pip, pea
576           IMPLICIT NONE
577     !-----------------------------------------------
578     ! Dummy arguments
579     !-----------------------------------------------
580           INTEGER, INTENT(INOUT) :: LAST_INDEX
581           INTEGER, INTENT(OUT) :: EMPTY_SPOT
582     !-----------------------------------------------
583     ! Local variables
584     !-----------------------------------------------
585           LOGICAL :: SPOT_FOUND
586           INTEGER :: LL
587     !-----------------------------------------------
588     
589           CALL INIT_ERR_MSG("PIC_FIND_EMPTY_SPOT")
590     
591           IF(LAST_INDEX.EQ.MAX_PIP) THEN
592              WRITE(ERR_MSG,2001)
593     
594              CALL FLUSH_ERR_MSG(abort = .true.)
595              call mfix_exit(mype)
596           ENDIF
597           SPOT_FOUND = .false.
598     
599           DO LL = LAST_INDEX+1, MAX_PIP
600     
601              if(.NOT.PEA(LL,1)) THEN
602                 EMPTY_SPOT = LL
603                 LAST_INDEX = LL
604                 SPOT_FOUND = .true.
605                 EXIT
606              ENDIF
607           ENDDO
608     
609           IF(.NOT.SPOT_FOUND) THEN
610              WRITE(ERR_MSG,2002)
611              CALL FLUSH_ERR_MSG(abort = .true.)
612           ENDIF
613     
614      2001 FORMAT(/,5X,  &
615           & 'ERROR IN PIC_FIND_EMPTY_SPOT', /5X, &
616           & 'NO MORE EMPTY SPOT IN THE PARTICLE ARRAY TO ADD A NEW PARTICLE',/5X, &
617           & 'TERMINAL ERROR: STOPPING')
618     
619      2002 FORMAT(/,5X,  &
620           & 'ERROR IN PIC_FIND_EMPTY_SPOT', /5X, &
621           & 'COULD NOT FIND A SPOT FOR ADDING NEW PARTICLE',/5X, &
622           & 'INCREASE THE SIZE OF THE INITIAL ARRAYS', 5X, &
623           & 'TERMINAL ERROR: STOPPING')
624     
625           CALL FINL_ERR_MSG
626           END SUBROUTINE PIC_FIND_EMPTY_SPOT
627     
628     
629     
630     
631     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
632     !                                                                      C
633     !  Subroutine: PIC_REFLECT_PARTICLE                                  C
634     !  Purpose:                                                            C
635     !                                                                      C
636     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
637     
638           SUBROUTINE PIC_REFLECT_PART(LL, WALL_NORM)
639     
640           USE discretelement, only : dimn, DES_VEL_NEW
641           USE mfix_pic, only : MPPIC_COEFF_EN_WALL, &
642           & MPPIC_COEFF_ET_WALL
643           IMPLICIT NONE
644     !-----------------------------------------------
645     ! Dummy arguments
646     !-----------------------------------------------
647           INTEGER, INTENT(IN) :: LL
648           DOUBLE PRECISION, INTENT(IN) ::  WALL_NORM(DIMN)
649     !-----------------------------------------------
650     ! Local variables
651     !-----------------------------------------------
652           !magnitude of pre-collisional normal and tangential velocity components
653           DOUBLE PRECISION :: VEL_NORMMAG_APP
654     
655           !pre collisional normal and tangential velocity components in vector format
656           !APP ==> approach
657           DOUBLE PRECISION :: VEL_NORM_APP(DIMN), VEL_TANG_APP(DIMN)
658     
659     
660           !post collisional normal and tangential velocity components in vector format
661           !SEP ==> separation
662           DOUBLE PRECISION :: VEL_NORM_SEP(DIMN), VEL_TANG_SEP(DIMN)
663     
664           DOUBLE PRECISION :: COEFF_REST_EN, COEFF_REST_ET
665     !-----------------------------------------------
666     
667     
668           VEL_NORMMAG_APP = DOT_PRODUCT(WALL_NORM(1:DIMN), DES_VEL_NEW(:, LL))
669     
670     
671     !currently assuming that wall is at rest. Needs improvement for moving wall
672     
673           VEL_NORM_APP(1:DIMN) = VEL_NORMMAG_APP*WALL_NORM(1:DIMN)
674     
675           VEL_TANG_APP(:) = DES_VEL_NEW(:, LL) - VEL_NORM_APP(1:DIMN)
676     
677     
678     
679     !post collisional velocities
680     
681           !COEFF_REST_EN = REAL_EN_WALL(PIJK(LL,5))
682           COEFF_REST_EN = MPPIC_COEFF_EN_WALL
683           COEFF_REST_ET = MPPIC_COEFF_ET_WALL
684     
685           !if(ep_g(PIJK(LL,4)).lt.0.42) coeff_rest_en = 1.05
686           VEL_NORM_SEP(1:DIMN) = -COEFF_REST_EN*VEL_NORM_APP(1:DIMN)
687     
688           VEL_TANG_SEP(1:DIMN) =  COEFF_REST_ET*VEL_TANG_APP(1:DIMN)
689     
690           DES_VEL_NEW(:, LL) = VEL_NORM_SEP(1:DIMN) + VEL_TANG_SEP(1:DIMN)
691     
692           END SUBROUTINE PIC_REFLECT_PART
693     
694     
695     
696     
697     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
698     !                                                                      C
699     !  Subroutine: PIC_FIND_NEW_CELL                                     C
700     !  Purpose:                                                            C
701     !                                                                      C
702     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
703     
704           SUBROUTINE PIC_FIND_NEW_CELL(LL)
705           USE discretelement
706           use mpi_utility
707           USE cutcell
708           USE functions
709     
710           IMPLICIT NONE
711     
712     !-----------------------------------------------
713     ! Dummy arguments
714     !-----------------------------------------------
715           INTEGER, INTENT(IN) :: LL
716     !-----------------------------------------------
717     ! Local variables
718     !-----------------------------------------------
719           DOUBLE PRECISION :: XPOS, YPOS, ZPOS
720           INTEGER :: I, J, K
721     !-----------------------------------------------
722     
723           I = PIJK(LL,1)
724           J = PIJK(LL,2)
725           K = PIJK(LL,3)
726           XPOS = DES_POS_NEW(1,LL)
727           YPOS = DES_POS_NEW(2,LL)
728           IF (DIMN .EQ. 3) THEN
729              ZPOS = DES_POS_NEW(3,LL)
730           ENDIF
731     
732           IF(XPOS >= XE(I-1) .AND. XPOS < XE(I)) THEN
733              PIJK(LL,1) = I
734           ELSEIF(XPOS >= XE(I)) THEN
735              PIJK(LL,1) = I+1
736           ELSE
737              PIJK(LL,1) = I-1
738           END IF
739     
740           IF(YPOS >= YN(J-1) .AND. YPOS < YN(J)) THEN
741              PIJK(LL,2) = J
742           ELSEIF(YPOS >= YN(J))THEN
743              PIJK(LL,2) = J+1
744           ELSE
745              PIJK(LL,2) = J-1
746           END IF
747     
748           IF(DIMN.EQ.2) THEN
749              PIJK(LL,3) = 1
750           ELSE
751              IF(ZPOS >= ZT(K-1) .AND. ZPOS < ZT(K)) THEN
752                 PIJK(LL,3) = K
753              ELSEIF(ZPOS >= ZT(K)) THEN
754                 PIJK(LL,3) = K+1
755              ELSE
756                 PIJK(LL,3) = K-1
757              END IF
758           ENDIF
759     
760           I = PIJK(LL,1)
761           J = PIJK(LL,2)
762           K = PIJK(LL,3)
763     
764           IF(I.EQ.IEND1+1) then
765              IF(XPOS >= XE(IEND1-1) .AND. XPOS <= XE(IEND1)) PIJK(LL,1) = IEND1
766           ENDIF
767     
768           IF(J.EQ.JEND1+1) then
769              IF(YPOS >= YN(JEND1-1) .AND. YPOS <= YN(JEND1)) PIJK(LL,2) = JEND1
770           ENDIF
771     
772           IF(DIMN.EQ.3.AND.K.EQ.KEND1+1) THEN
773              IF(ZPOS >= ZT(KEND1-1) .AND. ZPOS <= ZT(KEND1)) PIJK(LL,3) = KEND1
774           ENDIF
775     
776           PIJK(LL,4) = FUNIJK(PIJK(LL,1), PIJK(LL,2),PIJK(LL,3))
777     
778           END SUBROUTINE PIC_FIND_NEW_CELL
779     
780     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
781     !                                                                      C
782     !  SUBROUTINE: CHECK_IF_PARCEL_OVELAPS_STL                             C
783     !  Authors: Rahul Garg                               Date: 21-Mar-2014 C
784     !                                                                      C
785     !  Purpose: This subroutine is special written to check if a particle  C
786     !          overlaps any of the STL faces. The routine exits on         C
787     !          detecting an overlap. It is called after initial            C
788     !          generation of lattice configuration to remove out of domain C
789     !          particles                                                   C
790     !                                                                      C
791     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
792           SUBROUTINE CHECK_IF_PARCEL_OVELAPS_STL(POSITION, &
793           OVERLAP_EXISTS)
794     
795           USE constant
796           USE cutcell
797           USE des_stl_functions
798           USE desgrid
799           USE discretelement, only: dimn
800           USE functions
801           USE geometry
802           USE indices
803           USE mpi_utility
804           USE param1
805           USE run
806           USE stl
807           Implicit none
808     
809           DOUBLE PRECISION, INTENT(IN) :: POSITION(DIMN)
810           LOGICAL, INTENT(OUT) :: OVERLAP_EXISTS
811     
812           INTEGER I, J, K, IJK, NF, FOCUS_PARTICLE
813     
814           INTEGER :: COUNT_FAC, COUNT, NEIGH_CELLS, &
815           NEIGH_CELLS_NONNAT, &
816           LIST_OF_CELLS(27), CELL_ID, I_CELL, J_CELL, K_CELL, cell_count , &
817           IMINUS1, IPLUS1, JMINUS1, JPLUS1, KMINUS1, KPLUS1
818     
819           double precision :: velocity(dimn)
820           !reference point and direction of the line
821           double precision, dimension(dimn) :: ref_line,  dir_line
822           !reference point and normal of the plane
823           double precision, dimension(dimn) :: ref_plane, norm_plane
824           !line is parameterized as p = p_ref + t * dir_line, t is line_param
825           double precision :: line_t
826           !logical for determining if a point is on the triangle or not
827           logical :: ontriangle
828           !The line and plane intersection point
829           double precision, dimension(dimn) :: point_onplane
830     
831           FOCUS_PARTICLE = -1
832     
833           OVERLAP_EXISTS = .false.
834     
835           I_CELL = MIN(DG_IEND2,MAX(DG_ISTART2,IOFPOS(POSITION(1))))
836           J_CELL = MIN(DG_JEND2,MAX(DG_JSTART2,JOFPOS(POSITION(2))))
837           K_CELL = 1
838           IF(DO_K) K_CELL =MIN(DG_KEND2,MAX(DG_KSTART2,KOFPOS(POSITION(3))))
839     
840           CELL_ID = DG_FUNIJK(I_CELL, J_CELL, K_CELL)
841           IF (NO_NEIGHBORING_FACET_DES(CELL_ID)) RETURN
842     
843           LIST_OF_CELLS(:) = -1
844           NEIGH_CELLS = 0
845           NEIGH_CELLS_NONNAT  = 0
846     
847           COUNT_FAC = LIST_FACET_AT_DES(CELL_ID)%COUNT_FACETS
848     
849           IF (COUNT_FAC.gt.0)   then
850              !first add the facets in the cell the particle currently resides in
851              NEIGH_CELLS = NEIGH_CELLS + 1
852              LIST_OF_CELLS(NEIGH_CELLS) = CELL_ID
853           ENDIF
854     
855           IPLUS1  =  MIN( I_CELL + 1, DG_IEND2)
856           IMINUS1 =  MAX( I_CELL - 1, DG_ISTART2)
857     
858           JPLUS1  =  MIN (J_CELL + 1, DG_JEND2)
859           JMINUS1 =  MAX( J_CELL - 1, DG_JSTART2)
860     
861           KPLUS1  =  MIN (K_CELL + 1, DG_KEND2)
862           KMINUS1 =  MAX( K_CELL - 1, DG_KSTART2)
863     
864           DO K = KMINUS1, KPLUS1
865              DO J = JMINUS1, JPLUS1
866                 DO I = IMINUS1, IPLUS1
867                    IJK = DG_FUNIJK(I,J,K)
868                    COUNT_FAC = LIST_FACET_AT_DES(IJK)%COUNT_FACETS
869                    IF(COUNT_FAC.EQ.0) CYCLE
870                    NEIGH_CELLS_NONNAT = NEIGH_CELLS_NONNAT + 1
871                    NEIGH_CELLS = NEIGH_CELLS + 1
872                    LIST_OF_CELLS(NEIGH_CELLS) = IJK
873                    !WRITE(*,'(A10, 4(2x,i5))') 'WCELL  = ', IJK, I,J,K
874                 ENDDO
875              ENDDO
876           ENDDO
877     
878     
879           DO CELL_COUNT = 1, NEIGH_CELLS
880              IJK = LIST_OF_CELLS(CELL_COUNT)
881     
882              DO COUNT = 1, LIST_FACET_AT_DES(IJK)%COUNT_FACETS
883                 NF = LIST_FACET_AT_DES(IJK)%FACET_LIST(COUNT)
884                 line_t  = -Undefined
885                 !-undefined, because non zero values will imply intersection
886                 !with the plane
887                 ontriangle = .false.
888     
889                 !parametrize a line as p = p_0 + t normal
890                 !and intersect with the triangular plane.
891                 !if t>0, then point is on the
892                 !non-fluid side of the plane, if the plane normal
893                 !is assumed to point toward the fluid side
894                 ref_line(1:dimn) = position(1:dimn)
895                 dir_line(1:dimn) = NORM_FACE(1:dimn,NF)
896                 !Since this is for checking static config, line's direction
897                 !is the same as plane's normal. For moving particles,
898                 !the line's normal will be along the point joining new
899                 !and old positions.
900     
901                 norm_plane(1:dimn) = NORM_FACE(1:dimn,NF)
902                 ref_plane(1:dimn)  = VERTEX(1, 1:dimn,NF)
903                 CALL intersectLnPlane(ref_line, dir_line, ref_plane, &
904                 norm_plane, line_t)
905                 if(line_t.gt.zero) then
906                    !this implies by orthogonal projection
907                    !that the point is on non-fluid side of the
908                    !facet
909                    point_onplane(1:dimn) = ref_line(1:dimn) + &
910                    line_t*dir_line(1:dimn)
911     !Now check through barycentric coordinates if
912     !this orthogonal projection lies on the facet or not
913     !If it does, then this point is deemed to be on the
914     !non-fluid side of the facet
915                    CALL checkPTonTriangle(point_onplane(1:dimn), &
916                    VERTEX(:,:,NF), ontriangle)
917     
918                    if(ontriangle) then
919                       OVERLAP_EXISTS = .true.
920                       !velocity = zero
921                       !write(*,*) 'over lap detected', line_t
922                       !call write_this_facet_and_parcel(NF, position, velocity)
923                       RETURN
924                    endif
925                 endif
926     
927              ENDDO
928     
929           end DO
930     
931           RETURN
932     
933           END SUBROUTINE CHECK_IF_PARCEL_OVELAPS_STL
934     
935     
936           SUBROUTINE write_this_facet_and_parcel(FID, position, velocity)
937           USE run
938           USE param1
939           USE discretelement
940           USE geometry
941           USE compar
942           USE constant
943           USE cutcell
944           USE funits
945           USE indices
946           USE physprop
947           USE parallel
948           USE stl
949           USE des_stl_functions
950           Implicit none
951           !facet id and particle id
952           double precision, intent(in), dimension(dimn) :: position, velocity
953           Integer, intent(in) :: fid
954           Integer :: stl_unit, vtp_unit , k
955           CHARACTER(LEN=100) :: stl_fname, vtp_fname
956           real :: temp_array(3)
957     
958           stl_unit = 1001
959           vtp_unit = 1002
960     
961           WRITE(vtp_fname,'(A,"_OFFENDING_PARTICLE",".vtp")') TRIM(RUN_NAME)
962           WRITE(stl_fname,'(A,"_STL_FACE",".stl")') TRIM(RUN_NAME)
963     
964           open(vtp_unit, file = vtp_fname, form='formatted')
965           open(stl_unit, file = stl_fname, form='formatted')
966     
967           write(vtp_unit,"(a)") '<?xml version="1.0"?>'
968           write(vtp_unit,"(a,es24.16,a)") '<!-- time =',s_time,'s -->'
969           write(vtp_unit,"(a,a)") '<VTKFile type="PolyData"',&
970                ' version="0.1" byte_order="LittleEndian">'
971           write(vtp_unit,"(3x,a)") '<PolyData>'
972           write(vtp_unit,"(6x,a,i10.10,a,a)")&
973                '<Piece NumberOfPoints="',1,'" NumberOfVerts="0" ',&
974                'NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0">'
975           write(vtp_unit,"(9x,a)")&
976                '<PointData Scalars="Diameter" Vectors="Velocity">'
977           write(vtp_unit,"(12x,a)")&
978                '<DataArray type="Float32" Name="Diameter" format="ascii">'
979           write (vtp_unit,"(15x,es13.6)") (1.d0)
980           write(vtp_unit,"(12x,a)") '</DataArray>'
981     
982           temp_array = zero
983           temp_array(1:DIMN) = velocity(1:dimn)
984           write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
985                'Name="Velocity" NumberOfComponents="3" format="ascii">'
986           write (vtp_unit,"(15x,3(es13.6,3x))")&
987                ((temp_array(k)),k=1,3)
988           write(vtp_unit,"(12x,a,/9x,a)") '</DataArray>','</PointData>'
989           ! skip cell data
990           write(vtp_unit,"(9x,a)") '<CellData></CellData>'
991     
992           temp_array = zero
993           temp_array(1:dimn) = position(1:dimn)
994           write(vtp_unit,"(9x,a)") '<Points>'
995           write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
996                'Name="Position" NumberOfComponents="3" format="ascii">'
997           write (vtp_unit,"(15x,3(es13.6,3x))")&
998                ((temp_array(k)),k=1,3)
999           write(vtp_unit,"(12x,a,/9x,a)")'</DataArray>','</Points>'
1000           ! Write tags for data not included (vtp format style)
1001           write(vtp_unit,"(9x,a,/9x,a,/9x,a,/9x,a)")'<Verts></Verts>',&
1002                '<Lines></Lines>','<Strips></Strips>','<Polys></Polys>'
1003           write(vtp_unit,"(6x,a,/3x,a,/a)")&
1004                '</Piece>','</PolyData>','</VTKFile>'
1005     
1006           !Now write the facet info
1007     
1008           write(stl_unit,*)'solid vcg'
1009     
1010           write(stl_unit,*) '   facet normal ', NORM_FACE(1:3,FID)
1011           write(stl_unit,*) '      outer loop'
1012           write(stl_unit,*) '         vertex ', VERTEX(1,1:3,FID)
1013           write(stl_unit,*) '         vertex ', VERTEX(2,1:3,FID)
1014           write(stl_unit,*) '         vertex ', VERTEX(3,1:3,FID)
1015           write(stl_unit,*) '      endloop'
1016           write(stl_unit,*) '   endfacet'
1017     
1018           write(stl_unit,*)'endsolid vcg'
1019     
1020           close(vtp_unit, status = 'keep')
1021           close(stl_unit, status = 'keep')
1022           write(*,*) 'wrote a facet and a parcel. now waiting'
1023           read(*,*)
1024         end SUBROUTINE write_this_facet_and_parcel
1025     
1026