File: N:\mfix\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     ! Max number of particles in process
16           use discretelement, only: MAX_PIP
17     ! The I/J/K/IJK indicies of the fluid cell
18           use discretelement, only: PIJK
19     ! Number of particles in fluid cells
20           use discretelement, only: PINC
21     ! Run time flag indicating DEM or PIC solids.
22           use run, only: DEM_SOLIDS, PIC_SOLIDS
23     ! Flag for fluid cells
24           use functions, only: FLUID_AT
25     
26           use functions, only: IS_NORMAL
27           use mpi_utility
28           use error_manager
29     
30           IMPLICIT NONE
31     
32     ! Local Variables:
33     !----------------------------------------------------------------------!
34     ! Loop indicies:
35           INTEGER :: L, I, J, K, IJK
36     ! Integer error flag.
37           INTEGER :: IER
38     ! Flag to remove rogue DEM particles
39           LOGICAL, PARAMETER :: REMOVE_ROGUE_PARTICLES = .FALSE.
40     
41     ! Initialize local variables.
42           IER = 0
43     
44     ! Set an error flag if any errors are found. Preform a global collection
45     ! to sync error flags. If needed, reort errors.
46     !.......................................................................
47           DO L = 1, MAX_PIP
48     ! skipping particles that do not exist
49              IF(IS_NORMAL(L)) THEN
50     
51                 I = PIJK(L,1)
52                 J = PIJK(L,2)
53                 K = PIJK(L,3)
54                 IJK = PIJK(L,4)
55     
56                 IF((I > IEND1 .OR. I < ISTART1) .OR. &
57                    (J > JEND1 .OR. J < JSTART1) .OR. &
58                    (K > KEND1 .OR. K < KSTART1)) THEN
59     
60     ! This particle cannot say in the current cell.
61                    PINC(IJK) = PINC(IJK) - 1
62     ! PIC parcels can be relocated.
63                    IF(PIC_SOLIDS) THEN
64                       CALL RECOVER_PARCEL(L)
65     ! DEM particles can be removed.
66                    ELSEIF(REMOVE_ROGUE_PARTICLES) THEN
67                       CALL DELETE_PARTICLE(L)
68     ! Otherwise, this is a fatal error.
69                    ELSE
70                       IER = 1
71                    ENDIF
72     
73                 ENDIF
74              ENDIF
75           ENDDO
76     
77           IF(DEM_SOLIDS) THEN
78              IF(REMOVE_ROGUE_PARTICLES) RETURN
79              CALL GLOBAL_ALL_SUM(IER)
80              IF(IER > 0) CALL CHECK_CELL_MOVEMENT_DEM
81     !     ELSE
82     !        CALL CHECK_CELL_MOVEMENT_PIC
83           ENDIF
84     
85     
86           RETURN
87           END SUBROUTINE CHECK_CELL_MOVEMENT
88     
89     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
90     !                                                                      !
91     !  Subroutine: CHECK_CELL_MOVEMENT_DEM                                 !
92     !                                                                      !
93     !  Purpose: Report which DEM particles have moved into ghost cells.    !
94     !  This is a dead-end routine. Once called, the simulation will exit.  !
95     !                                                                      !
96     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
97           SUBROUTINE CHECK_CELL_MOVEMENT_DEM
98     
99     ! Global Variables:
100     !---------------------------------------------------------------------//
101     ! The global ID of a particle.
102           use discretelement, only: iGlobal_ID
103     ! Max number of particles in process
104           use discretelement, only: MAX_PIP
105     ! The I/J/K/IJK indicies of the fluid cell
106           use discretelement, only: PIJK
107     ! Particle positions and velocities
108           use discretelement, only: DES_POS_NEW, DES_VEL_NEW
109     
110           use functions, only: IS_NORMAL
111           use mpi_utility
112           USE error_manager
113     
114           IMPLICIT NONE
115     
116     ! Local Variables:
117     !----------------------------------------------------------------------!
118     ! Loop indicies:.
119           INTEGER :: L, I, J, K
120     ! Integer error flag
121           INTEGER :: IER
122     
123     
124           CALL INIT_ERR_MSG("CHECK_CELL_MOVEMENT_DEM")
125           CALL OPEN_PE_LOG(IER)
126     
127           WRITE(ERR_MSG, 1100)
128           CALL FLUSH_ERR_MSG(FOOTER=.FALSE.)
129     
130      1100 FORMAT('Error 1100: Particles detected in a ghost cell:',/' ')
131     
132           DO L = 1, MAX_PIP
133     ! skipping particles that do not exist
134              IF(.NOT.IS_NORMAL(L)) CYCLE
135     
136     ! assigning local aliases for particle i, j, k fluid grid indices
137              I = PIJK(L,1)
138              J = PIJK(L,2)
139              K = PIJK(L,3)
140     
141              IF (I.GT.IEND1 .OR. I.LT.ISTART1) THEN
142                 WRITE(ERR_MSG, 1101) trim(iVal(iGlobal_ID(L))),'I',        &
143                    trim(iVal(I)),'X',DES_POS_NEW(L,1),'X',DES_VEL_NEW(L,1)
144                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
145              ENDIF
146     
147              IF(J.GT.JEND1 .OR. J.LT.JSTART1) THEN
148                 WRITE(ERR_MSG, 1101) trim(iVal(iGlobal_id(L))),'J',        &
149                    trim(iVal(J)),'Y',DES_POS_NEW(L,2),'Y',DES_VEL_NEW(L,2)
150                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
151              ENDIF
152     
153              IF (DO_K .AND. (K.GT.KEND1 .OR. K.LT.KSTART1)) THEN
154                 WRITE(ERR_MSG, 1101) trim(iVal(iGlobal_ID(L))),'K',        &
155                    trim(iVal(K)),'Z',DES_POS_NEW(L,3),'Z',DES_VEL_NEW(L,3)
156                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
157              ENDIF
158           ENDDO
159     
160      1101 FORMAT('Particle ',A,' moved into cell with ',A,' index ',A,/    &
161              3x,A,'-Position: ',g11.4,6x,A,'-Velocity:',g11.4,/' ')
162     
163           WRITE(ERR_MSG, 1102)
164           CALL FLUSH_ERR_MSG(HEADER=.FALSE.)
165      1102 FORMAT('This is a fatal error. A particle output file (vtp) ',   &
166              'will be written',/'to aid debugging.')
167     
168     
169           CALL WRITE_DES_DATA
170           CALL MFIX_EXIT(myPE)
171     
172           END SUBROUTINE CHECK_CELL_MOVEMENT_DEM
173     
174     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
175     !                                                                      !
176     !  Subroutine: RECOVER_PARCEL                                          !
177     !                                                                      !
178     !  Purpose: Try to recover the parcel. First move it back to its       !
179     !  previous position, then do a full search to bin the parcel in a     !
180     !  fluid cell. Delete the parcel if it still is located outside of a   !
181     !  fluid cell.                                                         !
182     !                                                                      !
183     !  Possible Future Work:                                               !
184     !                                                                      !
185     !  1. Redistribute particle's weight among other particles in the      !
186     !     domain to conserve mass.                                         !
187     !                                                                      !
188     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
189           SUBROUTINE RECOVER_PARCEL(NP)
190     
191     ! Global Variables:
192     !---------------------------------------------------------------------//
193     ! The (max) number of particles in process
194           use discretelement, only: PIP, MAX_PIP
195     ! The I/J/K/IJK indicies of the fluid cell
196           use discretelement, only: PIJK
197     ! The number of particles in cell IJK
198           use discretelement, only: PINC
199     ! Particle position and velocity
200           use discretelement, only: DES_POS_NEW, DES_VEL_NEW
201     ! The East/North/Top face loctions of fluid cells
202           use discretelement, only: XE, YN, ZT
203     ! Particle velocities
204           use discretelement, only: DTSOLID
205     ! Fluid grid cell dimensions and mesh size
206           USE geometry, only: IMIN2, IMAX2
207           USE geometry, only: JMIN2, JMAX2
208           USE geometry, only: KMIN2, KMAX2
209     ! The East/North/Top face location of a given I/J/K index.
210           use discretelement, only: XE, YN, ZT
211     ! Fixed array sizes in the I/J/K direction
212           use param, only: DIMENSION_I, DIMENSION_J, DIMENSION_K
213     
214     
215           use cutcell
216           use error_manager
217           use functions, only: funijk, fluid_at
218           use mpi_utility
219           use discretelement, only: iglobal_id
220     
221           IMPLICIT NONE
222     
223     
224     ! Dummy Arguments
225     !----------------------------------------------------------------------!
226     ! Parcel ID
227           INTEGER, INTENT(IN) :: NP
228     
229     
230     ! Local Variables:
231     !----------------------------------------------------------------------!
232     ! particle no.
233           INTEGER :: I, J, K, IJK
234     ! Integer error flag.
235           INTEGER :: IER
236     ! Local parameter to print verbose messages about particles.
237           LOGICAL, PARAMETER :: lDEBUG = .FALSE.
238     
239           DOUBLE PRECISION :: oPOS(3)
240     !.......................................................................
241     
242           CALL INIT_ERR_MSG("RECOVER_PARCEL")
243           IF(lDEBUG) CALL OPEN_PE_LOG(IER)
244     
245           oPOS = DES_POS_NEW(NP,:)
246     
247     ! Reflect the parcle.
248           DES_VEL_NEW(NP,:) = -DES_VEL_NEW(NP,:)
249     
250     ! Move the particle back to the previous position.
251           DES_POS_NEW(NP,:) = DES_POS_NEW(NP,:) + &
252              DES_VEL_NEW(NP,:) * DTSOLID
253     
254     ! Rebin the particle to the fluid grid.
255           CALL PIC_SEARCH(I,DES_POS_NEW(NP,1),XE,DIMENSION_I,IMIN2,IMAX2)
256           CALL PIC_SEARCH(J,DES_POS_NEW(NP,2),YN,DIMENSION_J,JMIN2,JMAX2)
257           CALL PIC_SEARCH(K,DES_POS_NEW(NP,3),ZT,DIMENSION_K,KMIN2,KMAX2)
258     
259     ! Calculate the fluid cell index.
260           IJK = FUNIJK(I,J,K)
261           IF(FLUID_AT(IJK)) THEN
262     
263     ! Assign PIJK(L,1:4)
264              PIJK(NP,1) = I
265              PIJK(NP,2) = J
266              PIJK(NP,3) = K
267              PIJK(NP,4) = IJK
268     
269              PINC(IJK) = PINC(IJK) + 1
270           ELSE
271     
272              write(*,*) 'Still not cool -->', iGLOBAL_ID(NP)
273     
274     !         write(*,*) 'POS:',oPOS
275     !         call write_des_data
276     !         stop 'killer stop'
277              CALL DELETE_PARTICLE(NP)
278     
279           ENDIF
280     
281     
282      1100 FORMAT('Warninge 1100: Particle ',A,' was recovered from a ',    &
283              'ghost cell.',2/2x,'Moved into cell with ',A1,' index: ',A,   &
284              /2x,A1,'-Position OLD:',g11.4,/2x,A1,'-Position NEW:',g11.4,  &
285              /2x,A1,'-Velocity:',g11.4)
286     
287      1110 FORMAT('Warninge 1110: Particle ',A,' was deleted from a ',      &
288              'ghost cell.',2/2x,'Moved into cell with ',A1,' index: ',A,   &
289              /2x,A1,'-Position OLD:',g11.4,/2x,A1,'-Position NEW:',g11.4,  &
290              /2X,A1,'-Velocity:',g11.4,/2x,'Fluid Cell: ',A,/2x,           &
291              'Cut cell? ',L1,/2x,'Fluid at? ',L1)
292     
293           CALL FINL_ERR_MSG
294     
295           RETURN
296           END SUBROUTINE RECOVER_PARCEL
297     
298     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
299     !                                                                      !
300     !  Subroutine: CHECK_CELL_MOVEMENT_PIC                                 !
301     !                                                                      !
302     !  Purpose: Report which PIC particles have moved into ghost cells.    !
303     !  Unlike DEM particles, this routine either deletes particles that    !
304     !  are out of the domain, or reassigns the index to try and recover    !
305     !  the particle.                                                       !
306     !                                                                      !
307     !  Notes:                                                              !
308     !                                                                      !
309     !  PIC particles may end up in ghost cells if the cell adjacent to the !
310     !  ghost cell is a cut-cell the particle was not detected outside the  !
311     !  system because of tolerances.                                       !
312     !                                                                      !
313     !  Future Work:                                                        !
314     !                                                                      !
315     !  1. Redistribute particle's weight among other particles in the      !
316     !     domain to conserve mass.                                         !
317     !                                                                      !
318     !  2. Rather than deactivating the particle, reflect the particle      !
319     !     inside the domain using the ghost cell bc's                      !
320     !                                                                      !
321     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
322           SUBROUTINE CHECK_CELL_MOVEMENT_PIC
323     
324     ! Global Variables:
325     !---------------------------------------------------------------------//
326     ! The (max) number of particles in process
327           use discretelement, only: PIP, MAX_PIP
328     ! The I/J/K/IJK indicies of the fluid cell
329           use discretelement, only: PIJK
330     ! The number of particles in cell IJK
331           use discretelement, only: PINC
332     ! Particle positions New/Previous
333           use discretelement, only: DES_POS_NEW, DES_POS_OLD
334     ! The East/North/Top face loctions of fluid cells
335           use discretelement, only: XE, YN, ZT
336     ! Particle velocities
337           use discretelement, only: DES_VEL_NEW
338     ! Flag: identifies a fluid cell as a cut cell.
339           use cutcell, only: CUT_CELL_AT
340     
341           use mpi_utility
342           use error_manager
343           use functions
344     
345           IMPLICIT NONE
346     
347     ! Local Variables:
348     !----------------------------------------------------------------------!
349     ! particle no.
350           INTEGER :: L, I, J, K, IJK
351     ! Integer error flag.
352           INTEGER :: IER
353     ! Position
354           DOUBLE PRECISION :: lPOS
355     ! Number of deleted particles found on the local process
356           INTEGER :: lDELETED, gDELETED
357     ! Number of recovered particles found on the local process
358           INTEGER :: lRECOVERED, gRECOVERED
359     ! Local parameter to print verbose messages about particles.
360           LOGICAL, PARAMETER :: lDEBUG = .FALSE.
361     
362     !.......................................................................
363     
364     
365           CALL INIT_ERR_MSG("CHECK_CELL_MOVEMENT_PIC")
366           IF(lDEBUG) CALL OPEN_PE_LOG(IER)
367     
368     
369     ! Initialize the counters for deleted/recovered parcels.
370           lDELETED = 0
371           lRECOVERED = 0
372     
373           DO L = 1, MAX_PIP
374     ! skipping particles that do not exist
375              IF(.NOT.IS_NORMAL(L)) CYCLE
376     
377     ! assigning local aliases for particle i, j, k fluid grid indices
378              I = PIJK(L,1)
379              J = PIJK(L,2)
380              K = PIJK(L,3)
381     
382     ! this ijk is still an old value as it has not been updated
383              IJK=PIJK(L,4)
384     
385     ! in MPPIC a particle can lie on the surface of the wall as only the
386     ! centers are tracked.
387              IF (I > IEND1 .OR. I < ISTART1) THEN
388     
389                 lPOS = DES_POS_NEW(L,1)
390                 IF(I.EQ.IEND1+1 .AND. &
391                    (lPOS >= XE(IEND1-1) .AND. lPOS <= XE(IEND1)) )THEN
392     
393                    lRECOVERED = lRECOVERED + 1
394                    PIJK(L,1) = IEND1
395     
396                    IF(lDEBUG) THEN
397                       WRITE(ERR_MSG,1100) trim(iVal(L)),'I',trim(iVal(I)), &
398                       'X',DES_POS_OLD(L,1),'X',lPOS,'X',DES_VEL_NEW(L,1)
399                       CALL FLUSH_ERR_MSG
400                    ENDIF
401                 ELSE
402     
403                    lDELETED = lDELETED + 1
404                    CALL SET_NONEXISTENT(L)
405                    PINC(IJK) = PINC(IJK) - 1
406     
407                    IF(lDEBUG) THEN
408                       WRITE(ERR_MSG,1110) trim(iVal(L)),'I',trim(iVal(I)), &
409                       'X',DES_POS_OLD(L,1),'X',lPOS,'X',DES_VEL_NEW(L,1),  &
410                       trim(iVal(IJK)), CUT_CELL_AT(IJK), FLUID_AT(IJK)
411                       CALL FLUSH_ERR_MSG
412                    ENDIF
413                    CYCLE
414                 ENDIF
415              ENDIF
416     
417              IF(J.GT.JEND1 .OR. J.LT.JSTART1) THEN
418                 lPOS = DES_POS_NEW(L,2)
419                 IF(J.EQ.JEND1+1.AND.&
420                   (lPOS >= YN(JEND1-1) .AND. lPOS <= YN(JEND1)) ) THEN
421     
422                    lRECOVERED = lRECOVERED + 1
423                    PIJK(L,2) = JEND1
424     
425                    IF(lDEBUG) THEN
426                       WRITE(ERR_MSG,1100) trim(iVal(L)),'J',trim(iVal(J)), &
427                       'Y',DES_POS_OLD(L,2),'Y',lPOS,'Y',DES_VEL_NEW(L,2)
428                       CALL FLUSH_ERR_MSG
429                    ENDIF
430     
431                 ELSE
432     
433                    lDELETED = lDELETED + 1
434                    CALL SET_NONEXISTENT(L)
435                    PINC(IJK) = PINC(IJK) - 1
436     
437                    IF(lDEBUG) THEN
438                       WRITE(ERR_MSG,1110) trim(iVal(L)),'J',trim(iVal(J)), &
439                       'Y',DES_POS_OLD(L,2),'Y',lPOS,'Y',DES_VEL_NEW(L,2),  &
440                       trim(iVal(IJK)), CUT_CELL_AT(IJK), FLUID_AT(IJK)
441                       CALL FLUSH_ERR_MSG
442                    ENDIF
443                    CYCLE
444                 ENDIF
445              ENDIF
446     
447              IF(DO_K .AND. (K > KEND1 .OR. K < KSTART1)) THEN
448                 lPOS = DES_POS_NEW(L,3)
449                 IF(K == KEND1+1 .AND. &
450                   (lPOS >= ZT(KEND1-1) .AND. lPOS <= ZT(KEND1)) ) THEN
451     
452                    lRECOVERED = lRECOVERED + 1
453                    PIJK(L,3) = KEND1
454     
455                    IF(lDEBUG) THEN
456                       WRITE(ERR_MSG,1100) trim(iVal(L)),'K',trim(iVal(K)), &
457                       'Z',DES_POS_OLD(L,3),'Z',lPOS,'Z',DES_VEL_NEW(L,3)
458                       CALL FLUSH_ERR_MSG
459                    ENDIF
460                 ELSE
461     
462                    lDELETED = lDELETED + 1
463                    CALL SET_NONEXISTENT(L)
464                    PINC(IJK) = PINC(IJK) - 1
465     
466                    IF(lDEBUG) THEN
467                       WRITE(ERR_MSG,1110) trim(iVal(L)),'K',trim(iVal(K)), &
468                       'Z',DES_POS_OLD(L,3),'Z',lPOS,'Z',DES_VEL_NEW(L,3),  &
469                       trim(iVal(IJK)), CUT_CELL_AT(IJK), FLUID_AT(IJK)
470                       CALL FLUSH_ERR_MSG
471                    ENDIF
472                    CYCLE
473                 ENDIF
474              ENDIF
475           ENDDO
476     
477      1100 FORMAT('Warninge 1100: Particle ',A,' was recovered from a ',    &
478              'ghost cell.',2/2x,'Moved into cell with ',A1,' index: ',A,   &
479              /2x,A1,'-Position OLD:',g11.4,/2x,A1,'-Position NEW:',g11.4,  &
480              /2x,A1,'-Velocity:',g11.4)
481     
482      1110 FORMAT('Warninge 1110: Particle ',A,' was deleted from a ',      &
483              'ghost cell.',2/2x,'Moved into cell with ',A1,' index: ',A,   &
484              /2x,A1,'-Position OLD:',g11.4,/2x,A1,'-Position NEW:',g11.4,  &
485              /2X,A1,'-Velocity:',g11.4,/2x,'Fluid Cell: ',A,/2x,           &
486              'Cut cell? ',L1,/2x,'Fluid at? ',L1)
487     
488     ! Update the number of particles
489           PIP = PIP - lDELETED
490     
491     ! Send the numbers to the IO process.
492           CALL GLOBAL_SUM(lRECOVERED, gRECOVERED, PE_IO)
493           CALL GLOBAL_SUM(lDELETED, gDELETED, PE_IO)
494     
495           IF(gRECOVERED + gDELETED > 0) THEN
496              WRITE(ERR_MSG,1115) trim(iVal(gDELETED + gRECOVERED)),        &
497                 trim(iVal(gDELETED)), trim(iVal(gRECOVERED))
498              CALL FLUSH_ERR_MSG
499           ENDIF
500     
501      1115 FORMAT('Warning 1115: ',A,' particles detected outside the ',    &
502              'domain.',/2x,A,' particles were deleted.',/2x,A,' particles',&
503              ' were recovered.')
504     
505           CALL FINL_ERR_MSG
506     
507           RETURN
508           END SUBROUTINE CHECK_CELL_MOVEMENT_PIC
509