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