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