check_cell_movement.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CHECK_CELL_MOVEMENT                                     !
!                                                                      !
!  Purpose: Check to see if particles have moved into ghost cells.     !
!                                                                      !
!  Note: This routine needs a global communicator to identify errors.  !
!  The collection could get expensive so the call frequency of this    !
!  routine should probably be reduced.                                 !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE check_cell_movement

! Global Variables:
!---------------------------------------------------------------------//
! Max number of particles in process
      use discretelement, only: max_pip
! The I/J/K/IJK indicies of the fluid cell
      use discretelement, only: pijk
! Number of particles in fluid cells
      use discretelement, only: pinc
! Run time flag indicating DEM or PIC solids.
      use run, only: dem_solids, pic_solids
! Flag for fluid cells
      use functions, only: fluid_at

      use functions, only: is_normal
      use mpi_utility
      use error_manager

      IMPLICIT NONE

! Local Variables:
!----------------------------------------------------------------------!
! Loop indicies:
      INTEGER :: L, I, J, K, IJK
! Integer error flag.
      INTEGER :: IER
! Flag to remove rogue DEM particles
      LOGICAL, PARAMETER :: REMOVE_ROGUE_PARTICLES = .false.

! Initialize local variables.
      ier = 0

! Set an error flag if any errors are found. Preform a global collection
! to sync error flags. If needed, reort errors.
!.......................................................................
      DO l = 1, max_pip
! skipping particles that do not exist
         IF(is_normal(l)) THEN

            i = pijk(l,1)
            j = pijk(l,2)
            k = pijk(l,3)
            ijk = pijk(l,4)

            IF((i > iend1 .OR. i < istart1) .OR. &
               (j > jend1 .OR. j < jstart1) .OR. &
               (k > kend1 .OR. k < kstart1)) THEN

! This particle cannot say in the current cell.
               pinc(ijk) = pinc(ijk) - 1
! PIC parcels can be relocated.
               IF(pic_solids) THEN
                  CALL recover_parcel(l)
! DEM particles can be removed.
               ELSEIF(remove_rogue_particles) THEN
                  CALL delete_particle(l)
! Otherwise, this is a fatal error.
               ELSE
                  ier = 1
               ENDIF

            ENDIF
         ENDIF
      ENDDO

      IF(dem_solids) THEN
         IF(remove_rogue_particles) RETURN
         CALL global_all_sum(ier)
         IF(ier > 0) CALL check_cell_movement_dem
!     ELSE
!        CALL CHECK_CELL_MOVEMENT_PIC
      ENDIF


      RETURN
      END SUBROUTINE check_cell_movement

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CHECK_CELL_MOVEMENT_DEM                                 !
!                                                                      !
!  Purpose: Report which DEM particles have moved into ghost cells.    !
!  This is a dead-end routine. Once called, the simulation will exit.  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE check_cell_movement_dem

! Global Variables:
!---------------------------------------------------------------------//
! The global ID of a particle.
      use discretelement, only: iglobal_id
! Max number of particles in process
      use discretelement, only: max_pip
! The I/J/K/IJK indicies of the fluid cell
      use discretelement, only: pijk
! Particle positions and velocities
      use discretelement, only: des_pos_new, des_vel_new

      use functions, only: is_normal
      use mpi_utility
      USE error_manager

      IMPLICIT NONE

! Local Variables:
!----------------------------------------------------------------------!
! Loop indicies:.
      INTEGER :: L, I, J, K
! Integer error flag
      INTEGER :: IER


      CALL init_err_msg("CHECK_CELL_MOVEMENT_DEM")
      CALL open_pe_log(ier)

      WRITE(err_msg, 1100)
      CALL flush_err_msg(footer=.false.)

 1100 FORMAT('Error 1100: Particles detected in a ghost cell:',/' ')

      DO l = 1, max_pip
! skipping particles that do not exist
         IF(.NOT.is_normal(l)) cycle

! assigning local aliases for particle i, j, k fluid grid indices
         i = pijk(l,1)
         j = pijk(l,2)
         k = pijk(l,3)

         IF (i.GT.iend1 .OR. i.LT.istart1) THEN
            WRITE(err_msg, 1101) trim(ival(iglobal_id(l))),'I',        &
               trim(ival(i)),'X',des_pos_new(l,1),'X',des_vel_new(l,1)
            CALL flush_err_msg(header=.false., footer=.false.)
         ENDIF

         IF(j.GT.jend1 .OR. j.LT.jstart1) THEN
            WRITE(err_msg, 1101) trim(ival(iglobal_id(l))),'J',        &
               trim(ival(j)),'Y',des_pos_new(l,2),'Y',des_vel_new(l,2)
            CALL flush_err_msg(header=.false., footer=.false.)
         ENDIF

         IF (do_k .AND. (k.GT.kend1 .OR. k.LT.kstart1)) THEN
            WRITE(err_msg, 1101) trim(ival(iglobal_id(l))),'K',        &
               trim(ival(k)),'Z',des_pos_new(l,3),'Z',des_vel_new(l,3)
            CALL flush_err_msg(header=.false., footer=.false.)
         ENDIF
      ENDDO

 1101 FORMAT('Particle ',a,' moved into cell with ',a,' index ',a,/    &
         3x,a,'-Position: ',g11.4,6x,a,'-Velocity:',g11.4,/' ')

      WRITE(err_msg, 1102)
      CALL flush_err_msg(header=.false.)
 1102 FORMAT('This is a fatal error. A particle output file (vtp) ',   &
         'will be written',/'to aid debugging.')


      CALL write_des_data
      CALL mfix_exit(mype)

      END SUBROUTINE check_cell_movement_dem

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: RECOVER_PARCEL                                          !
!                                                                      !
!  Purpose: Try to recover the parcel. First move it back to its       !
!  previous position, then do a full search to bin the parcel in a     !
!  fluid cell. Delete the parcel if it still is located outside of a   !
!  fluid cell.                                                         !
!                                                                      !
!  Possible Future Work:                                               !
!                                                                      !
!  1. Redistribute particle's weight among other particles in the      !
!     domain to conserve mass.                                         !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE recover_parcel(NP)

! Global Variables:
!---------------------------------------------------------------------//
! The (max) number of particles in process
      use discretelement, only: pip, max_pip
! The I/J/K/IJK indicies of the fluid cell
      use discretelement, only: pijk
! The number of particles in cell IJK
      use discretelement, only: pinc
! Particle position and velocity
      use discretelement, only: des_pos_new, des_vel_new
! The East/North/Top face loctions of fluid cells
      use discretelement, only: xe, yn, zt
! Particle velocities
      use discretelement, only: dtsolid
! Fluid grid cell dimensions and mesh size
      USE geometry, only: imin2, imax2
      USE geometry, only: jmin2, jmax2
      USE geometry, only: kmin2, kmax2
! The East/North/Top face location of a given I/J/K index.
      use discretelement, only: xe, yn, zt
! Fixed array sizes in the I/J/K direction
      use param, only: dimension_i, dimension_j, dimension_k


      use cutcell
      use error_manager
      use functions, only: funijk, fluid_at
      use mpi_utility
      use discretelement, only: iglobal_id

      IMPLICIT NONE


! Dummy Arguments
!----------------------------------------------------------------------!
! Parcel ID
      INTEGER, INTENT(IN) :: NP


! Local Variables:
!----------------------------------------------------------------------!
! particle no.
      INTEGER :: I, J, K, IJK
! Integer error flag.
      INTEGER :: IER
! Local parameter to print verbose messages about particles.
      LOGICAL, PARAMETER :: lDEBUG = .false.

      DOUBLE PRECISION :: oPOS(3)
!.......................................................................

      CALL init_err_msg("RECOVER_PARCEL")
      IF(ldebug) CALL open_pe_log(ier)

      opos = des_pos_new(np,:)

! Reflect the parcle.
      des_vel_new(np,:) = -des_vel_new(np,:)

! Move the particle back to the previous position.
      des_pos_new(np,:) = des_pos_new(np,:) + &
         des_vel_new(np,:) * dtsolid

! Rebin the particle to the fluid grid.
      CALL pic_search(i,des_pos_new(np,1),xe,dimension_i,imin2,imax2)
      CALL pic_search(j,des_pos_new(np,2),yn,dimension_j,jmin2,jmax2)
      CALL pic_search(k,des_pos_new(np,3),zt,dimension_k,kmin2,kmax2)

! Calculate the fluid cell index.
      ijk = funijk(i,j,k)
      IF(fluid_at(ijk)) THEN

! Assign PIJK(L,1:4)
         pijk(np,1) = i
         pijk(np,2) = j
         pijk(np,3) = k
         pijk(np,4) = ijk

         pinc(ijk) = pinc(ijk) + 1
      ELSE

         write(*,*) 'Still not cool -->', iglobal_id(np)

!         write(*,*) 'POS:',oPOS
!         call write_des_data
!         stop 'killer stop'
         CALL delete_particle(np)

      ENDIF


 1100 FORMAT('Warninge 1100: Particle ',a,' was recovered from a ',    &
         'ghost cell.',2/2x,'Moved into cell with ',a1,' index: ',a,   &
         /2x,a1,'-Position OLD:',g11.4,/2x,a1,'-Position NEW:',g11.4,  &
         /2x,a1,'-Velocity:',g11.4)

 1110 FORMAT('Warninge 1110: Particle ',a,' was deleted from a ',      &
         'ghost cell.',2/2x,'Moved into cell with ',a1,' index: ',a,   &
         /2x,a1,'-Position OLD:',g11.4,/2x,a1,'-Position NEW:',g11.4,  &
         /2x,a1,'-Velocity:',g11.4,/2x,'Fluid Cell: ',a,/2x,           &
         'Cut cell? ',l1,/2x,'Fluid at? ',l1)

      CALL finl_err_msg

      RETURN
      END SUBROUTINE recover_parcel

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CHECK_CELL_MOVEMENT_PIC                                 !
!                                                                      !
!  Purpose: Report which PIC particles have moved into ghost cells.    !
!  Unlike DEM particles, this routine either deletes particles that    !
!  are out of the domain, or reassigns the index to try and recover    !
!  the particle.                                                       !
!                                                                      !
!  Notes:                                                              !
!                                                                      !
!  PIC particles may end up in ghost cells if the cell adjacent to the !
!  ghost cell is a cut-cell the particle was not detected outside the  !
!  system because of tolerances.                                       !
!                                                                      !
!  Future Work:                                                        !
!                                                                      !
!  1. Redistribute particle's weight among other particles in the      !
!     domain to conserve mass.                                         !
!                                                                      !
!  2. Rather than deactivating the particle, reflect the particle      !
!     inside the domain using the ghost cell bc's                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE check_cell_movement_pic

! Global Variables:
!---------------------------------------------------------------------//
! The (max) number of particles in process
      use discretelement, only: pip, max_pip
! The I/J/K/IJK indicies of the fluid cell
      use discretelement, only: pijk
! The number of particles in cell IJK
      use discretelement, only: pinc
! Particle positions New/Previous
      use discretelement, only: des_pos_new, des_pos_old
! The East/North/Top face loctions of fluid cells
      use discretelement, only: xe, yn, zt
! Particle velocities
      use discretelement, only: des_vel_new
! Flag: identifies a fluid cell as a cut cell.
      use cutcell, only: cut_cell_at

      use mpi_utility
      use error_manager
      use functions

      IMPLICIT NONE

! Local Variables:
!----------------------------------------------------------------------!
! particle no.
      INTEGER :: L, I, J, K, IJK
! Integer error flag.
      INTEGER :: IER
! Position
      DOUBLE PRECISION :: lPOS
! Number of deleted particles found on the local process
      INTEGER :: lDELETED, gDELETED
! Number of recovered particles found on the local process
      INTEGER :: lRECOVERED, gRECOVERED
! Local parameter to print verbose messages about particles.
      LOGICAL, PARAMETER :: lDEBUG = .false.

!.......................................................................


      CALL init_err_msg("CHECK_CELL_MOVEMENT_PIC")
      IF(ldebug) CALL open_pe_log(ier)


! Initialize the counters for deleted/recovered parcels.
      ldeleted = 0
      lrecovered = 0

      DO l = 1, max_pip
! skipping particles that do not exist
         IF(.NOT.is_normal(l)) cycle

! assigning local aliases for particle i, j, k fluid grid indices
         i = pijk(l,1)
         j = pijk(l,2)
         k = pijk(l,3)

! this ijk is still an old value as it has not been updated
         ijk=pijk(l,4)

! in MPPIC a particle can lie on the surface of the wall as only the
! centers are tracked.
         IF (i > iend1 .OR. i < istart1) THEN

            lpos = des_pos_new(l,1)
            IF(i.EQ.iend1+1 .AND. &
               (lpos >= xe(iend1-1) .AND. lpos <= xe(iend1)) )THEN

               lrecovered = lrecovered + 1
               pijk(l,1) = iend1

               IF(ldebug) THEN
                  WRITE(err_msg,1100) trim(ival(l)),'I',trim(ival(i)), &
                  'X',des_pos_old(l,1),'X',lpos,'X',des_vel_new(l,1)
                  CALL flush_err_msg
               ENDIF
            ELSE

               ldeleted = ldeleted + 1
               CALL set_nonexistent(l)
               pinc(ijk) = pinc(ijk) - 1

               IF(ldebug) THEN
                  WRITE(err_msg,1110) trim(ival(l)),'I',trim(ival(i)), &
                  'X',des_pos_old(l,1),'X',lpos,'X',des_vel_new(l,1),  &
                  trim(ival(ijk)), cut_cell_at(ijk), fluid_at(ijk)
                  CALL flush_err_msg
               ENDIF
               cycle
            ENDIF
         ENDIF

         IF(j.GT.jend1 .OR. j.LT.jstart1) THEN
            lpos = des_pos_new(l,2)
            IF(j.EQ.jend1+1.AND.&
              (lpos >= yn(jend1-1) .AND. lpos <= yn(jend1)) ) THEN

               lrecovered = lrecovered + 1
               pijk(l,2) = jend1

               IF(ldebug) THEN
                  WRITE(err_msg,1100) trim(ival(l)),'J',trim(ival(j)), &
                  'Y',des_pos_old(l,2),'Y',lpos,'Y',des_vel_new(l,2)
                  CALL flush_err_msg
               ENDIF

            ELSE

               ldeleted = ldeleted + 1
               CALL set_nonexistent(l)
               pinc(ijk) = pinc(ijk) - 1

               IF(ldebug) THEN
                  WRITE(err_msg,1110) trim(ival(l)),'J',trim(ival(j)), &
                  'Y',des_pos_old(l,2),'Y',lpos,'Y',des_vel_new(l,2),  &
                  trim(ival(ijk)), cut_cell_at(ijk), fluid_at(ijk)
                  CALL flush_err_msg
               ENDIF
               cycle
            ENDIF
         ENDIF

         IF(do_k .AND. (k > kend1 .OR. k < kstart1)) THEN
            lpos = des_pos_new(l,3)
            IF(k == kend1+1 .AND. &
              (lpos >= zt(kend1-1) .AND. lpos <= zt(kend1)) ) THEN

               lrecovered = lrecovered + 1
               pijk(l,3) = kend1

               IF(ldebug) THEN
                  WRITE(err_msg,1100) trim(ival(l)),'K',trim(ival(k)), &
                  'Z',des_pos_old(l,3),'Z',lpos,'Z',des_vel_new(l,3)
                  CALL flush_err_msg
               ENDIF
            ELSE

               ldeleted = ldeleted + 1
               CALL set_nonexistent(l)
               pinc(ijk) = pinc(ijk) - 1

               IF(ldebug) THEN
                  WRITE(err_msg,1110) trim(ival(l)),'K',trim(ival(k)), &
                  'Z',des_pos_old(l,3),'Z',lpos,'Z',des_vel_new(l,3),  &
                  trim(ival(ijk)), cut_cell_at(ijk), fluid_at(ijk)
                  CALL flush_err_msg
               ENDIF
               cycle
            ENDIF
         ENDIF
      ENDDO

 1100 FORMAT('Warninge 1100: Particle ',a,' was recovered from a ',    &
         'ghost cell.',2/2x,'Moved into cell with ',a1,' index: ',a,   &
         /2x,a1,'-Position OLD:',g11.4,/2x,a1,'-Position NEW:',g11.4,  &
         /2x,a1,'-Velocity:',g11.4)

 1110 FORMAT('Warninge 1110: Particle ',a,' was deleted from a ',      &
         'ghost cell.',2/2x,'Moved into cell with ',a1,' index: ',a,   &
         /2x,a1,'-Position OLD:',g11.4,/2x,a1,'-Position NEW:',g11.4,  &
         /2x,a1,'-Velocity:',g11.4,/2x,'Fluid Cell: ',a,/2x,           &
         'Cut cell? ',l1,/2x,'Fluid at? ',l1)

! Update the number of particles
      pip = pip - ldeleted

! Send the numbers to the IO process.
      CALL global_sum(lrecovered, grecovered, pe_io)
      CALL global_sum(ldeleted, gdeleted, pe_io)

      IF(grecovered + gdeleted > 0) THEN
         WRITE(err_msg,1115) trim(ival(gdeleted + grecovered)),        &
            trim(ival(gdeleted)), trim(ival(grecovered))
         CALL flush_err_msg
      ENDIF

 1115 FORMAT('Warning 1115: ',a,' particles detected outside the ',    &
         'domain.',/2x,a,' particles were deleted.',/2x,a,' particles',&
         ' were recovered.')

      CALL finl_err_msg

      RETURN
      END SUBROUTINE check_cell_movement_pic