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