mass_inflow_dem.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: DES_MASS_INLET                                          !
!  Author: J.Musser                                   Date: 13-Jul-09  !
!                                                                      !
!  Purpose:  This routine fills in the necessary information for new   !
!  particles entereing the system.                                     !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE mass_inflow_dem

! Modules
!---------------------------------------------------------------------//
      use bc
      use des_allocate
      use des_bc
      use derived_types, only: dg_pic
      use discretelement
      use functions
      use mpi_utility, only: bcast
      implicit none

! Local variables
!---------------------------------------------------------------------//
      INTEGER :: IP, LS, M, NP, IJK, LC
      INTEGER :: BCV, BCV_I
      LOGICAL :: CHECK_FOR_ERRORS, OWNS
! I/J/K index of fluid cell containing the new particle.
      INTEGER :: IJKP(3)
      DOUBLE PRECISION :: DIST, POS(3)
! Random numbers -shared by all ranks-
      DOUBLE PRECISION :: RAND(3)
!......................................................................!

      check_for_errors = .false.

      DO bcv_i = 1, dem_bcmi
         bcv = dem_bcmi_map(bcv_i)

         DO lc=dem_bcmi_ijkstart(bcv_i), dem_bcmi_ijkend(bcv_i)
           ijk = dem_bcmi_ijk(lc)

            DO ls= 1,dg_pic(ijk)%ISIZE
               np = dg_pic(ijk)%P(ls)
               IF(is_exiting(np) .or. is_exiting_ghost(np)) cycle
               SELECT CASE (bc_plane(bcv))
               CASE('N'); dist = des_pos_new(np,2) - yn(bc_j_s(bcv))
               CASE('S'); dist = yn(bc_j_s(bcv)-1) - des_pos_new(np,2)
               CASE('E'); dist = des_pos_new(np,1) - xe(bc_i_w(bcv))
               CASE('W'); dist = xe(bc_i_w(bcv)-1) - des_pos_new(np,1)
               CASE('T'); dist = des_pos_new(np,3) - zt(bc_k_b(bcv))
               CASE('B'); dist = zt(bc_k_b(bcv)-1) - des_pos_new(np,3)
               END SELECT
! The particle is still inside the domain
               IF(dist > des_radius(np)) THEN
                  IF(is_entering(np)) CALL set_normal(np)
                  IF(is_entering_ghost(np)) CALL set_ghost(np)
               ENDIF
            ENDDO
         ENDDO

! All ranks generate random numbers but PE_IO BCASTS its values so that
! the calculated position (with random wiggles) is the same on all ranks.
         CALL random_number(rand)
         CALL bcast(rand)

! Check if any particles need seeded this time step.
         IF(dem_mi_time(bcv_i) > s_time) cycle

         ls = 1

! Loop over the particles being injected
         ploop: DO ip = 1, pi_count(bcv_i)

! Increment the global max particle ID (all ranks).
            imax_global_id = imax_global_id + 1

! Determine the location and phase of the incoming particle.
            CALL seed_new_particle(bcv, bcv_i, rand, m, pos, ijkp, owns)

! Only the rank receiving the new particle needs to continue
            IF(.NOT.owns) cycle ploop

! Increment the number of particle on the processor by one. If the max
! number of particles is exceeded, set the error flag and cycle.
            pip = pip + 1
            CALL particle_grow(pip)
            max_pip = max(pip,max_pip)

! Find the first free space in the particle existance array.
            np_lp: DO np = ls, max_pip
               IF(is_nonexistent(np)) THEN
                  ls = np
                  EXIT np_lp
               ENDIF
            ENDDO np_lp

! Set the particle's global ID.
            iglobal_id(ls) = imax_global_id

! Set the properties of the new particle.
            CALL set_new_particle_props(bcv, m, ls, pos, ijkp)

         ENDDO ploop

! Update the time for seeding the next particle.
         dem_mi_time(bcv_i) = s_time + pi_factor(bcv_i)*dtsolid
! Set the flag for error checking.
         check_for_errors = .true.
      ENDDO

      IF(check_for_errors) THEN
      ENDIF

      RETURN
      END SUBROUTINE mass_inflow_dem

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: SEED_NEW_PARTICLE                                       !
!  Author: J.Musser                                   Date: 14-Aug-09  !
!                                                                      !
!  Purpose:  This routine uses the classification information to place !
!  a new particle in the proper location.                              !
!                                                                      !
!  Comments:                                                           !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE seed_new_particle(lBCV, lBCV_I, lRAND, lM, lPOS, &
         lijkp, lowns)

! Modules
!---------------------------------------------------------------------//
      USE bc
      USE compar
      USE des_bc
      USE desgrid
      USE discretelement
      USE funits
      USE geometry
      USE param1
      USE physprop
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! The associated bc no.
      INTEGER, INTENT(IN) :: lBCV, lBCV_I
! Random numbers
      DOUBLE PRECISION, INTENT(IN) :: lRAND(3)
! Phase of incoming particle.
      INTEGER, INTENT(OUT) :: lM
! Position of incoming particle.
      DOUBLE PRECISION, INTENT(OUT) :: lPOS(3)
! I/J/K index of fluid cell containing the new particle.
      INTEGER, INTENT(OUT) :: lIJKP(3)
! Logical indicating that the current rank is the owner
      LOGICAL, INTENT(OUT) :: lOWNS

! Local variables
!---------------------------------------------------------------------//
! the associated bc no.
!      INTEGER :: BCV
! Random position
      DOUBLE PRECISION RAND1, RAND2
! Array index of vacant position
      INTEGER :: VACANCY
! Number of array position.
      INTEGER :: OCCUPANTS
      INTEGER :: RAND_I
      INTEGER :: lI, lJ, lK
      DOUBLE PRECISION :: WINDOW
!......................................................................!

!      IF(PARTICLE_PLCMNT(lBCV_I) == 'ORDR')THEN
         vacancy = dem_mi(lbcv_i)%VACANCY
         occupants = dem_mi(lbcv_i)%OCCUPANTS
         dem_mi(lbcv_i)%VACANCY = mod(vacancy,occupants) + 1
!      ELSE
!         call bcast(lpar_rad)
!         call bcast(lpar_pos)
!         call bcast(m)
!      ENDIF


! Assign the phase of the incoming particle.
      IF(dem_mi(lbcv_i)%POLYDISPERSE) THEN
         rand_i = ceiling(dble(numfrac_limit)*lrand(1))
         lm = dem_bc_poly_layout(lbcv_i, rand_i)
      ELSE
         lm = dem_bc_poly_layout(lbcv_i,1)
      ENDIF

      window = dem_mi(lbcv_i)%WINDOW
      rand1 = half*d_p0(lm) + (window - d_p0(lm))*lrand(1)
      rand2 = half*d_p0(lm) + (window - d_p0(lm))*lrand(2)


! Set the physical location and I/J/K location of the particle.
      SELECT CASE(bc_plane(lbcv))
      CASE('N','S')

         lpos(1) = dem_mi(lbcv_i)%P(vacancy) + rand1
         lpos(3) = dem_mi(lbcv_i)%Q(vacancy) + rand2
         lpos(2) = dem_mi(lbcv_i)%OFFSET

         lijkp(1) = dem_mi(lbcv_i)%W(vacancy)
         lijkp(3) = dem_mi(lbcv_i)%H(vacancy)
         lijkp(2) = dem_mi(lbcv_i)%L

      CASE('E','W')

         lpos(2) = dem_mi(lbcv_i)%P(vacancy) + rand1
         lpos(3) = dem_mi(lbcv_i)%Q(vacancy) + rand2
         lpos(1) = dem_mi(lbcv_i)%OFFSET

         lijkp(2) = dem_mi(lbcv_i)%W(vacancy)
         lijkp(3) = dem_mi(lbcv_i)%H(vacancy)
         lijkp(1) = dem_mi(lbcv_i)%L

      CASE('T','B')

         lpos(1) = dem_mi(lbcv_i)%P(vacancy) + rand1
         lpos(2) = dem_mi(lbcv_i)%Q(vacancy) + rand2
         lpos(3) = dem_mi(lbcv_i)%OFFSET

         lijkp(1) = dem_mi(lbcv_i)%W(vacancy)
         lijkp(2) = dem_mi(lbcv_i)%H(vacancy)
         lijkp(3) = dem_mi(lbcv_i)%L

      END SELECT

! Easier and cleaner to clear out the third component at the end.
      IF(no_k) lpos(3) = zero

      li = iofpos(lpos(1))
      lj = jofpos(lpos(2))
      lk = kofpos(lpos(3))

      lowns = ((dg_istart <= li) .AND. (li <= dg_iend) .AND.           &
         (dg_jstart <= lj) .AND. (lj <= dg_jend))

      IF(do_k) lowns = lowns .AND. (dg_kstart<=lk) .AND. (lk<=dg_kend)

      RETURN
      END SUBROUTINE seed_new_particle


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: SET_NEW_PARTICLE_PROPS                                  !
!  Author: J.Musser                                   Date: 14-Aug-09  !
!                                                                      !
!  Purpose:  Set the properties of the new particle.                   !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE set_new_particle_props(lBCV, lM, lNP, lPOS, lIJKP)

! Modules
!---------------------------------------------------------------------//
      USE bc
      USE compar
      use constant, only: pi
      USE des_bc
      use desgrid
      use des_rxns
      USE des_thermo
      USE discretelement
      use functions
      USE funits
      USE geometry
      use indices
      USE param1
      USE physprop, only: d_p0, ro_s0, c_ps0
      USE physprop, only: nmax
      use run, only: energy_eq
      use run, only: any_species_eq
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! The associated bc no.
      INTEGER, INTENT(IN) :: lBCV
! Phase of incoming particle.
      INTEGER, INTENT(IN) :: lM
! Index of new particle
      INTEGER, INTENT(IN) :: lNP
! Position of incoming particle.
      DOUBLE PRECISION, INTENT(IN) :: lPOS(3)
! I/J/K index of fluid cell containing the new particle.
      INTEGER, INTENT(IN) :: lIJKP(3)

! Local variables
!---------------------------------------------------------------------//
! I/J/K index of DES grid cell
      INTEGER :: lI, lJ, lK
!......................................................................!


! The particle exists and is entering, not exiting nor a ghost particle
      IF (is_ghost(lnp)) THEN
         CALL set_entering_ghost(lnp)
      ELSE
         CALL set_entering(lnp)
      ENDIF

! Set the initial position values based on mass inlet class
      des_pos_new(lnp,:) = lpos(:)
      ppos(lnp,:) = lpos(:)
      des_vel_new(lnp,1) = bc_u_s(lbcv,lm)
      des_vel_new(lnp,2) = bc_v_s(lbcv,lm)
      des_vel_new(lnp,3) = bc_w_s(lbcv,lm)

! Set the initial velocity values
      IF (do_old) THEN
         des_pos_old(lnp,:) = lpos(:)
         des_vel_old(lnp,:) = des_vel_new(lnp,:)
         omega_old(lnp,:) = 0
      ENDIF

! Set the initial angular velocity values
      omega_new(lnp,:) = 0

! Set the initial angular position values
      IF(particle_orientation) orientation(1:3,lnp) = init_orientation

! Set the particle radius value
      des_radius(lnp) = half * d_p0(lm)

! Set the particle density value
      ro_sol(lnp) = ro_s0(lm)

! Store the I/J/K indices of the particle.
      pijk(lnp,1:3) = lijkp(:)
      pijk(lnp,4) = funijk(lijkp(1), lijkp(2), lijkp(3))

! Set the particle mass phase
      pijk(lnp,5) = lm

! Calculate the DES grid cell indices.
      li = min(dg_iend2,max(dg_istart2,iofpos(des_pos_new(lnp,1))))
      lj = min(dg_jend2,max(dg_jstart2,jofpos(des_pos_new(lnp,2))))
      IF(no_k) THEN
         lk = 1
      ELSE
         lk = min(dg_kend2,max(dg_kstart2,kofpos(des_pos_new(lnp,3))))
      ENDIF
! Store the triple
      dg_pijk(lnp) = dg_funijk(li,lj,lk)

! Calculate the new particle's Volume, Mass, OMOI
      pvol(lnp) = (4.0d0/3.0d0) * pi * des_radius(lnp)**3
      pmass(lnp) = pvol(lnp) * ro_sol(lnp)
      omoi(lnp) = 5.d0 / (2.d0 * pmass(lnp) * des_radius(lnp)**2)

! Clear the drag force
      IF(des_explicitly_coupled) drag_fc(lnp,:) = zero

! If solving the energy equations, set the temperature
      IF(any_species_eq .OR. energy_eq ) des_t_s(lnp) = bc_t_s(lbcv,lm)

! Set species mass fractions
      IF((energy_eq .AND. c_ps0(lm)/=undefined) .OR. any_species_eq)&
         des_x_s(lnp,1:nmax(lm)) = bc_x_s(lbcv,lm,1:nmax(lm))

! Calculate time dependent physical properties
      CALL des_physical_prop(lnp, .false.)


      RETURN
      END SUBROUTINE set_new_particle_props



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine:  DES_NEW_PARTICLE_TEST                                  !
!                                                                      !
!  Purpose:  This routine checks if a new particle placed using the    !
!  random inlet was placed in contact with an existing particle.  If   !
!  so a flag is set indicating contact, and the new particle is        !
!  repositioned within the inlet domain.                               !
!                                                                      !
!  Author: J.Musser                                   Date: 14-Aug-09  !
!                                                                      !
!  Purpose: This routine has to be modified for parallel version       !
!           the parameter now accepts the lpar_rad and lpar_pos and    !
!           tests if it touches any particles                          !
!  Comments:                                                           !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE des_new_particle_test(BCV_I,ppar_rad,ppar_pos,TOUCHING)

! Modules
!---------------------------------------------------------------------//
      USE compar
      USE constant
      USE des_bc
      USE derived_types, only: pic
      USE discretelement
      USE functions
      USE funits
      USE geometry
      USE indices
      USE param1
      USE physprop
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! index of boundary condition
      INTEGER, INTENT(IN) :: BCV_I
      DOUBLE PRECISION, INTENT(IN) :: ppar_pos(dimn)
      DOUBLE PRECISION, INTENT(IN) :: ppar_rad
      LOGICAL, INTENT(INOUT) :: TOUCHING

! Local variables
!---------------------------------------------------------------------//
! particle number id of a potential overlapping/contacting particle
      INTEGER NP2
! total number of particles in current ijk cell and loop counter
      INTEGER NPG, LL
! i, j, k indices along boundary used for loop counters
      INTEGER I, J, K, IJK
! for parallel processing
      integer listart,liend,ljstart,ljend,lkstart,lkend

      DOUBLE PRECISION  DISTVEC(dimn), DIST, R_LM
!......................................................................!

      touching = .false.

! For parallel processing the arrays has to be limited
!      select case (des_mi_class(bcv_i))
!      case ('XW','XE', 'YZw','YZe')
!         listart = gs_array(bcv_i,1)
!         liend = gs_array(bcv_i,2)
!         ljstart = max(gs_array(bcv_i,3),jstart)
!         ljend = min(gs_array(bcv_i,4),jend)
!         lkstart = max(gs_array(bcv_i,5),jstart)
!         lkend = min(gs_array(bcv_i,6),jend)
!      case ('YN','YS', 'XZn','XZs')
!         listart = max(gs_array(bcv_i,1),istart)
!         liend = min(gs_array(bcv_i,2),iend)
!         ljstart = gs_array(bcv_i,3)
!         ljend = gs_array(bcv_i,4)
!         lkstart = max(gs_array(bcv_i,5),jstart)
!         lkend = min(gs_array(bcv_i,6),jend)
!      case ('ZT','ZB', 'XYt','XYb')
!         listart = max(gs_array(bcv_i,1),istart)
!         liend = min(gs_array(bcv_i,2),iend)
!         ljstart = max(gs_array(bcv_i,3),jstart)
!         ljend = min(gs_array(bcv_i,4),jend)
!         lkstart = gs_array(bcv_i,5)
!         lkend = gs_array(bcv_i,6)
!      end select

      DO k = lkstart,lkend
      DO j = ljstart,ljend
      DO i = listart,liend
!      DO K = GS_ARRAY(BCV_I,5), GS_ARRAY(BCV_I,6)
!         DO J = GS_ARRAY(BCV_I,3), GS_ARRAY(BCV_I,4)
!           DO I =  GS_ARRAY(BCV_I,1), GS_ARRAY(BCV_I,2)
             ijk = funijk(i,j,k)
             IF(ASSOCIATED(pic(ijk)%P)) THEN
               npg =  SIZE(pic(ijk)%P)
               DO ll = 1, npg
                  np2 = pic(ijk)%P(ll)
                  distvec(:) = ppar_pos(:) - des_pos_new(np2,:)
                  dist = dot_product(distvec,distvec)
                  r_lm = ppar_rad + des_radius(np2)
                  IF(dist .LE. r_lm*r_lm) touching = .true.
               ENDDO
             ENDIF
           ENDDO
         ENDDO
       ENDDO

      RETURN
      END SUBROUTINE des_new_particle_test