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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !  Subroutine: PARTICLES_IN_CELL                                       !
3     !                                                                      !
4     !  Purpose:                                                            !
5     !     - For each particle find the computational fluid cell            !
6     !       containing the particle center.                                !
7     !     - Calculate the bulk density in each computational fluid         !
8     !       cell.                                                          !
9     !     - Calculate the volume average solids velocity in each           !
10     !       computational fluid cell.                                      !
11     !     - For parallel processing indices are altered                    !
12     !                                                                      !
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
14           SUBROUTINE PARTICLES_IN_CELL
15     
16     ! The number of particles on the current process.
17           use discretelement, only: PIP, MAX_PIP
18     ! The I/J/K, IJK, and phase index of each particle
19           use discretelement, only: PIJK
20     ! The number and list of particles in each fluid cell IJK.
21           use derived_types, only: PIC
22           use discretelement, only: PINC
23     ! The East/North/Top face location of a given I/J/K index.
24           use discretelement, only: XE, YN, ZT
25     ! Flag for 2D simulations.
26           use geometry, only: NO_K
27     ! The start and end indices of IJK loops
28           use compar, only: IJKStart3, IJKEnd3
29     ! The Upper and Loper indices covered by the current process.
30           use compar, only: ISTART3, IEND3
31           use compar, only: JSTART3, JEND3
32           use compar, only: KSTART3, KEND3
33     ! Fluid grid cell dimensions and mesh size
34           USE geometry, only: IMIN2, IMAX2
35           USE geometry, only: JMIN2, JMAX2
36           USE geometry, only: KMIN2, KMAX2
37     ! Fixed array sizes in the I/J/K direction
38           use param, only: DIMENSION_I, DIMENSION_J, DIMENSION_K
39           use param, only: DIMENSION_3
40     ! Function to conpute IJK from I/J/K
41           use functions, only: FUNIJK
42     
43           !     The accumulated number of particles in each IJK.
44           use generate_particles, only: particle_count
45     
46           use discretelement, only: DES_POS_NEW
47           use functions, only: IS_NONEXISTENT, IS_GHOST, IS_ENTERING_GHOST, IS_EXITING_GHOST
48     
49           IMPLICIT NONE
50     
51     ! Local Variables
52     !---------------------------------------------------------------------//
53     ! particle no.
54           INTEGER L
55     ! accounted for particles
56           INTEGER PC
57     ! ijk indices
58           INTEGER I, J, K, IJK
59     ! variables that count/store the number of particles in i, j, k cell
60           INTEGER:: npic, pos
61     !......................................................................!
62     
63     ! following quantities are reset every call to particles_in_cell
64           PINC(:) = 0
65     
66     !      allocate(PARTICLE_COUNT(DIMENSION_3))
67     ! Use an incremental approach to determine the new particle location.
68     !-----------------------------------------------------------------------
69     !!$omp parallel default(shared) private(L, I, J, K, IJK)
70     !!$omp do reduction(+:PINC) schedule (guided,50)
71     
72           DO L = 1, MAX_PIP
73     ! skipping particles that do not exist
74              IF(IS_NONEXISTENT(L)) CYCLE
75     
76              I = PIJK(L,1)
77              IF(I <= ISTART3 .OR. I >= IEND3) THEN
78                 CALL PIC_SEARCH(I, DES_POS_NEW(L,1), XE,                   &
79                    DIMENSION_I, IMIN2, IMAX2)
80              ELSE
81                 IF((DES_POS_NEW(L,1) >= XE(I-1)) .AND.                     &
82                    (DES_POS_NEW(L,1) <  XE(I))) THEN
83                    I = I
84                 ELSEIF((DES_POS_NEW(L,1) >= XE(I)) .AND.                   &
85                    (DES_POS_NEW(L,1) < XE(I+1))) THEN
86                   I = I+1
87                 ELSEIF((DES_POS_NEW(L,1) >= XE(I-2)) .AND.                 &
88                    (DES_POS_NEW(L,1) < XE(I-1))) THEN
89                    I = I-1
90                 ELSE
91                    CALL PIC_SEARCH(I, DES_POS_NEW(L,1), XE,                &
92                       DIMENSION_I, IMIN2, IMAX2)
93                 ENDIF
94              ENDIF
95     
96              J = PIJK(L,2)
97              IF(J <= JSTART3 .OR. J >= JEND3) THEN
98                 CALL PIC_SEARCH(J, DES_POS_NEW(L,2), YN,                   &
99                    DIMENSION_J, JMIN2, JMAX2)
100              ELSE
101                 IF((DES_POS_NEW(L,2) >= YN(J-1)) .AND.                     &
102                    (DES_POS_NEW(L,2) < YN(J))) THEN
103                    J = J
104                 ELSEIF((DES_POS_NEW(L,2) >= YN(J)) .AND.                   &
105                    (DES_POS_NEW(L,2) < YN(J+1))) THEN
106                    J = J+1
107                 ELSEIF((DES_POS_NEW(L,2) >= YN(J-2)) .AND.                 &
108                    (DES_POS_NEW(L,2) < YN(J-1)))THEN
109                    J = J-1
110                 ELSE
111                    CALL PIC_SEARCH(J, DES_POS_NEW(L,2), YN,                &
112                       DIMENSION_J, JMIN2, JMAX2)
113                 ENDIF
114              ENDIF
115     
116     
117              IF(NO_K) THEN
118                 K = 1
119              ELSE
120                 K = PIJK(L,3)
121                 IF(K <= KSTART3 .OR. K >= KEND3) THEN
122                    CALL PIC_SEARCH(K, DES_POS_NEW(L,3), ZT,                &
123                       DIMENSION_K, KMIN2, KMAX2)
124                 ELSE
125                    IF((DES_POS_NEW(L,3) >= ZT(K-1)) .AND.                  &
126                       (DES_POS_NEW(L,3) < ZT(K))) THEN
127                       K = K
128                     ELSEIF((DES_POS_NEW(L,3) >= ZT(K)) .AND.               &
129                       (DES_POS_NEW(L,3) < ZT(K+1))) THEN
130                       K = K+1
131                    ELSEIF((DES_POS_NEW(L,3) >= ZT(K-2)) .AND.              &
132                       (DES_POS_NEW(L,3) < ZT(K-1))) THEN
133                       K = K-1
134                    ELSE
135                       CALL PIC_SEARCH(K, DES_POS_NEW(L,3), ZT,             &
136                          DIMENSION_K, KMIN2, KMAX2)
137                    ENDIF
138                 ENDIF
139              ENDIF
140     
141     ! Calculate the fluid cell index.
142              IJK = FUNIJK(I,J,K)
143     
144     ! Assign PIJK(L,1:4)
145              PIJK(L,1) = I
146              PIJK(L,2) = J
147              PIJK(L,3) = K
148              PIJK(L,4) = IJK
149     
150     ! Increment the number of particles in cell IJK
151              IF(.NOT.IS_GHOST(L) .AND. .NOT.IS_ENTERING_GHOST(L) .AND. &
152                 .NOT.IS_EXITING_GHOST(L)) PINC(IJK) = PINC(IJK) + 1
153     
154           ENDDO
155     !!$omp end parallel
156     
157           CALL CHECK_CELL_MOVEMENT
158     
159     ! Assigning the variable PIC(IJK)%p(:). For each computational fluid
160     ! cell compare the number of current particles in the cell to what was
161     ! in the cell previously. If different reallocate. Store the particle
162     ! ids
163     ! ---------------------------------------------------------------->>>
164     !!$omp parallel do if(ijkend3 .ge. 2000) default(shared)           &
165     !!$omp private(ijk,npic) !schedule (guided,50)
166           DO IJK = IJKSTART3, IJKEND3
167     
168     ! checking all cells (including ghost cells); updating entering/exiting
169     ! particle regions
170              NPIC =  PINC(IJK)
171              IF (ASSOCIATED(PIC(IJK)%p)) THEN
172                 IF (NPIC.NE.SIZE(PIC(IJK)%p)) THEN
173                    DEALLOCATE(PIC(IJK)%p)
174                    IF (NPIC.GT.0) ALLOCATE(PIC(IJK)%p(NPIC))
175                 ENDIF
176              ELSE
177                 IF (NPIC.GT.0) ALLOCATE(PIC(IJK)%p(NPIC))
178              ENDIF
179           ENDDO
180     !!$omp end parallel do
181     
182           PARTICLE_COUNT(:) = 1
183           PC = 1
184           DO L = 1, MAX_PIP
185     ! exiting loop if reached max number of particles in processor
186              IF(PC.GT.PIP) exit
187     ! skipping indices with no particles (non-existent particles)
188              IF(IS_NONEXISTENT(L)) CYCLE
189     ! incrementing particle account when particle exists
190              PC = PC+1
191     ! skipping ghost particles
192              IF(IS_GHOST(L) .OR. IS_ENTERING_GHOST(L) .OR. IS_EXITING_GHOST(L)) CYCLE
193              IJK = PIJK(L,4)
194              POS = PARTICLE_COUNT(IJK)
195              PIC(IJK)%P(POS) = L
196              PARTICLE_COUNT(IJK) = PARTICLE_COUNT(IJK) + 1
197           ENDDO
198     
199     !      deallocate(PARTICLE_COUNT)
200     
201           RETURN
202           END SUBROUTINE PARTICLES_IN_CELL
203     
204     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
205     !                                                                      !
206     !  Subroutine: INIT_PARTICLES_IN_CELL                                  !
207     !                                                                      !
208     !  Purpose:                                                            !
209     !     - For each particle find the computational fluid cell            !
210     !       containing the particle center.                                !
211     !     - Calculate the bulk density in each computational fluid         !
212     !       cell.                                                          !
213     !     - Calculate the volume average solids velocity in each           !
214     !       computational fluid cell.                                      !
215     !     - For parallel processing indices are altered                    !
216     !                                                                      !
217     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
218           SUBROUTINE INIT_PARTICLES_IN_CELL
219     
220           use discretelement, only: PIJK, PINC
221           USE discretelement, only: DES_POS_NEW
222           USE discretelement, only: MAX_PIP
223           USE discretelement, only: XE, YN, ZT
224           USE functions, only: IS_NONEXISTENT, IS_GHOST, IS_ENTERING_GHOST, IS_EXITING_GHOST
225           use mpi_funs_des, only: des_par_exchange
226     
227     ! Number of particles in the I/J/K direction
228           use param, only: DIMENSION_I, DIMENSION_J, DIMENSION_K
229     
230           use mpi_utility
231           use sendrecv
232     
233           USE error_manager
234           USE desgrid, only: desgrid_pic
235     
236           IMPLICIT NONE
237     !-----------------------------------------------
238     ! Local Variables
239     !-----------------------------------------------
240     ! particle no.
241           INTEGER :: L
242     ! ijk indices
243           INTEGER :: I, J, K, IJK
244     
245           CALL INIT_ERR_MSG("INIT_PARTICLES_IN_CELL")
246     
247     ! following quantities are reset every call to particles_in_cell
248           PINC(:) = 0
249     
250     ! Bin the particles to the DES grid.
251           CALL DESGRID_PIC(.TRUE.)
252     ! Call exchange particles - this will exchange particle crossing
253     ! boundaries as well as updates ghost particles information
254           CALL DES_PAR_EXCHANGE
255     
256     ! Assigning PIJK(L,1), PIJK(L,2) and PIJK(L,3) the i, j, k indices
257     ! of particle L (locating particle on fluid grid). Also determine
258     ! composite ijk index. If first_pass, also assigning PIJK(L,5) the
259     ! solids phase index of particle.
260     ! ---------------------------------------------------------------->>>
261           DO L = 1, MAX_PIP
262     ! skipping particles that do not exist
263              IF(IS_NONEXISTENT(L)) CYCLE
264     
265     ! Use a brute force technique to determine the particle locations in
266     ! the Eulerian fluid grid.
267     
268              CALL PIC_SEARCH(I, DES_POS_NEW(L,1), XE,                      &
269                 DIMENSION_I, IMIN2, IMAX2)
270              PIJK(L,1) = I
271     
272              CALL PIC_SEARCH(J, DES_POS_NEW(L,2), YN,                      &
273                 DIMENSION_J, JMIN2, JMAX2)
274              PIJK(L,2) = J
275     
276              IF(NO_K) THEN
277                 K=1
278                 PIJK(L,3) = 1
279              ELSE
280                 CALL PIC_SEARCH(K, DES_POS_NEW(L,3), ZT,                   &
281                    DIMENSION_K, KMIN2, KMAX2)
282                 PIJK(L,3) = K
283              ENDIF
284     
285     ! Assigning PIJK(L,4) now that particles have been located on the fluid
286              IJK = FUNIJK(I,J,K)
287              PIJK(L,4) = IJK
288     
289     ! Enumerate the number of 'real' particles in the ghost cell.
290              IF(.NOT.IS_GHOST(L) .AND. .NOT.IS_ENTERING_GHOST(L) .AND. &
291                 .NOT.IS_EXITING_GHOST(L)) PINC(IJK) = PINC(IJK) + 1
292           ENDDO
293     
294     ! Bin the particles to the DES grid.
295           CALL DESGRID_PIC(.TRUE.)
296     ! Calling exchange particles - this will exchange particle crossing
297     ! boundaries as well as updates ghost particles information
298     ! unclear why this needs to be called again.
299           CALL DES_PAR_EXCHANGE
300     
301           CALL FINL_ERR_MSG
302     
303           RETURN
304           END SUBROUTINE INIT_PARTICLES_IN_CELL
305     
306     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
307     !                                                                      !
308     !  Subroutine: PIC_SEARCH                                              !
309     !                                                                      !
310     !  Purpose: Identify the I (or J or K) index of the fluid cell that    !
311     !  contains the particle centroid.                                     !
312     !                                                                      !
313     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
314           SUBROUTINE PIC_SEARCH(IDX, lPOS, ENT_POS, lDIMN, lSTART, lEND)
315     
316           IMPLICIT NONE
317     !-----------------------------------------------
318     ! Local Variables
319     !-----------------------------------------------
320     ! Index being searched for (I, J, or K)
321           INTEGER, INTENT(OUT) :: IDX
322     ! Particle x,y,z position
323           DOUBLE PRECISION, INTENT(IN) :: lPOS
324     ! Dimension of ENT_POS array
325           INTEGER, INTENT(IN) :: lDIMN
326     ! East, North, or Top cell face location
327           DOUBLE PRECISION, INTENT(IN) :: ENT_POS(0:lDIMN)
328     ! Search bounds (by rank)
329           INTEGER, INTENT(IN) :: lSTART, lEND
330     
331           DO IDX = lSTART,lEND
332              IF(lPOS >= ENT_POS(IDX-1) .AND. lPOS < ENT_POS(IDX)) EXIT
333           ENDDO
334     
335           RETURN
336           END SUBROUTINE PIC_SEARCH
337