calc_collision_wall_mod.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: CALC_COLLISION_WALL                                    C
!  Author: Rahul Garg                               Date: 1-Dec-2013   C
!                                                                      C
!  Purpose: subroutines for particle-wall collisions when cutcell is   C
!           used. Also contains rehack of routines for cfslide and     C
!           cffctow which might be different from the stand alone      C
!           routines. Eventually all the DEM routines will be          C
!           consolidated.                                              C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      MODULE calc_collision_wall

      PRIVATE
      PUBLIC :: calc_dem_force_with_wall_stl,&
      & calc_dem_thermo_with_wall_stl

      CONTAINS

!VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV!
!                                                                      !
!                                                                      !
!                                                                      !
!                                                                      !
!                                                                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_dem_force_with_wall_stl

      USE run
      USE param1
      USE desgrid
      USE discretelement
      USE geometry
      USE compar
      USE constant
      USE indices
      USE stl
      USE stl_functions_des
      USE functions
      Implicit none

      INTEGER :: LL
      INTEGER :: IJK, NF
      DOUBLE PRECISION ::OVERLAP_N, SQRT_OVERLAP

      DOUBLE PRECISION :: V_REL_TRANS_NORM, DISTSQ, RADSQ, CLOSEST_PT(dimn)
! local normal and tangential forces
      DOUBLE PRECISION :: NORMAL(dimn), VREL_T(dimn), DIST(dimn), DISTMOD
      DOUBLE PRECISION, DIMENSION(DIMN) :: FTAN, FNORM, OVERLAP_T

      LOGICAL :: DES_LOC_DEBUG
      INTEGER :: CELL_ID, cell_count
      INTEGER :: PHASELL

      DOUBLE PRECISION :: TANGENT(dimn)
      DOUBLE PRECISION :: FTMD, FNMD
! local values used spring constants and damping coefficients
      DOUBLE PRECISION ETAN_DES_W, ETAT_DES_W, KN_DES_W, KT_DES_W

      double precision :: MAG_OVERLAP_T

      double precision :: line_t
! flag to tell if the orthogonal projection of sphere center to
! extended plane detects an overlap

      DOUBLE PRECISION :: MAX_DISTSQ, DISTAPART, FORCE_COH, R_LM
      INTEGER :: MAX_NF, axis
      DOUBLE PRECISION, DIMENSION(3) :: PARTICLE_MIN, PARTICLE_MAX, POS_TMP

      des_loc_debug = .false. ;      debug_des = .false.
      focus_particle = -1

!$omp parallel default(none) private(LL,ijk,MAG_OVERLAP_T,             &
!$omp    cell_id,radsq,particle_max,particle_min,tangent,              &
!$omp    axis,nf,closest_pt,dist,r_lm,distapart,force_coh,distsq,      &
!$omp    line_t,max_distsq,max_nf,normal,distmod,overlap_n,VREL_T,     &
!$omp    v_rel_trans_norm,phaseLL,sqrt_overlap,kn_des_w,kt_des_w,      &
!$omp    etan_des_w,etat_des_w,fnorm,overlap_t,ftan,ftmd,fnmd,pos_tmp) &
!$omp shared(max_pip,focus_particle,debug_des,                         &
!$omp    pijk,dg_pijk,i_of,j_of,k_of,des_pos_new,    &
!$omp    des_radius,facets_at_dg,vertex,  &
!$omp    hert_kwn,hert_kwt,kn_w,kt_w,des_coll_model_enum,mew_w,tow,    &
!$omp    des_etan_wall,des_etat_wall,dtsolid,fc,norm_face,             &
!$omp    wall_collision_facet_id,wall_collision_PFT,use_cohesion,      &
!$omp    van_der_waals,wall_hamaker_constant,wall_vdw_outer_cutoff,    &
!$omp    wall_vdw_inner_cutoff,asperities,surface_energy)
!$omp do
      DO ll = 1, max_pip

         IF(ll.EQ.focus_particle) debug_des = .true.

! skipping non-existent particles or ghost particles
! make sure the particle is not classified as a new 'entering' particle
! or is already marked as a potential exiting particle
         IF(.NOT.is_normal(ll)) cycle

         cell_id = dg_pijk(ll)

! If no neighboring facet in the surrounding 27 cells, then exit
         IF(facets_at_dg(cell_id)%COUNT < 1) THEN
            wall_collision_facet_id(:,ll) = -1
            wall_collision_pft(:,:,ll) = 0.0d0
            cycle
         ENDIF

! Check particle LL for wall contacts
         radsq = des_radius(ll)*des_radius(ll)

         particle_max(:) = des_pos_new( ll,:) + des_radius(ll)
         particle_min(:) = des_pos_new( ll,:) - des_radius(ll)

         DO cell_count = 1, facets_at_dg(cell_id)%count

            axis = facets_at_dg(cell_id)%dir(cell_count)

            nf = facets_at_dg(cell_id)%id(cell_count)

! Compute particle-particle VDW cohesive short-range forces
            IF(use_cohesion .AND. van_der_waals) THEN

               CALL closestptpointtriangle(des_pos_new(ll,:),          &
                  vertex(:,:,nf), closest_pt(:))

               dist(:) = closest_pt(:) - des_pos_new(ll,:)
               distsq = dot_product(dist, dist)
               r_lm = 2*des_radius(ll)

               IF(distsq < (r_lm+wall_vdw_outer_cutoff)**2) THEN
                  IF(distsq > (wall_vdw_inner_cutoff+r_lm)**2) THEN
                     distapart = (sqrt(distsq)-r_lm)
                     force_coh = wall_hamaker_constant*des_radius(ll) /&
                        (12d0*distapart**2)*(asperities/(asperities +  &
                        des_radius(ll)) + one/(one+asperities/         &
                        distapart)**2)
                  ELSE
                     force_coh = 2d0*pi*surface_energy*des_radius(ll)* &
                        (asperities/(asperities+des_radius(ll)) + one/ &
                        (one+asperities/wall_vdw_inner_cutoff)**2 )
                  ENDIF
                  fc(ll,:) = fc(ll,:) + dist(:)*force_coh/sqrt(distsq)
               ENDIF
            ENDIF

            if (facets_at_dg(cell_id)%min(cell_count) >    &
               particle_max(axis)) then
               call remove_collision(ll, nf, wall_collision_facet_id)
               cycle
            endif

            if (facets_at_dg(cell_id)%max(cell_count) <    &
               particle_min(axis)) then
               call remove_collision(ll, nf, wall_collision_facet_id)
               cycle
            endif

! Checking all the facets is time consuming due to the expensive
! separating axis test. Remove this facet from contention based on
! a simple orthogonal projection test.

! Parametrize a line as p = p_0 + t normal and intersect with the
! triangular plane. If t>0, then point is on the non-fluid side of
! the plane. If the plane normal is assumed to point toward the fluid.

! -undefined, because non zero values will imply the sphere center
! is on the non-fluid side of the plane. Since the testing
! is with extended plane, this could very well happen even
! when the particle is well inside the domain (assuming the plane
! normal points toward the fluid). See the pic below. So check
! only when line_t is negative

!                 \   Solid  /
!                  \  Side  /
!                   \      /
!                    \    /
! Wall 1, fluid side  \  /  Wall 2, fluid side
!                      \/
!                        o particle
!
! line_t will be positive for wall 1 (incorrectly indicating center
! is outside the domain) and line_t will be negative for wall 2.
!
! Therefore, only stick with this test when line_t is negative and let
! the separating axis test take care of the other cases.

! Since this is for checking static config, line's direction is the
! same as plane's normal. For moving particles, the line's normal will
! be along the point joining new and old positions.

            line_t = dot_product(vertex(1,:,nf) - des_pos_new(ll,:),&
               norm_face(:,nf))

! k - rad >= tol_orth, where k = -line_t, then orthogonal
! projection is false. Substituting for k
! => line_t + rad <= -tol_orth
! choosing tol_orth = 0.01% of des_radius = 0.0001*des_radius

! Orthogonal projection will detect false positives even
! when the particle does not overlap the triangle.
! However, if the orthogonal projection shows no overlap, then
! that is a big fat negative and overlaps are not possible.
            if((line_t.le.-1.0001d0*des_radius(ll))) then  ! no overlap
               call remove_collision(ll,nf,wall_collision_facet_id)
               cycle
            ENDIF

            pos_tmp = des_pos_new(ll,:)
            CALL closestptpointtriangle(pos_tmp,             &
               vertex(:,:,nf), closest_pt(:))

            dist(:) = closest_pt(:) - des_pos_new(ll,:)
            distsq = dot_product(dist, dist)

            IF(distsq .GE. radsq - small_number) THEN !No overlap exists
               call remove_collision(ll,nf,wall_collision_facet_id)
               cycle
            ENDIF

            max_distsq = distsq
            max_nf = nf

! Assign the collision normal based on the facet with the
! largest overlap.
            normal(:) = dist(:)/sqrt(distsq)

! Facet's normal is correct normal only when the intersection is with
! the face. When the intersection is with edge or vertex, then the
! normal is based on closest pt and sphere center. The definition above
! of the normal is generic enough to account for differences between
! vertex, edge, and facet.

! Calculate the particle/wall overlap.
            distmod = sqrt(max_distsq)
            overlap_n = des_radius(ll) - distmod

! Calculate the translational relative velocity
            CALL cfrelvel_wall(ll, v_rel_trans_norm,vrel_t,           &
               normal, distmod)

! Calculate the spring model parameters.
            phasell = pijk(ll,5)

! Hertz vs linear spring-dashpot contact model
            IF (des_coll_model_enum .EQ. hertzian) THEN
               sqrt_overlap = sqrt(overlap_n)
               kn_des_w = hert_kwn(phasell)*sqrt_overlap
               kt_des_w = hert_kwt(phasell)*sqrt_overlap
               sqrt_overlap = sqrt(sqrt_overlap)
               etan_des_w = des_etan_wall(phasell)*sqrt_overlap
               etat_des_w = des_etat_wall(phasell)*sqrt_overlap
            ELSE
               kn_des_w = kn_w
               kt_des_w = kt_w
               etan_des_w = des_etan_wall(phasell)
               etat_des_w = des_etat_wall(phasell)
            ENDIF

! Calculate the normal contact force
            fnorm(:) = -(kn_des_w * overlap_n * normal(:) + &
               etan_des_w * v_rel_trans_norm * normal(:))

! Calculate the tangential displacement.
            overlap_t(:) = dtsolid*vrel_t(:) + get_collision(ll,       &
               nf, wall_collision_facet_id, wall_collision_pft)
            mag_overlap_t = sqrt(dot_product(overlap_t, overlap_t))

! Check for Coulombs friction law and limit the maximum value of the
! tangential force on a particle in contact with a wall.
            IF(mag_overlap_t > 0.0) THEN
! Tangential froce from spring.
               ftmd = kt_des_w*mag_overlap_t
! Max force before the on set of frictional slip.
               fnmd = mew_w*sqrt(dot_product(fnorm,fnorm))
! Direction of tangential force.
               tangent = overlap_t/mag_overlap_t
               IF(ftmd < fnmd) THEN
                  ftan = -ftmd * tangent
               ELSE
                  ftan = -fnmd * tangent
                  overlap_t = (fnmd/kt_des_w) * tangent
               ENDIF
            ELSE
               ftan = 0.0
            ENDIF
! Add in the tangential dashpot damping force
            ftan = ftan - etat_des_w*vrel_t(:)

! Save the tangential displacement.
            CALL update_collision(overlap_t, ll, nf,                   &
               wall_collision_facet_id, wall_collision_pft)

! Add the collision force to the total forces acting on the particle.
            fc(ll,:) = fc(ll,:) + fnorm(:) + ftan(:)

! Add the torque force to the toal torque acting on the particle.
            tow(ll,:) = tow(ll,:) + distmod*des_crossprdct(normal,ftan)

         ENDDO

      ENDDO
!$omp end do
!$omp end parallel

      RETURN

       contains

!......................................................................!
!  Function: GET_COLLISION                                             !
!                                                                      !
!  Purpose: Return the integrated (t0->t) tangential displacement.     !
!......................................................................!
      FUNCTION get_collision(LLL,FACET_ID,WALL_COLLISION_FACET_ID,     &
          wall_collision_pft)

      use stl_dbg_des, only: write_this_stl
      use stl_dbg_des, only: write_stls_this_dg

      use error_manager

      IMPLICIT NONE

      DOUBLE PRECISION :: GET_COLLISION(dimn)
      INTEGER, INTENT(IN) :: LLL,FACET_ID
      INTEGER, INTENT(INOUT) :: WALL_COLLISION_FACET_ID(:,:)
      DOUBLE PRECISION, INTENT(INOUT) :: WALL_COLLISION_PFT(:,:,:)
      INTEGER :: CC, FREE_INDEX, LC, dgIJK

      INTEGER :: lSIZE1, lSIZE2, lSIZE3
      INTEGER, ALLOCATABLE :: tmpI2(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: tmpR3(:,:,:)

      free_index = -1

      do cc = 1, collision_array_max
         if (facet_id == wall_collision_facet_id(cc,lll)) then
            get_collision(:) = wall_collision_pft(:,cc,lll)
            return
         else if (-1 == wall_collision_facet_id(cc,lll)) then
            free_index = cc
         endif
      enddo

! Overwrite old data. This is needed because a particle moving from
! one dg cell to another may no longer 'see' an STL before it moved
! out of contact range. Therefore, the 'remove_collision' function
! does not get called to cleanup the stale data.
      if(-1 == free_index) then
         dgijk=dg_pijk(lll)
         cc_lp: do cc=1, collision_array_max
            do lc=1, facets_at_dg(dgijk)%count
               if(wall_collision_facet_id(cc,lll) == &
                  facets_at_dg(dgijk)%id(lc))  cycle cc_lp
            enddo
            free_index = cc
            exit cc_lp
         enddo cc_lp
      endif

! Last resort... grow the collision array
      if(-1 == free_index) then
         free_index=collision_array_max+1
         collision_array_max = 2*collision_array_max
         CALL grow_wall_collision(collision_array_max)
      endif

      wall_collision_facet_id(free_index,lll) = facet_id
      wall_collision_pft(:,free_index,lll) = zero
      get_collision(:) = wall_collision_pft(:,free_index,lll)
      return

      END FUNCTION get_collision


!......................................................................!
!  Subroutine: GROW_WALL_COLLISION                                     !
!                                                                      !
!  Purpose: Return the integrated (t0->t) tangential displacement.     !
!......................................................................!
      SUBROUTINE grow_wall_collision(NEW_SIZE)

      use discretelement

      use stl_dbg_des, only: write_this_stl
      use stl_dbg_des, only: write_stls_this_dg

      use error_manager

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: NEW_SIZE
      INTEGER :: lSIZE1, lSIZE2, lSIZE3
      INTEGER, ALLOCATABLE :: tmpI2(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: tmpR3(:,:,:)

      lsize1 = size(wall_collision_facet_id,1)
      lsize2 = size(wall_collision_facet_id,2)

      allocate(tmpi2(new_size, lsize2))
      tmpi2(1:lsize1,:) = wall_collision_facet_id(1:lsize1,:)
      call move_alloc(tmpi2, wall_collision_facet_id)
      wall_collision_facet_id(lsize1+1:new_size,:) = -1

      lsize1 = size(wall_collision_pft,1)
      lsize2 = size(wall_collision_pft,2)
      lsize3 = size(wall_collision_pft,3)

      allocate(tmpr3(lsize1, new_size, lsize3))
      tmpr3(:,1:lsize2,:) = wall_collision_pft(:,1:lsize2,:)
      call move_alloc(tmpr3, wall_collision_pft)

      RETURN
      END SUBROUTINE grow_wall_collision




!......................................................................!
!  Function: UPDATE_COLLISION                                          !
!                                                                      !
!  Purpose: Update the integrated (t0->t) tangential displacement.     !
!......................................................................!
      SUBROUTINE update_collision(PFT, LLL, FACET_ID,                  &
         wall_collision_facet_id, wall_collision_pft)

      use error_manager
      implicit none

      DOUBLE PRECISION, INTENT(IN) :: PFT(dimn)
      INTEGER, INTENT(IN) :: LLL,FACET_ID
      INTEGER, INTENT(IN) :: WALL_COLLISION_FACET_ID(:,:)
      DOUBLE PRECISION, INTENT(INOUT) :: WALL_COLLISION_PFT(:,:,:)
      INTEGER :: CC

      do cc = 1, collision_array_max
         if (facet_id == wall_collision_facet_id(cc,lll)) then
            wall_collision_pft(:,cc,lll) = pft(:)
            return
         endif
      enddo

      CALL init_err_msg("CALC_COLLISION_WALL_MOD: UPDATE_COLLISION")
      WRITE(err_msg, 1100)
      CALL flush_err_msg(abort=.true.)

 1100 FORMAT('Error: COLLISION_ARRAY_MAX too small. ')

      END SUBROUTINE update_collision

!......................................................................!
!  Function: REMOVE_COLLISION                                          !
!                                                                      !
!  Purpose: Clear the integrated (t0->t) tangential displacement once  !
!  the collision is over (contact ended).                              !
!......................................................................!
      SUBROUTINE remove_collision(LLL,FACET_ID,WALL_COLLISION_FACET_ID)

      use error_manager

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: LLL,FACET_ID
      INTEGER, INTENT(INOUT) :: WALL_COLLISION_FACET_ID(:,:)
      INTEGER :: CC

      DO cc = 1, collision_array_max
         IF (facet_id == wall_collision_facet_id(cc,lll)) THEN
            wall_collision_facet_id(cc,lll) = -1
            RETURN
         ENDIF
      ENDDO

      END SUBROUTINE remove_collision

      END SUBROUTINE calc_dem_force_with_wall_stl


!----------------------------------------------------------------------!
!                                                                      !
!  Subroutine: CFRELVEL_WALL                                           !
!                                                                      !
!  Purpose: Calculate the normal and tangential components of the      !
!  relative velocity between a particle and wall contact.              !
!                                                                      !
!  Comments: Only the magnitude of the normal component is returned    !
!  whereas the full tangential vector is returned.                     !
!----------------------------------------------------------------------!
      SUBROUTINE cfrelvel_wall(LL, VRN, VRT, NORM, DIST)

! Particle translational velocity
      use discretelement, only: des_vel_new
! Particle rotational velocity
      use discretelement, only: omega_new
! Spatial array size (parameter)
      use discretelement, only: dimn
! Function for calculating the cross prodcut
      use discretelement, only: des_crossprdct

      IMPLICIT NONE

! Dummy arguments:
!---------------------------------------------------------------------//
! Particle index.
      INTEGER, INTENT(IN) :: LL
! Magnitude of the total relative translational velocity.
      DOUBLE PRECISION, INTENT(OUT):: VRN
! Total relative translational velocity (vector).
      DOUBLE PRECISION, DIMENSION(DIMN), INTENT(OUT):: VRT
! Unit normal from particle center to closest point on stl (wall)
      DOUBLE PRECISION, DIMENSION(DIMN), INTENT(IN) :: NORM
! Distance between particle center and stl (wall).
      DOUBLE PRECISION, INTENT(IN) :: DIST

! Local variables
!---------------------------------------------------------------------//
! Additional relative translational motion due to rotation
      DOUBLE PRECISION, DIMENSION(DIMN) :: V_ROT
! Total relative velocity at contact point
      DOUBLE PRECISION, DIMENSION(DIMN) :: VRELTRANS

! Total relative velocity + rotational contribution
      v_rot = dist*omega_new(ll,:)
      vreltrans(:) =  des_vel_new(ll,:) + des_crossprdct(v_rot, norm)

! magnitude of normal component of relative velocity (scalar)
      vrn = dot_product(vreltrans,norm)

! total relative translational slip velocity at the contact point
! Equation (8) in Tsuji et al. 1992
      vrt(:) =  vreltrans(:) - vrn*norm(:)

      RETURN
      END SUBROUTINE cfrelvel_wall


!----------------------------------------------------------------//
! SUBROUTINE: CALC_DEM_THERMO_WITH_WALL_STL
! By: Aaron M.
! Purpose: Compute heat transfer to particles due to boundaries
!----------------------------------------------------------------//

      SUBROUTINE calc_dem_thermo_with_wall_stl

      USE bc
      USE compar, only: istart3, iend3, jstart3, jend3, kstart3, kend3
      USE cutcell, only: cut_cell_at, area_cut, blocked_cell_at, small_cell_at
      USE des_thermo
      USE des_thermo_cond
      USE desgrid
      USE discretelement
      USE exit, only: mfix_exit
      USE functions
      USE geometry, only: do_k, imax1, imax2, imin1, imin2
      USE geometry, only: no_k, jmax1, jmax2, jmin1, jmin2
      USE geometry, only: kmax1, kmax2, kmin1, kmin2, zlength
      USE geometry, only: dx, dy, dz
      USE param, only: dimension_i, dimension_j, dimension_k
      USE param1
      USE physprop
      USE run, only: units, time
      USE stl
      USE stl_functions_des

      IMPLICIT NONE

      INTEGER :: LL             ! Loop index for particle
      INTEGER :: CELL_ID        ! Desgrid cell index
      INTEGER :: CELL_COUNT     ! Loop index for facets
      INTEGER :: IJK_FLUID      ! IJK index for fluid cell
      INTEGER :: I_FLUID, J_FLUID, K_FLUID
      INTEGER :: I_Facet, J_Facet, K_Facet, IJK_Facet
      INTEGER :: I1,J1,K1
      INTEGER :: phase_LL       ! Phase index for particle

      INTEGER, PARAMETER :: MAX_CONTACTS = 12
      INTEGER :: iFacet         ! loop index for facets
      INTEGER :: count_Facets   ! counter for number of facets particle contacts
      DOUBLE PRECISION, DIMENSION(3,MAX_CONTACTS) :: NORM_FAC_CONTACT
      LOGICAL :: USE_FACET
      LOGICAL :: DOMAIN_BDRY

      INTEGER :: NF             ! Facet ID
      INTEGER :: AXIS           ! Facet direction
      DOUBLE PRECISION :: RLENS_SQ ! lens radius squared
      DOUBLE PRECISION :: RLENS ! lens radius
      DOUBLE PRECISION :: RPART ! Particle radius

      ! vectors for minimum / maximum possible particle contacts
      DOUBLE PRECISION, DIMENSION(3) :: PARTPOS_MIN, PARTPOS_MAX
      DOUBLE PRECISION, DIMENSION(3) :: POS_TMP ! Temp vector
      DOUBLE PRECISION, DIMENSION(3) :: DIST, NORMAL

      DOUBLE PRECISION, DIMENSION(3) :: CLOSEST_PT

      DOUBLE PRECISION :: line_t  ! Normal projection for part/wall
      DOUBLE PRECISION :: DISTSQ ! Separation from particle and facet (squared)
      DOUBLE PRECISION :: PROJ

      INTEGER :: BC_ID ! BC ID
      INTEGER :: IBC   ! BC Loop Index

      DOUBLE PRECISION :: TWALL
      DOUBLE PRECISION :: K_gas
      DOUBLE PRECISION :: TPART
      DOUBLE PRECISION :: OVERLAP
      DOUBLE PRECISION :: QSWALL, AREA

      LOGICAL, SAVE :: OUTPUT_WARNING = .true.



    !  IF(.NOT.DES_CONTINUUM_COUPLED.or.DES_EXPLICITLY_COUPLED)THEN
    !     CALL PARTICLES_IN_CELL
    !  ENDIF
      DO ll = 1, max_pip
         ! Skip non-existent particles or ghost particles
         IF (.NOT.is_normal(ll)) cycle
         phase_ll = pijk(ll,5)
         IF(.NOT.calc_cond_des(phase_ll))cycle
         ! Get desgrid cell index
         cell_id = dg_pijk(ll)

         ! Skip cells that do not have neighboring facet
         IF (facets_at_dg(cell_id)%COUNT <1) cycle

         ! Store lens radius of particle
         rlens = (one+flpc)*des_radius(ll)
         rlens_sq = rlens*rlens

         rpart = des_radius(ll)

         ! Compute max/min particle locations
         partpos_max(:) = des_pos_new(ll,:) + rlens
         partpos_min(:) = des_pos_new(ll,:) - rlens

         ! Get fluid cell
         i_fluid = pijk(ll,1)
         j_fluid = pijk(ll,2)
         k_fluid = pijk(ll,3)
         ijk_fluid = pijk(ll,4)
         tpart = des_t_s(ll)

         ! Sometimes PIJK is not updated every solids timestep and
         ! therefore doesn't give the correct fluid cell that
         ! contains the particles.  If so, determine actual fluid
         ! cell that contains the particle
         IF(.NOT.des_continuum_coupled.or.des_explicitly_coupled)THEN
            IF(i_fluid <= istart3 .OR. i_fluid >= iend3) THEN
               CALL pic_search(i_fluid, des_pos_new(ll,1), xe,       &
               dimension_i, imin2, imax2)
            ELSE
               IF((des_pos_new(ll,1) >= xe(i_fluid-1)) .AND.         &
               (des_pos_new(ll,1) <  xe(i_fluid))) THEN
               i_fluid = i_fluid
               ELSEIF((des_pos_new(ll,1) >= xe(i_fluid)) .AND.       &
                  (des_pos_new(ll,1) < xe(i_fluid+1))) THEN
                  i_fluid = i_fluid+1
               ELSEIF((des_pos_new(ll,1) >= xe(i_fluid-2)) .AND.     &
                  (des_pos_new(ll,1) < xe(i_fluid-1))) THEN
                  i_fluid = i_fluid-1
               ELSE
                  CALL pic_search(i_fluid, des_pos_new(ll,1), xe,    &
                  dimension_i, imin2, imax2)
               ENDIF
            ENDIF !(I)
            IF(j_fluid <= jstart3 .OR. j_fluid >= jend3) THEN
               CALL pic_search(j_fluid, des_pos_new(ll,2), yn,       &
               dimension_j, jmin2, jmax2)
            ELSE
               IF((des_pos_new(ll,2) >= yn(j_fluid-1)) .AND.         &
                  (des_pos_new(ll,2) <  yn(j_fluid))) THEN
                  j_fluid = j_fluid
               ELSEIF((des_pos_new(ll,2) >= yn(j_fluid)) .AND.       &
                  (des_pos_new(ll,2) < yn(j_fluid+1))) THEN
                  j_fluid = j_fluid+1
               ELSEIF((des_pos_new(ll,2) >= yn(j_fluid-2)) .AND.     &
                  (des_pos_new(ll,2) < yn(j_fluid-1))) THEN
                  j_fluid = j_fluid-1
               ELSE
                  CALL pic_search(j_fluid, des_pos_new(ll,2), yn,    &
                  dimension_j, jmin2, jmax2)
               ENDIF
            ENDIF !(J)

            IF(no_k) THEN
               k_fluid = 1
            ELSE
               IF(k_fluid <= kstart3 .OR. k_fluid >= kend3) THEN
               CALL pic_search(k_fluid, des_pos_new(ll,3), zt,         &
                  dimension_k, kmin2, kmax2)
               ELSE
                  IF((des_pos_new(ll,3) >= zt(k_fluid-1)) .AND.        &
                     (des_pos_new(ll,3) < zt(k_fluid))) THEN
                     k_fluid = k_fluid
                  ELSEIF((des_pos_new(ll,3) >= zt(k_fluid)) .AND.      &
                     (des_pos_new(ll,3) < zt(k_fluid+1))) THEN
                     k_fluid = k_fluid+1
                  ELSEIF((des_pos_new(ll,3) >= zt(k_fluid-2)) .AND.    &
                     (des_pos_new(ll,3) < zt(k_fluid-1))) THEN
                     k_fluid = k_fluid-1
                  ELSE
                     CALL pic_search(k_fluid, des_pos_new(ll,3), zt,   &
                     dimension_k, kmin2, kmax2)
                  ENDIF
               ENDIF
            ENDIF !(K)
            ijk_fluid = pijk(ll,4)
         ENDIF ! (NOT CONTINUUM_COUPLED OR EXPLICIT)


! Initialize counter for number of facets
         count_facets = 0

         ! Loop over potential facets
         DO cell_count = 1, facets_at_dg(cell_id)%count
            ! Get direction (axis) and facet id (NF)
            axis = facets_at_dg(cell_id)%dir(cell_count)
            nf = facets_at_dg(cell_id)%id(cell_count)

            ! Check to see if facet is out of reach of particle
            if (facets_at_dg(cell_id)%min(cell_count) >    &
               partpos_max(axis)) then
               cycle
            endif
            if (facets_at_dg(cell_id)%max(cell_count) <    &
               partpos_min(axis)) then
               cycle
            endif

            ! Compute projection (normal distance) from wall to particle
            line_t = dot_product(vertex(1,:,nf) - des_pos_new(ll,:),&
            &        norm_face(:,nf))

            ! If normal exceeds particle radius, particle is not in contact
            if((line_t.lt.-rlens))cycle

            ! Compute closest point on facet
            pos_tmp(:) = des_pos_new(ll,:)
            CALL closestptpointtriangle(pos_tmp, vertex(:,:,nf), &
            &    closest_pt(:))
            ! Compute position vector from particle to closest point on facet
            dist(:) = closest_pt(:)-pos_tmp(:)
            distsq = dot_product(dist,dist)

            ! Skip particles that are more than lens radius from facet
            IF(distsq .GE. (rlens_sq-small_number))cycle

            ! Only do heat transfer for particles that are directly above facet
            ! Heat transfer routines not generalized for edge contacts, but
            ! those should normally yield negligible heat transfer contributions

            ! Normalize distance vector and compare to facet normal
            normal(:)=-dist(:)/sqrt(distsq)
            proj = sqrt(abs(dot_product(normal, norm_face(:,nf))))
            IF(abs(one-proj).gt.1.0d-6)cycle

            ! Get overlap
            overlap = rpart - sqrt(distsq)

            ! Initialize area for facet (for post-proc. flux)
            area = zero

            ! Initialize BDRY FLAG
            domain_bdry = .false.
            ! Get BC_ID
            bc_id = bc_id_stl_face(nf)

            ! BC_ID is set to 0 in pre-proc if stl is a domain boundary
            IF(bc_id.eq.0)then
               i1=i_fluid
               j1=j_fluid
               k1=k_fluid
               domain_bdry = .true.

               ! Domain boundary, figure out real boundary ID
               IF(norm_face(1,nf).ge.0.9999)THEN
                  ! WEST face
                  i1=imin2
                  IF(do_k)THEN
                     area = dy(j_fluid)*dz(k_fluid)
                  ELSE
                     area = dy(j_fluid)*zlength
                  ENDIF
               ELSEIF(norm_face(1,nf).le.-0.9999)THEN
                  ! EAST FACE
                  i1=imax2
                  IF(do_k)THEN
                     area = dy(j_fluid)*dz(k_fluid)
                  ELSE
                     area = dy(j_fluid)*zlength
                  ENDIF
               ELSEIF(norm_face(2,nf).ge.0.9999)THEN
                  ! SOUTH FACE
                  j1=jmin2
                  IF(do_k)THEN
                     area = dx(i_fluid)*dz(k_fluid)
                  ELSE
                     area = dx(i_fluid)*zlength
                  ENDIF
               ELSEIF(norm_face(2,nf).le.-0.9999)THEN
                  ! NORTH FACE
                  j1=jmax2
                  IF(do_k)THEN
                     area = dx(i_fluid)*dz(k_fluid)
                  ELSE
                     area = dx(i_fluid)*zlength
                  ENDIF
               ELSEIF(norm_face(3,nf).ge.0.9999)THEN
                  ! BOTTOM FACE
                  k1=kmin2
                  area = dx(i_fluid)*dy(j_fluid)

               ELSEIF(norm_face(3,nf).le.-0.9999)THEN
                  ! TOP FACE
                  k1=kmax2
                  area = dx(i_fluid)*dy(j_fluid)
               ELSE
                  WRITE( *,*)'PROBLEM, COULD NOT FIND DOMAIN BOUNDARY'
                  WRITE(*,*)' In calc_thermo_des_wall_stl'
                  call mfix_exit(1)

               ENDIF


               ! Loop through defined BCs to see which one particle neighbors
               DO ibc = 1, dimension_bc
                  IF(.NOT.bc_defined(ibc))cycle
                  IF (i1.ge.bc_i_w(ibc).and.i1.le.bc_i_e(ibc).and.&
                      j1.ge.bc_j_s(ibc).and.j1.le.bc_j_n(ibc).and.&
                      k1.ge.bc_k_b(ibc).and.k1.le.bc_k_t(ibc))THEN
                      bc_id = ibc
                      exit
                  ENDIF
               ENDDO
               IF(bc_id.eq.0)then
                  IF(output_warning)THEN
                     write(*,*)'WARNING: Could not find BC'
                     write(*,*)'Check input deck to make sure domain boundaries &
                     are defined'
                     write(*,*)'SUPPRESSING FUTURE OUTPUT'
                     write(*,*)'DES_POS_NEW'
                     write(*,'(3(F12.5, 3X))')(des_pos_new(ll,ibc),ibc=1,3)
                     write(*,*)'I,J,K'
                     write(*,*)i1,j1,k1
                     write(*,*)'CLOSEST PT'
                     write(*,'(3(F12.5, 3X))')(closest_pt(ibc),ibc=1,3)
                     write(*,*)'NORM_FACE'
                     write(*,'(3(F12.5, 3X))')(norm_face(ibc,nf),ibc=1,3)
                     output_warning = .false.
                  ENDIF
                  cycle  !
               ENDIF

            ENDIF !Domain Boundary (facet ID was 0)

            IF (bc_type_enum(bc_id) == no_slip_wall .OR. &
               bc_type_enum(bc_id) == free_slip_wall .OR. &
               bc_type_enum(bc_id) == par_slip_wall .OR. &
               bc_type_enum(bc_id) == cg_nsw .OR. &
               bc_type_enum(bc_id) == cg_fsw .OR. &
               bc_type_enum(bc_id) == cg_psw) THEN



               ! CHECK TO MAKE SURE FACET IS UNIQUE
               use_facet=.true.
               DO ifacet=1,count_facets
                  ! DO CHECK BY ENSURING NORMAL VECTOR IS NEARLY PARALLEL
                  proj = sqrt(abs(dot_product(normal, norm_fac_contact(:,ifacet))))
                  IF(abs(one-proj).gt.1.0d-6)THEN
                     use_facet=.false.
                     EXIT
                  ENDIF
               ENDDO
               IF(.NOT.use_facet)cycle

               ! FACET IS UNIQUE
               count_facets=count_facets+1
               norm_fac_contact(:,count_facets)=normal(:)

! Do heat transfer
               ! GET WALL TEMPERATURE
               twall = bc_tw_s(bc_id,phase_ll)

               ! GET GAS THERMAL CONDUCTIVITY
               if(k_g0.eq.undefined)then
                  ! Compute gas conductivity as is done in calc_k_g
                  ! But use average of particle and wall temperature to be gas temperature

                  k_gas = 6.02d-5*sqrt(half*(twall+tpart)/300.d0) ! cal/(s.cm.K)
                  ! 1 cal = 4.183925D0 J
                  IF (units == 'SI') k_gas = 418.3925d0*k_gas !J/s.m.K
               else
                  k_gas=k_g0
               endif
               IF(twall.eq.undefined)cycle
               qswall = des_conduction_wall(overlap,k_s0(phase_ll), &
               &        k_s0(phase_ll),k_gas,twall, tpart, rpart, &
               &        rlens, phase_ll)


               q_source(ll) = q_source(ll)+qswall

               ! BELOW CODE IS ONLY NECESSARY FOR OUTPUTING
               ! DATA.  Need to know fluid cell that contact
               ! point resides in so that wall flux can be
               ! output correctly.

               i_facet = i_fluid
               j_facet = j_fluid
               k_facet = k_fluid

               ! This checks to see if the contact point was NOT on a domain boundary

               IF(.NOT.domain_bdry)THEN
                  IF(closest_pt(1) >= xe(i_fluid))THEN
                     i_facet = min(imax1, i_fluid+1)
                  ELSEIF(closest_pt(1) < xe(i_fluid-1))THEN
                     i_facet = max(imin1, i_fluid-1)
                  ENDIF

                  IF(closest_pt(2) >= yn(j_fluid))THEN
                     j_facet = min(jmax1, j_fluid+1)
                  ELSEIF(closest_pt(2) < yn(j_fluid-1))THEN
                     j_facet = max(jmin1, j_fluid-1)
                  ENDIF
                  IF(do_k)THEN
                     IF(closest_pt(3) >= zt(k_fluid))THEN
                        k_facet = min(kmax1, k_fluid+1)
                     ELSEIF(closest_pt(3) < zt(k_fluid-1))THEN
                        k_facet = max(kmin1, k_fluid-1)
                     ENDIF
                  ENDIF
                  ijk_facet=funijk(i_facet,j_facet,k_facet)
                  area=area_cut(ijk_facet)

               ! AREA is left undefined if contact was with cut-cell surface
                  IF (area.eq.zero)then
                     i_facet = i_fluid
                     j_facet = j_fluid
                     k_facet = k_fluid
                     ijk_facet=funijk(i_facet,j_facet,k_facet)
                     IF(norm_face(1,nf).ge.0.9999)THEN
                     ! WEST face
                        IF(do_k)THEN
                           area = dy(j_facet)*dz(k_facet)
                        ELSE
                           area = dy(j_facet)*zlength
                        ENDIF
                     ELSEIF(norm_face(1,nf).le.-0.9999)THEN
                     ! EAST FACE
                        IF(do_k)THEN
                           area = dy(j_facet)*dz(k_facet)
                        ELSE
                           area = dy(j_facet)*zlength
                        ENDIF
                     ELSEIF(norm_face(2,nf).ge.0.9999)THEN
                     ! SOUTH FACE
                        IF(do_k)THEN
                           area = dx(i_facet)*dz(k_facet)
                        ELSE
                           area = dx(i_facet)*zlength
                        ENDIF
                     ELSEIF(norm_face(2,nf).le.-0.9999)THEN
                     ! NORTH FACE
                        IF(do_k)THEN
                           area = dx(i_facet)*dz(k_facet)
                        ELSE
                           area = dx(i_facet)*zlength
                        ENDIF
                     ELSEIF(norm_face(3,nf).ge.0.9999)THEN
                     ! BOTTOM FACE
                        area = dx(i_facet)*dy(j_facet)
                     ELSEIF(norm_face(3,nf).le.-0.9999)THEN
                     ! TOP FACE
                        area = dx(i_facet)*dy(j_facet)
                     ENDIF
                  ENDIF ! Area==0 (because cut-cell facet exists in non cut-cell)
               ENDIF ! NOT DOMAIN_BDRY
               ijk_facet = funijk(i_facet,j_facet,k_facet)
               ! AM error check
               IF(.NOT.fluid_at(ijk_facet).AND. &
               &  .NOT.blocked_cell_at(ijk_facet))THEN
                  write(*,*)'ERROR: Cell containing facet is not a fluid &
                  &  cell or a blocked cell'
                  write(*,*)fluid_at(ijk_facet), blocked_cell_at(ijk_facet)
                  write(*,*)'PART POS',(des_pos_new(ll,ibc),ibc=1,3)
                  write(*,*)'FACET NORM',(norm_face(ibc,nf),ibc=1,3)
                  write(*,*)'BC_ID', bc_id
                  write(*,*)'I,J,K (Facet)', i_facet,j_facet,k_facet

                  call mfix_exit(1)
               ENDIF

               IF(fluid_at(ijk_facet))THEN
                  des_qw_cond(ijk_facet,phase_ll) = &
                     des_qw_cond(ijk_facet, phase_ll) + qswall/area
               ENDIF

            ENDIF ! WALL BDRY
         ENDDO                  ! CELL_COUNT (facets)
      ENDDO  ! LL
      RETURN

      END SUBROUTINE calc_dem_thermo_with_wall_stl

      end module calc_collision_wall