interpolate_pic.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_PS_GRAD_PIC                                        !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose:                                                            !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE interpolate_pic

      use particle_filter, only: des_interp_scheme_enum
      use particle_filter, only: des_interp_garg

      IMPLICIT NONE

      SELECT CASE(des_interp_scheme_enum)
      CASE(des_interp_garg) ; CALL interpolate_pic0
      CASE DEFAULT; CALL interpolate_pic1
      END SELECT

      RETURN
      END SUBROUTINE interpolate_pic



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: INTERPOLATE_PIC_NONE                                    !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose:                                                            !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE interpolate_pic1

      use discretelement, only: max_pip, pijk
      USE geometry, only: do_k
      use fldvar, only: ep_g, u_s, v_s, w_s
      use functions, only: fluid_at
      use functions, only: is_nonexistent
      use mfix_pic, only: avgvel => avgsolvel_p
      use mfix_pic, only: ps_grad, ps_force_pic
      use mfix_pic, only: epg_p
! Flag to use interpolation
      use particle_filter, only: des_interp_on
! Interpolation cells and weights
      use particle_filter, only: filter_cell, filter_weight
      use physprop, only: mmax

      implicit none

! Local Variables
!---------------------------------------------------------------------//
! Interpolation weight
      DOUBLE PRECISION :: WEIGHT
! temporary variables used to calculate pressure at scalar cell edge
      INTEGER :: IJK, NP, LC
! Loop bound for filter
      INTEGER :: LP_BND
!......................................................................!

! Loop bounds for interpolation.
      lp_bnd = merge(27,9,do_k)


      DO np=1,max_pip

! Skip parcels that don't exist or are in a non-fluid cell.
         IF(is_nonexistent(np) .OR. .NOT.fluid_at(pijk(np,4))) THEN
            avgvel(:,np) = 0.0d0
            ps_grad(:,np) = 0.0d0
            epg_p(np) = 1.0d0

! Interpolated solids pressure and average solids velocity.
         ELSEIF(des_interp_on) THEN
            avgvel(:,np) = 0.0d0
            ps_grad(:,np) = 0.0d0
            epg_p(np) = 0.0d0
            DO lc=1,lp_bnd
               ijk = filter_cell(lc,np)
               weight = filter_weight(lc,np)
! Gas phase velocity.
               avgvel(1,np) = avgvel(1,np) + u_s(ijk,mmax+1)*weight
               avgvel(2,np) = avgvel(2,np) + v_s(ijk,mmax+1)*weight
               avgvel(3,np) = avgvel(3,np) + w_s(ijk,mmax+1)*weight
! Particle normal stress.
               ps_grad(:,np) = ps_grad(:,np)+ps_force_pic(:,ijk)*weight
! TO BE REMOVED
               epg_p(np) = epg_p(np) + ep_g(ijk)*weight
            ENDDO

! Centroid method. 
        ELSE
            ijk = pijk(np,4)
! Gas phase velocity.
            avgvel(1,np) = u_s(ijk,mmax+1)
            avgvel(2,np) = v_s(ijk,mmax+1)
            avgvel(3,np) = w_s(ijk,mmax+1)
! Particle normal stress.
            ps_grad(:,np) = ps_force_pic(:,ijk)
! TO BE REMOVED
            epg_p(np) = ep_g(ijk)
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE interpolate_pic1


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_PS_GRAD_PIC                                        !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose:                                                            !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE interpolate_pic0

      USE bc
      USE compar
      USE constant
      USE cutcell
      USE derived_types, only: pic
      USE discretelement
      USE fldvar, only: ep_g, u_g, v_g, w_g
      USE fun_avg
      USE functions
      USE geometry
      USE indices
      USE interpolation
      USE mfix_pic
      USE parallel
      USE param
      USE param1
      USE physprop
      USE run
      USE sendrecv

      implicit none

      ! general i, j, k indices
      INTEGER I, J, K, IJK, IPJK, IJPK, IJKP, IDIM, M

      ! temporary variables used to calculate pressure at scalar cell edge
      DOUBLE PRECISION avg_factor, VOL_TOT_VEC(3), VOL_TOT_SCAL

      integer :: korder, iw,ie,js,jn,kb,ktp, onew, pcell(3), cur_ijk, NP, nindx

      integer :: ii,jj,kk, ipjpk, ijpkp, ipjkp, ipjpkp

      double precision :: vol_ijk, vol_ipjk, vol_ijpk, vol_ipjpk
      double precision :: vol_ijkp, vol_ipjkp, vol_ijpkp, vol_ipjpkp


      CALL set_interpolation_scheme(2)

      korder = merge( 1, 2, no_k) !1+(DIMN-2)

! avg_factor=0.25 (in 3D) or =0.5 (in 2D)
      !avg_factor = 0.25d0*(dimn-2) + 0.5d0*(3-dimn)
      avg_factor = merge(0.5d0, 0.25d0, no_k)

      do ijk = ijkstart3,ijkend3

         if(.not.fluid_at(ijk) .or. pinc(ijk).eq.0) cycle
         i = i_of(ijk)
         j = j_of(ijk)
         k = k_of(ijk)

         pcell(1) = i-1
         pcell(2) = j-1
         pcell(3) = merge(1, k-1, no_k) ! =k-1 (in 3d) or =1 (in 2d)

         CALL set_interpolation_stencil(pcell,iw,ie,js,jn,kb,          &
            ktp,interp_scheme,merge(2,3,no_k),ordernew = onew)

! Compute velocity at grid nodes and set the geometric stencil
         DO k = 1, merge(1, onew, no_k)!  (3-DIMN)*1+(DIMN-2)*ONEW
         DO j = 1,onew
         DO i = 1,onew

            ii = iw + i-1
            jj = js + j-1
            kk = kb + k-1

            cur_ijk = funijk_map_c(ii,jj,kk)

            ipjk    = funijk_map_c(ii+1,jj,kk)
            ijpk    = funijk_map_c(ii,jj+1,kk)
            ipjpk   = funijk_map_c(ii+1,jj+1,kk)

            vol_ijk = zero
            vol_ipjk = zero
            vol_ijpk = zero
            vol_ipjpk = zero

            vol_ijkp = zero
            vol_ipjkp = zero
            vol_ijpkp = zero
            vol_ipjpkp = zero

            IF(fluid_at(cur_ijk)) vol_ijk   = vol(cur_ijk)
            IF(fluid_at(ipjk))    vol_ipjk  = vol(ipjk)
            IF(fluid_at(ijpk))    vol_ijpk  = vol(ijpk)
            IF(fluid_at(ipjpk))   vol_ipjpk = vol(ipjpk)

            IF(do_k) THEN
               ijkp   = funijk_map_c(ii,jj,kk+1)
               ijpkp  = funijk_map_c(ii,jj+1,kk+1)
               ipjkp  = funijk_map_c(ii+1,jj,kk+1)
               ipjpkp = funijk_map_c(ii+1,jj+1,kk+1)

               IF(fluid_at(ijkp))   vol_ijkp   = vol(ijkp)
               IF(fluid_at(ipjkp))  vol_ipjkp  = vol(ipjkp)
               IF(fluid_at(ijpkp))  vol_ijpkp  = vol(ijpkp)
               IF(fluid_at(ipjpkp)) vol_ipjpkp = vol(ipjpkp)
            ENDIF

            gstencil(i,j,k,1) = xe(ii)
            gstencil(i,j,k,2) = yn(jj)
            gstencil(i,j,k,3) = merge(dz(1), zt(kk), no_k)

            vol_tot_scal = zero

            vol_tot_scal = vol_ijk + vol_ipjk + vol_ijpk + vol_ipjpk + &
               vol_ijkp + vol_ipjkp + vol_ijpkp + vol_ipjpkp

            vol_tot_vec = zero

            vol_tot_vec(1) = vol(cur_ijk) + vol(ijpk)
            vol_tot_vec(2) = vol(cur_ijk) + vol(ipjk)

            DO m = mmax+1,des_mmax+mmax
               vel_sol_stencil(i,j,k,1,m) = &
                  pic_u_s(cur_ijk,m)*vol(cur_ijk) + &
                  pic_u_s(ijpk,m)*vol(ijpk)

               vel_sol_stencil(i,j,k,2,m) = &
                  pic_v_s(cur_ijk,m)*vol(cur_ijk) + &
                  pic_v_s(ipjk,m)*vol(ipjk)
            ENDDO

            sstencil(i,j,k) = ep_g(cur_ijk)*vol_ijk + &
               ep_g(ipjk)*vol_ipjk + ep_g(ijpk)*vol_ijpk + &
               ep_g(ipjpk)*vol_ipjpk

            psgradstencil(i,j,k,1) = &
               ps_force_pic(1,cur_ijk)*vol(cur_ijk) + &
               ps_force_pic(1,ijpk)*vol(ijpk)

            psgradstencil(i,j,k,2) = &
               ps_force_pic(2,cur_ijk)*vol(cur_ijk) + &
               ps_force_pic(2,ipjk)*vol(ipjk)

            vstencil(i,j,k,1) = &
               u_g(cur_ijk)*vol(cur_ijk) + u_g(ijpk)*vol(ijpk)

            vstencil(i,j,k,2) = &
               v_g(cur_ijk)*vol(cur_ijk) + v_g(ipjk)*vol(ipjk)

            IF(do_k) THEN
               vol_tot_vec(1) = vol_tot_vec(1) + vol(ijkp) + vol(ijpkp)
               vol_tot_vec(2) = vol_tot_vec(2) + vol(ijkp) + vol(ipjkp)
               vol_tot_vec(3) = &
                  vol(cur_ijk) + vol(ipjk) + vol(ijpk) + vol(ipjpk)

               sstencil(i,j,k) = sstencil(i,j,k) + &
                  ep_g(ijkp)*vol_ijkp + ep_g(ipjkp)*vol_ipjkp + &
                  ep_g(ijpkp)*vol_ijpkp + ep_g(ipjpkp)*vol_ipjpkp

               psgradstencil(i,j,k,1) = psgradstencil(i,j,k,1) + &
                  ps_force_pic(1,ijkp)*vol(ijkp) + &
                  ps_force_pic(1,ijpkp)*vol(ijpkp)

               psgradstencil(i,j,k,2) = psgradstencil(i,j,k,2) + &
                  ps_force_pic(2,ijkp)*vol(ijkp) + &
                  ps_force_pic(2,ipjkp)*vol(ipjkp)

               psgradstencil(i,j,k,3) = &
                  ps_force_pic(3,cur_ijk)*vol(cur_ijk)+ &
                  ps_force_pic(3,ijpk)*vol(ijpk)+       &
                  ps_force_pic(3,ipjk)*vol(ipjk)+       &
                  ps_force_pic(3,ipjpk)*vol(ipjpk)

               vstencil(i,j,k,1) = vstencil(i,j,k,1) + &
                  u_g(ijkp)*vol(ijkp) + u_g(ijpkp)*vol(ijpkp)

               vstencil(i,j,k,2) = vstencil(i,j,k,2) + &
                  v_g(ijkp)*vol(ijkp) + v_g(ipjkp)*vol(ipjkp)

               vstencil(i,j,k,3) = w_g(cur_ijk)*vol(cur_ijk)+&
                  w_g(ijpk)*vol(ijpk) + w_g(ipjk)*vol(ipjk) + &
                  w_g(ipjpk)*vol(ipjpk)

               DO m = mmax+1, des_mmax+mmax
                  vel_sol_stencil(i,j,k,1,m) = &
                     vel_sol_stencil(i,j,k,1,m) + &
                     pic_u_s(ijkp,m)*vol(ijkp) + & 
                     pic_u_s(ijpkp,m)*vol(ijpkp)

                  vel_sol_stencil(i,j,k,2, m) = &
                     vel_sol_stencil(i,j,k,2,m) + &
                     pic_v_s(ijkp,m)*vol(ijkp) +&
                     pic_v_s(ipjkp,m)*vol(ipjkp)

                  vel_sol_stencil(i,j,k,3, m) = &
                     pic_w_s(cur_ijk,m)*vol(cur_ijk) + &
                     pic_w_s(ijpk,m)*vol(ijpk) + &
                     pic_w_s(ipjk,m)*vol(ipjk) + &
                     pic_w_s(ipjpk,m)*vol(ipjpk)
               ENDDO
            ELSE
               psgradstencil(i,j,k,3) = 0.d0
               vel_sol_stencil(i,j,k,3, mmax+1:des_mmax+mmax) = 0.d0
               vstencil(i,j,k,3) = 0.d0
            ENDIF




            DO idim = 1, merge(2,3,no_k)
               IF(vol_tot_vec(idim).GT.zero)  THEN
                  psgradstencil(i,j,k,idim) = &
                     psgradstencil(i,j,k,idim)/vol_tot_vec(idim)

                  vel_sol_stencil(i,j,k,idim, mmax+1:des_mmax+mmax) = &
                     vel_sol_stencil(i,j,k,idim, mmax+1:des_mmax+mmax)/&
                     vol_tot_vec(idim)

                  vstencil(i,j,k,idim) = &
                     vstencil(i,j,k,idim)/vol_tot_vec(idim)
               ENDIF
            ENDDO


            IF(vol_tot_scal.GT.zero) &
               sstencil(i,j,k) = sstencil(i,j,k)/vol_tot_scal

         enddo
         enddo
         enddo



! Loop through particles in the cell
         do nindx = 1,pinc(ijk)
            np = pic(ijk)%p(nindx)
            m = pijk(np,5)

! Interpolate particle stress to particle position
            if(no_k) then !2-D
               call interpolator(gstencil(1:onew,1:onew,1,1:dimn), &
               psgradstencil(1:onew,1:onew,1,1:dimn), &
               des_pos_new(np,1:dimn),ps_grad(1:dimn,np),  &
               onew,interp_scheme,weightp)
            else !3-D, diff in psgradstencil size
               call interpolator(gstencil(1:onew,1:onew,1:onew,1:dimn), &
               psgradstencil(1:onew,1:onew,1:onew,1:dimn), &
               des_pos_new(np,1:dimn),ps_grad(1:dimn,np),  &
               onew,interp_scheme,weightp)
            endif

            do idim = 1,  merge(2,3,no_k)
               avgsolvel_p(idim,np) = array_dot_product(               &
                  vel_sol_stencil(:,:,:,idim,m),weightp(:,:,:))
            ENDDO

            epg_p(np) = array_dot_product(sstencil(:,:,:),weightp(:,:,:))

         ENDDO
      ENDDO


      END SUBROUTINE interpolate_pic0