particles_in_cell.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: PARTICLES_IN_CELL                                       !
!                                                                      !
!  Purpose:                                                            !
!     - For each particle find the computational fluid cell            !
!       containing the particle center.                                !
!     - Calculate the bulk density in each computational fluid         !
!       cell.                                                          !
!     - Calculate the volume average solids velocity in each           !
!       computational fluid cell.                                      !
!     - For parallel processing indices are altered                    !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE particles_in_cell

! The number of particles on the current process.
      use discretelement, only: pip, max_pip
! The I/J/K, IJK, and phase index of each particle
      use discretelement, only: pijk
! The number and list of particles in each fluid cell IJK.
      use derived_types, only: pic
      use discretelement, only: pinc
! The East/North/Top face location of a given I/J/K index.
      use discretelement, only: xe, yn, zt
! Flag for 2D simulations.
      use geometry, only: no_k
! The start and end indices of IJK loops
      use compar, only: ijkstart3, ijkend3
! The Upper and Loper indices covered by the current process.
      use compar, only: istart3, iend3
      use compar, only: jstart3, jend3
      use compar, only: kstart3, kend3
! Fluid grid cell dimensions and mesh size
      USE geometry, only: imin2, imax2
      USE geometry, only: jmin2, jmax2
      USE geometry, only: kmin2, kmax2
! Fixed array sizes in the I/J/K direction
      use param, only: dimension_i, dimension_j, dimension_k
      use param, only: dimension_3
! Function to conpute IJK from I/J/K
      use functions, only: funijk

      !     The accumulated number of particles in each IJK.
      use generate_particles, only: particle_count

      use discretelement, only: des_pos_new
      use functions, only: is_nonexistent, is_ghost, is_entering_ghost, is_exiting_ghost

      IMPLICIT NONE

! Local Variables
!---------------------------------------------------------------------//
! particle no.
      INTEGER L
! accounted for particles
      INTEGER PC
! ijk indices
      INTEGER I, J, K, IJK
! variables that count/store the number of particles in i, j, k cell
      INTEGER:: npic, pos
!......................................................................!

! following quantities are reset every call to particles_in_cell
      pinc(:) = 0

!      allocate(PARTICLE_COUNT(DIMENSION_3))
! Use an incremental approach to determine the new particle location.
!-----------------------------------------------------------------------
!!$omp parallel default(shared) private(L, I, J, K, IJK)
!!$omp do reduction(+:PINC) schedule (guided,50)

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

         i = pijk(l,1)
         IF(i <= istart3 .OR. i >= iend3) THEN
            CALL pic_search(i, des_pos_new(l,1), xe,                   &
               dimension_i, imin2, imax2)
         ELSE
            IF((des_pos_new(l,1) >= xe(i-1)) .AND.                     &
               (des_pos_new(l,1) <  xe(i))) THEN
               i = i
            ELSEIF((des_pos_new(l,1) >= xe(i)) .AND.                   &
               (des_pos_new(l,1) < xe(i+1))) THEN
              i = i+1
            ELSEIF((des_pos_new(l,1) >= xe(i-2)) .AND.                 &
               (des_pos_new(l,1) < xe(i-1))) THEN
               i = i-1
            ELSE
               CALL pic_search(i, des_pos_new(l,1), xe,                &
                  dimension_i, imin2, imax2)
            ENDIF
         ENDIF

         j = pijk(l,2)
         IF(j <= jstart3 .OR. j >= jend3) THEN
            CALL pic_search(j, des_pos_new(l,2), yn,                   &
               dimension_j, jmin2, jmax2)
         ELSE
            IF((des_pos_new(l,2) >= yn(j-1)) .AND.                     &
               (des_pos_new(l,2) < yn(j))) THEN
               j = j
            ELSEIF((des_pos_new(l,2) >= yn(j)) .AND.                   &
               (des_pos_new(l,2) < yn(j+1))) THEN
               j = j+1
            ELSEIF((des_pos_new(l,2) >= yn(j-2)) .AND.                 &
               (des_pos_new(l,2) < yn(j-1)))THEN
               j = j-1
            ELSE
               CALL pic_search(j, des_pos_new(l,2), yn,                &
                  dimension_j, jmin2, jmax2)
            ENDIF
         ENDIF


         IF(no_k) THEN
            k = 1
         ELSE
            k = pijk(l,3)
            IF(k <= kstart3 .OR. k >= kend3) THEN
               CALL pic_search(k, des_pos_new(l,3), zt,                &
                  dimension_k, kmin2, kmax2)
            ELSE
               IF((des_pos_new(l,3) >= zt(k-1)) .AND.                  &
                  (des_pos_new(l,3) < zt(k))) THEN
                  k = k
                ELSEIF((des_pos_new(l,3) >= zt(k)) .AND.               &
                  (des_pos_new(l,3) < zt(k+1))) THEN
                  k = k+1
               ELSEIF((des_pos_new(l,3) >= zt(k-2)) .AND.              &
                  (des_pos_new(l,3) < zt(k-1))) THEN
                  k = k-1
               ELSE
                  CALL pic_search(k, des_pos_new(l,3), zt,             &
                     dimension_k, kmin2, kmax2)
               ENDIF
            ENDIF
         ENDIF

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

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

! Increment the number of particles in cell IJK
         IF(.NOT.is_ghost(l) .AND. .NOT.is_entering_ghost(l) .AND. &
            .NOT.is_exiting_ghost(l)) pinc(ijk) = pinc(ijk) + 1

      ENDDO
!!$omp end parallel

      CALL check_cell_movement

! Assigning the variable PIC(IJK)%p(:). For each computational fluid
! cell compare the number of current particles in the cell to what was
! in the cell previously. If different reallocate. Store the particle
! ids
! ---------------------------------------------------------------->>>
!!$omp parallel do if(ijkend3 .ge. 2000) default(shared)           &
!!$omp private(ijk,npic) !schedule (guided,50)
      DO ijk = ijkstart3, ijkend3

! checking all cells (including ghost cells); updating entering/exiting
! particle regions
         npic =  pinc(ijk)
         IF (ASSOCIATED(pic(ijk)%p)) THEN
            IF (npic.NE.SIZE(pic(ijk)%p)) THEN
               DEALLOCATE(pic(ijk)%p)
               IF (npic.GT.0) ALLOCATE(pic(ijk)%p(npic))
            ENDIF
         ELSE
            IF (npic.GT.0) ALLOCATE(pic(ijk)%p(npic))
         ENDIF
      ENDDO
!!$omp end parallel do

      particle_count(:) = 1
      pc = 1
      DO l = 1, max_pip
! exiting loop if reached max number of particles in processor
         IF(pc.GT.pip) exit
! skipping indices with no particles (non-existent particles)
         IF(is_nonexistent(l)) cycle
! incrementing particle account when particle exists
         pc = pc+1
! skipping ghost particles
         IF(is_ghost(l) .OR. is_entering_ghost(l) .OR. is_exiting_ghost(l)) cycle
         ijk = pijk(l,4)
         pos = particle_count(ijk)
         pic(ijk)%P(pos) = l
         particle_count(ijk) = particle_count(ijk) + 1
      ENDDO

!      deallocate(PARTICLE_COUNT)

      RETURN
      END SUBROUTINE particles_in_cell

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: INIT_PARTICLES_IN_CELL                                  !
!                                                                      !
!  Purpose:                                                            !
!     - For each particle find the computational fluid cell            !
!       containing the particle center.                                !
!     - Calculate the bulk density in each computational fluid         !
!       cell.                                                          !
!     - Calculate the volume average solids velocity in each           !
!       computational fluid cell.                                      !
!     - For parallel processing indices are altered                    !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE init_particles_in_cell

      use discretelement, only: pijk, pinc
      USE discretelement, only: des_pos_new
      USE discretelement, only: max_pip
      USE discretelement, only: xe, yn, zt
      USE functions, only: is_nonexistent, is_ghost, is_entering_ghost, is_exiting_ghost
      use mpi_funs_des, only: des_par_exchange

! Number of particles in the I/J/K direction
      use param, only: dimension_i, dimension_j, dimension_k

      use mpi_utility
      use sendrecv

      USE error_manager
      USE desgrid, only: desgrid_pic

      IMPLICIT NONE
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! particle no.
      INTEGER :: L
! ijk indices
      INTEGER :: I, J, K, IJK

      CALL init_err_msg("INIT_PARTICLES_IN_CELL")

! following quantities are reset every call to particles_in_cell
      pinc(:) = 0

! Bin the particles to the DES grid.
      CALL desgrid_pic(.true.)
! Call exchange particles - this will exchange particle crossing
! boundaries as well as updates ghost particles information
      CALL des_par_exchange

! Assigning PIJK(L,1), PIJK(L,2) and PIJK(L,3) the i, j, k indices
! of particle L (locating particle on fluid grid). Also determine
! composite ijk index. If first_pass, also assigning PIJK(L,5) the
! solids phase index of particle.
! ---------------------------------------------------------------->>>
      DO l = 1, max_pip
! skipping particles that do not exist
         IF(is_nonexistent(l)) cycle

! Use a brute force technique to determine the particle locations in
! the Eulerian fluid grid.

         CALL pic_search(i, des_pos_new(l,1), xe,                      &
            dimension_i, imin2, imax2)
         pijk(l,1) = i

         CALL pic_search(j, des_pos_new(l,2), yn,                      &
            dimension_j, jmin2, jmax2)
         pijk(l,2) = j

         IF(no_k) THEN
            k=1
            pijk(l,3) = 1
         ELSE
            CALL pic_search(k, des_pos_new(l,3), zt,                   &
               dimension_k, kmin2, kmax2)
            pijk(l,3) = k
         ENDIF

! Assigning PIJK(L,4) now that particles have been located on the fluid
         ijk = funijk(i,j,k)
         pijk(l,4) = ijk

! Enumerate the number of 'real' particles in the ghost cell.
         IF(.NOT.is_ghost(l) .AND. .NOT.is_entering_ghost(l) .AND. &
            .NOT.is_exiting_ghost(l)) pinc(ijk) = pinc(ijk) + 1
      ENDDO

! Bin the particles to the DES grid.
      CALL desgrid_pic(.true.)
! Calling exchange particles - this will exchange particle crossing
! boundaries as well as updates ghost particles information
! unclear why this needs to be called again.
      CALL des_par_exchange

      CALL finl_err_msg

      RETURN
      END SUBROUTINE init_particles_in_cell

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: PIC_SEARCH                                              !
!                                                                      !
!  Purpose: Identify the I (or J or K) index of the fluid cell that    !
!  contains the particle centroid.                                     !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE pic_search(IDX, lPOS, ENT_POS, lDIMN, lSTART, lEND)

      IMPLICIT NONE
!-----------------------------------------------
! Local Variables
!-----------------------------------------------
! Index being searched for (I, J, or K)
      INTEGER, INTENT(OUT) :: IDX
! Particle x,y,z position
      DOUBLE PRECISION, INTENT(IN) :: lPOS
! Dimension of ENT_POS array
      INTEGER, INTENT(IN) :: lDIMN
! East, North, or Top cell face location
      DOUBLE PRECISION, INTENT(IN) :: ENT_POS(0:ldimn)
! Search bounds (by rank)
      INTEGER, INTENT(IN) :: lSTART, lEND

      DO idx = lstart,lend
         IF(lpos >= ent_pos(idx-1) .AND. lpos < ent_pos(idx)) EXIT
      ENDDO

      RETURN
      END SUBROUTINE pic_search