apply_wall_bc_pic.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: APPLY_WALL_BC_PIC                                       !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose: Detect collisions between PIC particles and STLs.          !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE apply_wall_bc_pic

! Global Variables:
!---------------------------------------------------------------------//
! Particle postions: Current/Previous
      use discretelement, only: des_pos_new
! Paricle velocities
      use discretelement, only: des_vel_new
! Particle radius
      use discretelement, only: des_radius
! Max number of particles on this process
      use discretelement, only: max_pip
! Map from particle to DES grid cell.
      use discretelement, only: dg_pijk, dtsolid
! Flag indicating that the index cell contains no STLs
      use stl, only: facets_at_dg
! STL Vertices
      use stl, only: vertex
! STL Facet normals
      use stl, only: norm_face
! Minimum velocity to offset gravitational forces
      use pic_bc, only: minvel, minvel_mag, oominvel_mag
! Solids time step size
      use discretelement, only: dtsolid
! Gravitational force vector
      use discretelement, only: grav

! Global Parameters:
!---------------------------------------------------------------------//
      use param1, only: zero, small_number, one, undefined
      use functions, only: is_nonexistent

! Module Procedures:
!---------------------------------------------------------------------//
      use stl_functions_des

      use error_manager

      IMPLICIT NONE

! Local Variables:
!---------------------------------------------------------------------//
! Loop counters.
      INTEGER NF, NP, LC1
! STL normal vector
      DOUBLE PRECISION :: NORM_PLANE(3)
! line is parameterized as p = p_ref + t * dir_line, t is line_param
      DOUBLE PRECISION :: LINE_T

! distance parcel travelled out of domain along the facet normal
      DOUBLE PRECISION :: DIST(3), tPOS(3)

      DOUBLE PRECISION :: DT_RBND, RADSQ
!......................................................................!

      double precision :: vv(3)

! Minimum velocity needed to offset gravity.
      minvel = -dtsolid*grav(:)
      minvel_mag = dot_product(minvel,minvel)
      oominvel_mag = 1.0d0/minvel_mag

      dt_rbnd = 1.01d0*dtsolid

      DO np = 1, max_pip

! Skip non-existent particles
         IF(is_nonexistent(np)) cycle

! If no neighboring facet in the surrounding 27 cells, then exit
         IF(facets_at_dg(dg_pijk(np))%COUNT < 1) cycle

         radsq = des_radius(np)*des_radius(np)

! Loop through the STLs associated with the DES grid cell.
         lc1_lp: DO lc1=1, facets_at_dg(dg_pijk(np))%COUNT

! Get the index of the neighbor facet
            nf = facets_at_dg(dg_pijk(np))%ID(lc1)

            norm_plane = norm_face(:,nf)

!            IF(dot_product(NORM_PLANE, DES_VEL_NEW(NP,:)) > 0.d0) &
!               CYCLE LC1_LP

            vv = vertex(1,:,nf)
            CALL intersectlnplane(des_pos_new(np,:), des_vel_new(np,:),&
                vv, norm_face(:,nf), line_t)

            IF(abs(line_t) <= dt_rbnd) THEN

! Project the parcel onto the facet.
               dist = line_t*des_vel_new(np,:)
               tpos = des_pos_new(np,:) + dist

! Avoid collisions with STLs that are not next to the parcel
               IF(hit_facet(vertex(:,:,nf), tpos)) THEN

! Correct the position of particles found too close to the STL.
                  IF(dot_product(dist,dist) <= radsq .OR. &
                     line_t <= zero) THEN
                     des_pos_new(np,:) = des_pos_new(np,:) + dist(:) + &
                        des_radius(np)*norm_plane
                  ENDIF

! Reflect the parcel.
                  CALL pic_reflect_part(np, norm_plane(:))
               ENDIF
            ENDIF

         ENDDO lc1_lp
      END DO

      RETURN

      contains


!......................................................................!
      LOGICAL FUNCTION hit_facet(VERTS, POINT)

      DOUBLE PRECISION, INTENT(IN) :: VERTS(3,3)
      DOUBLE PRECISION, INTENT(IN) :: POINT(3)

      DOUBLE PRECISION :: V0(3), V1(3), V2(3)
      DOUBLE PRECISION :: d00, d01, d02, d11, d12

      DOUBLE PRECISION :: OoDEN, u, v

      v0 = verts(3,:) - verts(1,:)
      v1 = verts(2,:) - verts(1,:)
      v2 = point - verts(1,:)

      d00 = dot_product(v0,v0)
      d01 = dot_product(v0,v1)
      d02 = dot_product(v0,v2)

      d11 = dot_product(v1,v1)
      d12 = dot_product(v1,v2)

      ooden = 1.0d0/(d00*d11 - d01*d01)
      u = (d11*d02 - d01*d12)*ooden
      v = (d00*d12 - d01*d02)*ooden

      hit_facet = ((u>=0.0d0) .AND. (v>=0.0d0) .AND. (u+v < 1.0d0))

      RETURN
      END FUNCTION hit_facet

      END SUBROUTINE apply_wall_bc_pic


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: PIC_REFLECT_PARTICLE                                    C
!  Purpose:                                                            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE pic_reflect_part(LL, WALL_NORM)

      USE discretelement, only : dimn, des_vel_new
      USE mfix_pic, only : mppic_coeff_en_wall, mppic_coeff_et_wall

! Minimum velocity to offset gravitational forces
      use pic_bc, only: minvel, minvel_mag, oominvel_mag


      use param1, only: small_number

      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
      INTEGER, INTENT(IN) :: LL
      DOUBLE PRECISION, INTENT(IN) ::  WALL_NORM(dimn)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
      !magnitude of pre-collisional normal and tangential velocity components
      DOUBLE PRECISION :: VEL_NORMMAG_APP

      !pre collisional normal and tangential velocity components in vector format
      !APP ==> approach
      DOUBLE PRECISION :: VEL_NORM_APP(dimn), VEL_TANG_APP(dimn)


      !post collisional normal and tangential velocity components in vector format
      !SEP ==> separation
      DOUBLE PRECISION :: VEL_NORM_SEP(dimn), VEL_TANG_SEP(dimn)

      DOUBLE PRECISION :: COEFF_REST_EN, COEFF_REST_ET

! Minimum particle velocity with respect to gravity. This ensures that
! walls "push back" to offset gravitational forces.
      DOUBLE PRECISION :: projGRAV(3)

      INTEGER :: LC

!-----------------------------------------------





      vel_normmag_app = dot_product(wall_norm(:), des_vel_new(ll,:))

! Currently assuming that wall is at rest. Needs improvement for moving wall

      vel_norm_app(:) = vel_normmag_app*wall_norm(:)
      vel_tang_app(:) = des_vel_new(ll,:) - vel_norm_app(:)


!post collisional velocities

      coeff_rest_en = mppic_coeff_en_wall
      coeff_rest_et = mppic_coeff_et_wall

      vel_norm_sep(:) = -coeff_rest_en*vel_norm_app(:)
      vel_tang_sep(:) =  coeff_rest_et*vel_tang_app(:)

      des_vel_new(ll,:) = vel_norm_sep(:) + vel_tang_sep(:)

      IF(dot_product(wall_norm, minvel) > 0.0d0)THEN

! Projection of rebound velocity onto mininum velocity.
         projgrav = dot_product(minvel, des_vel_new(ll,:))*oominvel_mag
         projgrav = minvel*projgrav

         IF(dot_product(projgrav, projgrav) < minvel_mag) THEN
            DO lc=1,3
               IF(minvel(lc) > small_number) THEN
                  des_vel_new(ll,lc)=max(des_vel_new(ll,lc),minvel(lc))
               ELSEIF(minvel(lc) < -small_number) THEN
                  des_vel_new(ll,lc)=min(des_vel_new(ll,lc),minvel(lc))
               ENDIF
            ENDDO

         ENDIF
      ENDIF

      RETURN
      END SUBROUTINE pic_reflect_part


!
      SUBROUTINE out_this_stl(LC)

      Use usr

! STL Vertices
      use stl, only: vertex
! STL Facet normals
      use stl, only: norm_face


      IMPLICIT NONE
!-----------------------------------------------
      integer, INTENT(IN) :: lc

      logical :: lExists
      character(len=8) :: IDX


      write(idx,"(I8.8)") lc
      inquire(file='geo_'//idx//'.stl', exist=lexists)

      IF(lexists) RETURN

      open(unit=555, file='geo_'//idx//'.stl', status='UNKNOWN')
      write(555,*) 'solid vcg'
      write(555,*) '   facet normal ', norm_face(:,lc)
      write(555,*) '      outer loop'
      write(555,*) '         vertex ', vertex(1,1:3,lc)
      write(555,*) '         vertex ', vertex(2,1:3,lc)
      write(555,*) '         vertex ', vertex(3,1:3,lc)
      write(555,*) '      endloop'
      write(555,*) '   endfacet'
      close(555)

      RETURN
      END SUBROUTINE out_this_stl




!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  SUBROUTINE: CHECK_IF_PARCEL_OVERLAPS_STL                            C
!  Authors: Rahul Garg                               Date: 21-Mar-2014 C
!                                                                      C
!  Purpose: This subroutine is special written to check if a particle  C
!          overlaps any of the STL faces. The routine exits on         C
!          detecting an overlap. It is called after initial            C
!          generation of lattice configuration to remove out of domain C
!          particles                                                   C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE check_if_parcel_overlaps_stl(POS, OVERLAP_EXISTS)

      USE discretelement, only: max_radius
      USE geometry, only: do_k

      USE stl_functions_des
      USE desgrid
      USE stl

      Implicit none

      DOUBLE PRECISION, INTENT(IN) :: POS(3)
      LOGICAL, INTENT(OUT) :: OVERLAP_EXISTS

      INTEGER I, J, K, IJK, NF, LC

! line is parameterized as p = p_ref + t * dir_line, t is line_param
      DOUBLE PRECISION :: LINE_T
      DOUBLE PRECISION :: DIST(3), RADSQ
      double precision :: vv(3)

      k = 1
      IF(do_k) k = min(dg_kend2,max(dg_kstart2,kofpos(pos(3))))
      j = min(dg_jend2,max(dg_jstart2,jofpos(pos(2))))
      i = min(dg_iend2,max(dg_istart2,iofpos(pos(1))))

      ijk = dg_funijk(i,j,k)
      IF(facets_at_dg(ijk)%COUNT < 1) RETURN

      overlap_exists = .true.

! Pad the max radius to keep parcels sufficiently away from walls.
      radsq = (1.1d0*max_radius)**2

! 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 side.
      DO lc = 1, facets_at_dg(ijk)%COUNT
         nf = facets_at_dg(ijk)%ID(lc)

         vv = vertex(1,:,nf)
         CALL intersectlnplane(pos, norm_face(:,nf), &
            vv, norm_face(:,nf), line_t)

! Orthogonal projection puts the point outside of the domain or less than
! one particle radius to the facet.
         dist = line_t*norm_face(:,nf)
         IF(line_t > zero .OR. dot_product(dist,dist) <= radsq) RETURN

      ENDDO

      overlap_exists = .false.

      RETURN
      END SUBROUTINE check_if_parcel_overlaps_stl


!----------------------------------------------------------------------!
!                                                                      !
!                                                                      !
!                                                                      !
!                                                                      !
!----------------------------------------------------------------------!
      SUBROUTINE write_this_facet_and_parcel(FID, position, velocity)
      USE run
      USE param1
      USE discretelement
      USE geometry
      USE compar
      USE constant
      USE cutcell
      USE funits
      USE indices
      USE physprop
      USE parallel
      USE stl
      USE stl_functions_des
      Implicit none
      !facet id and particle id
      double precision, intent(in), dimension(dimn) :: position, velocity
      Integer, intent(in) :: fid
      Integer :: stl_unit, vtp_unit , k
      CHARACTER(LEN=100) :: stl_fname, vtp_fname
      real :: temp_array(3)

      stl_unit = 1001
      vtp_unit = 1002

      WRITE(vtp_fname,'(A,"_OFFENDING_PARTICLE",".vtp")') trim(run_name)
      WRITE(stl_fname,'(A,"_STL_FACE",".stl")') trim(run_name)

      open(vtp_unit, file = vtp_fname, form='formatted',convert='big_endian')
      open(stl_unit, file = stl_fname, form='formatted',convert='big_endian')

      write(vtp_unit,"(a)") '<?xml version="1.0"?>'
      write(vtp_unit,"(a,es24.16,a)") '<!-- time =',s_time,'s -->'
      write(vtp_unit,"(a,a)") '<VTKFile type="PolyData"',&
           ' version="0.1" byte_order="LittleEndian">'
      write(vtp_unit,"(3x,a)") '<PolyData>'
      write(vtp_unit,"(6x,a,i10.10,a,a)")&
           '<Piece NumberOfPoints="',1,'" NumberOfVerts="0" ',&
           'NumberOfLines="0" NumberOfStrips="0" NumberOfPolys="0">'
      write(vtp_unit,"(9x,a)")&
           '<PointData Scalars="Diameter" Vectors="Velocity">'
      write(vtp_unit,"(12x,a)")&
           '<DataArray type="Float32" Name="Diameter" format="ascii">'
      write (vtp_unit,"(15x,es13.6)") (1.d0)
      write(vtp_unit,"(12x,a)") '</DataArray>'

      temp_array = zero
      temp_array(1:dimn) = velocity(1:dimn)
      write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
           'Name="Velocity" NumberOfComponents="3" format="ascii">'
      write (vtp_unit,"(15x,3(es13.6,3x))")&
           ((temp_array(k)),k=1,3)
      write(vtp_unit,"(12x,a,/9x,a)") '</DataArray>','</PointData>'
      ! skip cell data
      write(vtp_unit,"(9x,a)") '<CellData></CellData>'

      temp_array = zero
      temp_array(1:dimn) = position(1:dimn)
      write(vtp_unit,"(9x,a)") '<Points>'
      write(vtp_unit,"(12x,a,a)") '<DataArray type="Float32" ',&
           'Name="Position" NumberOfComponents="3" format="ascii">'
      write (vtp_unit,"(15x,3(es13.6,3x))")&
           ((temp_array(k)),k=1,3)
      write(vtp_unit,"(12x,a,/9x,a)")'</DataArray>','</Points>'
      ! Write tags for data not included (vtp format style)
      write(vtp_unit,"(9x,a,/9x,a,/9x,a,/9x,a)")'<Verts></Verts>',&
           '<Lines></Lines>','<Strips></Strips>','<Polys></Polys>'
      write(vtp_unit,"(6x,a,/3x,a,/a)")&
           '</Piece>','</PolyData>','</VTKFile>'

      !Now write the facet info

      write(stl_unit,*)'solid vcg'

      write(stl_unit,*) '   facet normal ', norm_face(1:3,fid)
      write(stl_unit,*) '      outer loop'
      write(stl_unit,*) '         vertex ', vertex(1,1:3,fid)
      write(stl_unit,*) '         vertex ', vertex(2,1:3,fid)
      write(stl_unit,*) '         vertex ', vertex(3,1:3,fid)
      write(stl_unit,*) '      endloop'
      write(stl_unit,*) '   endfacet'

      write(stl_unit,*)'endsolid vcg'

      close(vtp_unit, status = 'keep')
      close(stl_unit, status = 'keep')
      write(*,*) 'wrote a facet and a parcel. now waiting'
      read(*,*)
      end SUBROUTINE write_this_facet_and_parcel