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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: CHECK_CELL_MOVEMENT                                     !
4     !                                                                      !
5     !  Purpose: Check to see if particles have moved into ghost cells.     !
6     !                                                                      !
7     !  Note: This routine needs a global communicator to identify errors.  !
8     !  The collection could get expensive so the call frequency of this    !
9     !  routine should probably be reduced.                                 !
10     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
11           SUBROUTINE CHECK_CELL_MOVEMENT
12     
13     ! Global Variables:
14     !---------------------------------------------------------------------//
15     ! Flags for classifying particles
16           use discretelement, only: PEA
17     ! Max number of particles in process
18           use discretelement, only: MAX_PIP
19     ! The I/J/K/IJK indicies of the fluid cell
20           use discretelement, only: PIJK
21     ! Run time flag indicating DEM or PIC solids.
22           use run, only: DEM_SOLIDS, PIC_SOLIDS
23     
24           use mpi_utility
25           use error_manager
26     
27           IMPLICIT NONE
28     
29     ! Local Variables:
30     !----------------------------------------------------------------------!
31     ! Loop indicies:
32           INTEGER :: L, I, J, K
33     ! Integer error flag.
34           INTEGER :: IER
35     
36     
37     ! Initialize local variables.
38           IER = 0
39     
40     ! Set an error flag if any errors are found. Preform a global collection
41     ! to sync error flags. If needed, reort errors.
42     !.......................................................................
43     !!$omp parallel default(shared) private(L, I, J, K, IJK)
44     !!$omp do reduction(+:IER) schedule (guided,50)
45           DO L = 1, MAX_PIP
46     ! skipping particles that do not exist
47              IF(.NOT.PEA(L,1) .OR. any(PEA(L,2:4))) CYCLE
48     
49     ! assigning local aliases for particle i, j, k fluid grid indices
50              I = PIJK(L,1)
51              J = PIJK(L,2)
52              K = PIJK(L,3)
53     
54              IF(I > IEND1 .OR. I < ISTART1) IER = 1
55              IF(J > JEND1 .OR. J < JSTART1) IER = 1
56              IF(DO_K .AND. (K > KEND1 .OR. K < KSTART1)) IER = 1
57           ENDDO
58     !!$omp end parallel
59     
60           CALL GLOBAL_ALL_SUM(IER)
61           IF(IER == 0) RETURN
62     
63           IF(DEM_SOLIDS) CALL CHECK_CELL_MOVEMENT_DEM
64           IF(PIC_SOLIDS) CALL CHECK_CELL_MOVEMENT_PIC
65     
66           RETURN
67           END SUBROUTINE CHECK_CELL_MOVEMENT
68     
69     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
70     !                                                                      !
71     !  Subroutine: CHECK_CELL_MOVEMENT_DEM                                 !
72     !                                                                      !
73     !  Purpose: Report which DEM particles have moved into ghost cells.    !
74     !  This is a dead-end routine. Once called, the simulation will exit.  !
75     !                                                                      !
76     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
77           SUBROUTINE CHECK_CELL_MOVEMENT_DEM
78     
79     ! Global Variables:
80     !---------------------------------------------------------------------//
81     ! The global ID of a particle.
82           use discretelement, only: iGlobal_ID
83     ! Flags for classifying particles
84           use discretelement, only: PEA
85     ! Max number of particles in process
86           use discretelement, only: MAX_PIP
87     ! The I/J/K/IJK indicies of the fluid cell
88           use discretelement, only: PIJK
89     ! Particle positions and velocities
90           use discretelement, only: DES_POS_NEW, DES_VEL_NEW
91     
92           use mpi_utility
93           USE error_manager
94     
95           IMPLICIT NONE
96     
97     ! Local Variables:
98     !----------------------------------------------------------------------!
99     ! Loop indicies:.
100           INTEGER :: L, I, J, K
101     ! Integer error flag
102           INTEGER :: IER
103     
104     
105           CALL INIT_ERR_MSG("CHECK_CELL_MOVEMENT_DEM")
106           CALL OPEN_PE_LOG(IER)
107     
108           WRITE(ERR_MSG, 1100)
109           CALL FLUSH_ERR_MSG(FOOTER=.FALSE.)
110     
111      1100 FORMAT('Error 1100: Particles detected in a ghost cell:',/' ')
112     
113           DO L = 1, MAX_PIP
114     ! skipping particles that do not exist
115              IF(.NOT.PEA(L,1) .OR. any(PEA(L,2:4))) CYCLE
116     
117     ! assigning local aliases for particle i, j, k fluid grid indices
118              I = PIJK(L,1)
119              J = PIJK(L,2)
120              K = PIJK(L,3)
121     
122              IF (I.GT.IEND1 .OR. I.LT.ISTART1) THEN
123                 WRITE(ERR_MSG, 1101) trim(iVal(iGlobal_ID(L))),'I',        &
124                    trim(iVal(I)),'X',DES_POS_NEW(1,L),'X',DES_VEL_NEW(1,L)
125                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
126              ENDIF
127     
128              IF(J.GT.JEND1 .OR. J.LT.JSTART1) THEN
129                 WRITE(ERR_MSG, 1101) trim(iVal(iGlobal_id(L))),'J',        &
130                    trim(iVal(J)),'Y',DES_POS_NEW(2,L),'Y',DES_VEL_NEW(2,L)
131                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
132              ENDIF
133     
134              IF (DO_K .AND. (K.GT.KEND1 .OR. K.LT.KSTART1)) THEN
135                 WRITE(ERR_MSG, 1101) trim(iVal(iGlobal_ID(L))),'K',        &
136                    trim(iVal(K)),'Z',DES_POS_NEW(3,L),'Z',DES_VEL_NEW(3,L)
137                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
138              ENDIF
139           ENDDO
140     
141      1101 FORMAT('Particle ',A,' moved into cell with ',A,' index ',A,/ &
142              3x,A,'-Position: ',g11.4,6x,A,'-Velocity:',g11.4,/' ')
143     
144           WRITE(ERR_MSG, 1102)
145           CALL FLUSH_ERR_MSG(HEADER=.FALSE.)
146      1102 FORMAT('This is a fatal error. A particle output file (vtp) ',   &
147              'will be written',/'to aid debugging.')
148     
149     
150           CALL WRITE_DES_DATA
151           CALL MFIX_EXIT(myPE)
152     
153           END SUBROUTINE CHECK_CELL_MOVEMENT_DEM
154     
155     
156     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
157     !                                                                      !
158     !  Subroutine: CHECK_CELL_MOVEMENT_PIC                                 !
159     !                                                                      !
160     !  Purpose: Report which PIC particles have moved into ghost cells.    !
161     !  Unlike DEM particles, this routine either deletes particles that    !
162     !  are out of the domain, or reassigns the index to try and recover    !
163     !  the particle.                                                       !
164     !                                                                      !
165     !  Notes:                                                              !
166     !                                                                      !
167     !  PIC particles may end up in ghost cells if the cell adjacent to the !
168     !  ghost cell is a cut-cell the particle was not detected outside the  !
169     !  system because of tolerances.                                       !
170     !                                                                      !
171     !  Future Work:                                                        !
172     !                                                                      !
173     !  1. Redistribute particle's weight among other particles in the      !
174     !     domain to conserve mass.                                         !
175     !                                                                      !
176     !  2. Rather than deactivating the particle, reflect the particle      !
177     !     inside the domain using the ghost cell bc's                      !
178     !                                                                      !
179     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
180           SUBROUTINE CHECK_CELL_MOVEMENT_PIC
181     
182     ! Global Variables:
183     !---------------------------------------------------------------------//
184     ! Flags for classifying particles
185           use discretelement, only: PEA
186     ! The (max) number of particles in process
187           use discretelement, only: PIP, MAX_PIP
188     ! The I/J/K/IJK indicies of the fluid cell
189           use discretelement, only: PIJK
190     ! The number of particles in cell IJK
191           use discretelement, only: PINC
192     ! Particle positions New/Previous
193           use discretelement, only: DES_POS_NEW, DES_POS_OLD
194     ! The East/North/Top face loctions of fluid cells
195           use discretelement, only: XE, YN, ZT
196     ! Particle velocities
197           use discretelement, only: DES_VEL_NEW
198     ! Flag: identifies a fluid cell as a cut cell.
199           use cutcell, only: CUT_CELL_AT
200     
201           use mpi_utility
202           use error_manager
203           use functions
204     
205           IMPLICIT NONE
206     
207     ! Local Variables:
208     !----------------------------------------------------------------------!
209     ! particle no.
210           INTEGER :: L, I, J, K, IJK
211     ! Integer error flag.
212           INTEGER :: IER
213     ! Position
214           DOUBLE PRECISION :: lPOS
215     ! Number of deleted particles found on the local process
216           INTEGER :: lDELETED, gDELETED
217     ! Number of recovered particles found on the local process
218           INTEGER :: lRECOVERED, gRECOVERED
219     ! Local parameter to print verbose messages about particles.
220           LOGICAL, PARAMETER :: lDEBUG = .FALSE.
221     
222     !.......................................................................
223     
224     
225           CALL INIT_ERR_MSG("CHECK_CELL_MOVEMENT_PIC")
226           IF(lDEBUG) CALL OPEN_PE_LOG(IER)
227     
228     
229     ! Initialize the counters for deleted/recovered parcels.
230           lDELETED = 0
231           lRECOVERED = 0
232     
233           DO L = 1, MAX_PIP
234     ! skipping particles that do not exist
235              IF(.NOT.PEA(L,1) .OR. any(PEA(L,2:4))) CYCLE
236     
237     ! Assigning local aliases for particle position
238     
239     ! assigning local aliases for particle i, j, k fluid grid indices
240              I = PIJK(L,1)
241              J = PIJK(L,2)
242              K = PIJK(L,3)
243     
244     ! this ijk is still an old value as it has not been updated
245              IJK=PIJK(L,4)
246     
247     ! in MPPIC a particle can lie on the surface of the wall as only the
248     ! centers are tracked.
249              IF (I.GT.IEND1 .OR. I.LT.ISTART1) THEN
250     
251                 lPOS = DES_POS_NEW(1,L)
252                 IF(I.EQ.IEND1+1 .AND. &
253                    (lPOS >= XE(IEND1-1) .AND. lPOS <= XE(IEND1)) )THEN
254     
255                    lRECOVERED = lRECOVERED + 1
256                    PIJK(L,1) = IEND1
257     
258                    IF(lDEBUG) THEN
259                       WRITE(ERR_MSG,1100) trim(iVal(L)),'I',trim(iVal(I)), &
260                       'X',DES_POS_OLD(1,L),'X',lPOS,'X',DES_VEL_NEW(1,L)
261                       CALL FLUSH_ERR_MSG
262                    ENDIF
263                 ELSE
264     
265                    lDELETED = lDELETED + 1
266                    PEA(L,1) = .FALSE.
267                    PINC(IJK) = PINC(IJK) - 1
268     
269                    IF(lDEBUG) THEN
270                       WRITE(ERR_MSG,1110) trim(iVal(L)),'I',trim(iVal(I)), &
271                       'X',DES_POS_OLD(1,L),'X',lPOS,'X',DES_VEL_NEW(1,L),  &
272                       trim(iVal(IJK)), CUT_CELL_AT(IJK), FLUID_AT(IJK)
273                       CALL FLUSH_ERR_MSG
274                    ENDIF
275                    CYCLE
276                 ENDIF
277              ENDIF
278     
279              IF(J.GT.JEND1 .OR. J.LT.JSTART1) THEN
280                 lPOS = DES_POS_NEW(2,L)
281                 IF(J.EQ.JEND1+1.AND.&
282                   (lPOS >= YN(JEND1-1) .AND. lPOS <= YN(JEND1)) ) THEN
283     
284                    lRECOVERED = lRECOVERED + 1
285                    PIJK(L,2) = JEND1
286     
287                    IF(lDEBUG) THEN
288                       WRITE(ERR_MSG,1100) trim(iVal(L)),'J',trim(iVal(J)), &
289                       'Y',DES_POS_OLD(2,L),'Y',lPOS,'Y',DES_VEL_NEW(2,L)
290                       CALL FLUSH_ERR_MSG
291                    ENDIF
292     
293                 ELSE
294     
295                    lDELETED = lDELETED + 1
296                    PEA(L,1) = .FALSE.
297                    PINC(IJK) = PINC(IJK) - 1
298     
299                    IF(lDEBUG) THEN
300                       WRITE(ERR_MSG,1110) trim(iVal(L)),'J',trim(iVal(J)), &
301                       'Y',DES_POS_OLD(2,L),'Y',lPOS,'Y',DES_VEL_NEW(2,L),  &
302                       trim(iVal(IJK)), CUT_CELL_AT(IJK), FLUID_AT(IJK)
303                       CALL FLUSH_ERR_MSG
304                    ENDIF
305                    CYCLE
306                 ENDIF
307              ENDIF
308     
309              IF(DO_K .AND. (K > KEND1 .OR. K < KSTART1)) THEN
310                 lPOS = DES_POS_NEW(3,L)
311                 IF(K == KEND1+1 .AND. &
312                   (lPOS >= ZT(KEND1-1) .AND. lPOS <= ZT(KEND1)) ) THEN
313     
314                    lRECOVERED = lRECOVERED + 1
315                    PIJK(L,3) = KEND1
316     
317                    IF(lDEBUG) THEN
318                       WRITE(ERR_MSG,1100) trim(iVal(L)),'K',trim(iVal(K)), &
319                       'Z',DES_POS_OLD(3,L),'Z',lPOS,'Z',DES_VEL_NEW(3,L)
320                       CALL FLUSH_ERR_MSG
321                    ENDIF
322                 ELSE
323     
324                    lDELETED = lDELETED + 1
325                    PEA(L,1) = .FALSE.
326                    PINC(IJK) = PINC(IJK) - 1
327     
328                    IF(lDEBUG) THEN
329                       WRITE(ERR_MSG,1110) trim(iVal(L)),'K',trim(iVal(K)), &
330                       'Z',DES_POS_OLD(3,L),'Z',lPOS,'Z',DES_VEL_NEW(3,L),  &
331                       trim(iVal(IJK)), CUT_CELL_AT(IJK), FLUID_AT(IJK)
332                       CALL FLUSH_ERR_MSG
333                    ENDIF
334                    CYCLE
335                 ENDIF
336              ENDIF
337           ENDDO
338     
339      1100 FORMAT('Warninge 1100: Particle ',A,' was recovered from a ',    &
340              'ghost cell.',2/2x,'Moved into cell with ',A1,' index: ',A,   &
341              /2x,A1,'-Position OLD:',g11.4,/2x,A1,'-Position NEW:',g11.4,  &
342              /2x,A1,'-Velocity:',g11.4)
343     
344      1110 FORMAT('Warninge 1110: Particle ',A,' was deleted from a ',      &
345              'ghost cell.',2/2x,'Moved into cell with ',A1,' index: ',A,   &
346              /2x,A1,'-Position OLD:',g11.4,/2x,A1,'-Position NEW:',g11.4,  &
347              /2X,A1,'-Velocity:',g11.4,/2x,'Fluid Cell: ',A,/2x,           &
348              'Cut cell? ',L1,/2x,'Fluid at? ',L1)
349     
350     ! Update the number of particles
351           PIP = PIP - lDELETED
352     
353     ! Send the numbers to the IO process.
354           CALL GLOBAL_SUM(lRECOVERED, gRECOVERED, PE_IO)
355           CALL GLOBAL_SUM(lDELETED, gDELETED, PE_IO)
356     
357           IF(gRECOVERED + gDELETED > 0) THEN
358              WRITE(ERR_MSG,1115) trim(iVal(gDELETED + gRECOVERED)),        &
359                 trim(iVal(gDELETED)), trim(iVal(gRECOVERED))
360              CALL FLUSH_ERR_MSG
361           ENDIF
362     
363      1115 FORMAT('Warning 1115: ',A,' particles detected outside the ',    &
364              'domain.',/2x,A,' particles were deleted.',/2x,A,' particles',&
365              ' were recovered.')
366     
367           CALL FINL_ERR_MSG
368     
369           RETURN
370           END SUBROUTINE CHECK_CELL_MOVEMENT_PIC
371     
372     
373     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
374     !                                                                      !
375     !  Subroutine: REPORT_PIC_STATS                                        !
376     !                                                                      !
377     !  Purpose: Output stats about PIC simulation.                         !
378     !                                                                      !
379     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
380           SUBROUTINE REPORT_PIC_STATS
381     
382     ! Global Variables:
383     !---------------------------------------------------------------------//
384     ! Flag to report minimum EP_G
385           use mfix_pic, only: PIC_REPORT_MIN_EPG
386     ! Gas phase volume fraction
387           use fldvar, only: EP_G
388     ! Location of cell faces (East, North, Top)
389           use discretelement, only: XE, YN, ZT
390     
391           use mpi_utility
392           USE error_manager
393           USE functions
394     
395           IMPLICIT NONE
396     
397     ! Local Variables:
398     !----------------------------------------------------------------------!
399     ! Loop counters
400           INTEGER I, J, K, IJK, IPROC
401     
402           INTEGER :: EPg_MIN_loc(0:numpes-1, 4), EPg_MIN_loc2(1)
403           DOUBLE PRECISION :: EPg_MIN(0:numpes-1), EPg_min2
404     
405     !-----------------------------------------------
406     
407           CALL INIT_ERR_MSG("REPORT_PIC_STATS")
408     
409     
410           IF(PIC_REPORT_MIN_EPG) THEN
411     
412              EPG_MIN(:) = 0
413              EPG_MIN(mype) = LARGE_NUMBER
414     
415              EPG_MIN_LOC(:,:) = 0
416              EPG_MIN_LOC(mype,:) = -1
417     
418              DO K = KSTART1, KEND1
419                 DO J = JSTART1, JEND1
420                    DO I = ISTART1, IEND1
421                       IJK = funijk(I,J,K)
422     
423                       IF(EP_G(IJK) < EPG_MIN(mype)) THEN
424                          EPG_MIN_LOC(mype,1) = I
425                          EPG_MIN_LOC(mype,2) = J
426                          EPG_MIN_LOC(mype,3) = K
427                          EPG_MIN_LOC(mype,4) = IJK
428                          EPG_MIN(mype) = EP_G(IJK)
429                       ENDIF
430                    ENDDO
431                 ENDDO
432              ENDDO
433     
434              call GLOBAL_ALL_SUM(EPg_MIN)
435              CALL GLOBAL_ALL_SUM(EPg_MIN_loc)
436     
437              epg_min2     = MINVAL(epg_min(0:numpes-1))
438              epg_min_loc2 = MINLOC(epg_min(0:numpes-1)) - 1
439              !-1, since minloc goes from 1:size of the array.
440              !If not corrected by -1, then the proc id will be off by 1
441     
442              iproc = epg_min_loc2(1)
443     
444              I     = epg_min_loc(iproc, 1)
445              J     = epg_min_loc(iproc, 2)
446              K     = epg_min_loc(iproc, 3)
447              IJK   = epg_min_loc(iproc, 4)
448              WRITE(ERR_MSG,1014) EPG_MIN2, Iproc, I, J, K, IJK, &
449                 XE(I) - 0.5*DX(I), YN(J)-0.5*DY(J), ZT(K) - 0.5*DZ(K)
450     
451      1014       FORMAT( /, &
452                 &      5x,'EPGMIN                    = ', 2x,g17.8,/ &
453                 &      5x,'EPGMIN PROC RANK          = ', 2x, I10, / &
454                 &      5x,'EPGMIN (I, J, K, IJK)     = ', 3(2x,i5),2x,i10,/ &
455                 &      5x,'XMID, YMID, ZMID FOR CELL = ', 3(2x,g17.8))
456     
457                 call flush_err_msg(header = .false., footer = .false.)
458     
459           ENDIF
460     
461           CALL FINL_ERR_MSG
462     
463           RETURN
464           END SUBROUTINE REPORT_PIC_STATS
465     
466     
467