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