mass_outflow_dem.f

Go to the documentation of this file.

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

      use derived_types, only: dg_pic
      use discretelement
      use des_bc
      use bc
      use functions
      use param1, only: zero

      use mpi_utility, only: global_all_or

      implicit none

      LOGICAL, INTENT(INOUT) :: FORCE_NSEARCH

      INTEGER :: IJK
      INTEGER :: LC, LP, NP, M
      INTEGER :: BCV, BCV_I, IDX

      DOUBLE PRECISION :: SGN
      DOUBLE PRECISION :: DIST

      LOGICAL :: FREEZE_VEL
      DOUBLE PRECISION :: FREEZE(3)

      DO bcv_i = 1, dem_bcmo

         bcv = dem_bcmo_map(bcv_i)

         freeze_vel = (bc_type_enum(bcv) /= mass_outflow)

         SELECT CASE (bc_plane(bcv))
         CASE('N'); freeze = (/0.0d0, 1.0d0, 0.0d0/); idx=2; sgn=-1.0d0
         CASE('S'); freeze = (/0.0d0, 1.0d0, 0.0d0/); idx=2; sgn= 1.0d0
         CASE('E'); freeze = (/1.0d0, 0.0d0, 0.0d0/); idx=1; sgn=-1.0d0
         CASE('W'); freeze = (/1.0d0, 0.0d0, 0.0d0/); idx=1; sgn= 1.0d0
         CASE('T'); freeze = (/0.0d0, 0.0d0, 1.0d0/); idx=3; sgn=-1.0d0
         CASE('B'); freeze = (/0.0d0, 0.0d0, 1.0d0/); idx=3; sgn= 1.0d0
         END SELECT

         DO lc=dem_bcmo_ijkstart(bcv_i), dem_bcmo_ijkend(bcv_i)
            ijk = dem_bcmo_ijk(lc)
            DO lp= 1,dg_pic(ijk)%ISIZE

               np = dg_pic(ijk)%P(lp)

               IF(is_nonexistent(np)) cycle
               IF(is_any_ghost(np)) cycle
               IF(is_entering(np)) cycle

               SELECT CASE (bc_plane(bcv))
               CASE('S'); dist = yn(bc_j_s(bcv)-1) - des_pos_new(np,2)
               CASE('N'); dist = des_pos_new(np,2) - yn(bc_j_s(bcv))
               CASE('W'); dist = xe(bc_i_w(bcv)-1) - des_pos_new(np,1)
               CASE('E'); dist = des_pos_new(np,1) - xe(bc_i_w(bcv))
               CASE('B'); dist = zt(bc_k_b(bcv)-1) - des_pos_new(np,3)
               CASE('T'); dist = des_pos_new(np,3) - zt(bc_k_b(bcv))
               END SELECT

! The particle is still inside the domain
               IF(dist > des_radius(np)) THEN
                  CALL set_normal(np)

! Check if the particle is crossing over the outlet plane.
               ELSEIF(dist > zero) THEN

! The velocity is 'frozen' normal to the outflow plane. This approach
! is strict because complex BCs (via STLs) can let particles pop through
! the wall along the outlet.
                  IF(freeze_vel) THEN
! Only 'freeze' a particle's velocy if it has it moving out of the
! domain. Otherwise, particles flagged as exiting but moving away from
! the BC appear to moon-walk through the domain until it crashes.
                     IF(des_vel_new(np,idx)*sgn > 0.0d0) THEN
                        des_vel_new(np,:) = des_vel_new(np,:)*freeze(:)
! Set the flags for an exiting particle.
                        IF (is_ghost(np)) THEN
                           CALL set_exiting_ghost(np)
                        ELSE
                           CALL set_exiting(np)
                        ENDIF
                     ENDIF

! The user specified velocity is applied to the exiting particle. This
! only applies to mass outflows where the speed at which particles
! exit needs to be controled.
                  ELSE
                     m = pijk(np,5)
                     des_vel_new(np,1) = bc_u_s(bcv,m)
                     des_vel_new(np,2) = bc_v_s(bcv,m)
                     des_vel_new(np,3) = bc_w_s(bcv,m)
! Set the flags for an exiting particle.
                     IF (is_ghost(np)) THEN
                        CALL set_exiting_ghost(np)
                     ELSE
                        CALL set_exiting(np)
                     ENDIF
                  ENDIF

! Ladies and gentlemen, the particle has left the building.
               ELSE
                  CALL delete_particle(np)
                  force_nsearch = .true.
               ENDIF

            ENDDO
         ENDDO
      ENDDO

! Sync the search flag across all processes.
!      CALL GLOBAL_ALL_OR(FORCE_NSEARCH)

      RETURN
      END SUBROUTINE mass_outflow_dem

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: DELETE_PARTICLE                                         !
!  Author: J.Musser                                   Date: 13-Jul-09  !
!                                                                      !
!  Purpose:  This routine is used to check if a new particle has fully !
!  entered the domain.  If so, the flag classifying the particle as new!
!  is removed, allowing the particle to respond to contact forces from !
!  walls and other particles.                                          !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE delete_particle(NP)

      USE compar
      USE constant
      USE des_bc
      USE discretelement
      USE funits
      USE geometry
      USE indices
      USE param1
      USE physprop
      USE functions

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: NP

!-----------------------------------------------
! Local variables
!-----------------------------------------------

      iglobal_id(np) = -1
      CALL set_nonexistent(np)

      des_pos_new(np,:) = zero
      des_vel_new(np,:) = zero
      omega_new(np,:) = zero

      IF(particle_orientation) orientation(1:3,np) = init_orientation

      IF (do_old) THEN
         des_pos_old(np,:) = zero
         des_vel_old(np,:) = zero
         omega_old(np,:) = zero
      ENDIF

      des_radius(np) = zero
      pmass(np) = huge(0.0)
      pvol(np) = huge(0.0)
      ro_sol(np) = zero
      omoi(np) = zero

      fc(np,:) = zero
      tow(np,:) = zero

      ppos(np,:) = zero

      wall_collision_facet_id(:,np) = -1

      pip = pip - 1

      RETURN
      END SUBROUTINE delete_particle