calc_ps_grad_pic.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_PS_GRAD_PIC                                        !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose: Calculate the solids pressure graident.                    !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_ps_grad_pic

! Modules
!---------------------------------------------------------------------//
      use particle_filter, only: des_interp_scheme_enum
      use particle_filter, only: des_interp_garg
      use mfix_pic, only: pic_p_s
      use mfix_pic, only: ps_force_pic
      IMPLICIT NONE

!......................................................................!

      SELECT CASE(des_interp_scheme_enum)
      CASE(des_interp_garg); CALL calc_ps_grad_pic0
      CASE DEFAULT; CALL calc_grad_des(pic_p_s(:,1), ps_force_pic)
      END SELECT

      RETURN
      END SUBROUTINE calc_ps_grad_pic



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!  Subroutine: CALC_PS_GRAD_PIC0                                       !
!  Author: R. Garg                                                     !
!                                                                      !
!  Purpose: Calculate the solids pressure graident. This routine       !
!  stores the solid pressure graident at cell faces.                   !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_ps_grad_pic0

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      use compar, only: istart2, jstart2, kstart2
      use compar, only: istart3, jstart3, kstart3
      use compar, only: iend3, jend3, kend3
      use functions, only: fluid_at
      use indices, only: i_of, j_of, k_of
      use functions, only: ip_of, jp_of, kp_of
      use functions, only: funijk
      use functions, only: is_on_mype_owns
      use geometry, only: do_k, no_k
      use geometry, only: dx, dy, dz
      use geometry, only: imin2, jmin2, kmin2
      USE mfix_pic, only: pic_p_s, ps_force_pic
      use param1, only: zero
      use sendrecv, only: send_recv
      implicit none

! Local variables
!---------------------------------------------------------------------//
! general i, j, k indices
      INTEGER I, J, K, IJK, IPJK, IJPK, IJKP, IDIM
      integer :: I1, J1, K1
!......................................................................!

! Since EP_G is already shared across the processors, the pressure
! gradient calculation can be made a function call so that the extra
! communication of P_S can be avoided.

      DO ijk = ijkstart3, ijkend3
         ps_force_pic(:,ijk) = zero
         IF(.NOT.fluid_at(ijk)) cycle

         i = i_of(ijk)
         j = j_of(ijk)
         k = k_of(ijk)

         ipjk = ip_of(ijk)
         ijpk = jp_of(ijk)
         ijkp = kp_of(ijk)

! Calculate the solids pressure gradient at east face.
         IF(fluid_at(ipjk)) THEN
            ps_force_pic(1,ijk) = 2.0d0 *                              &
               (pic_p_s(ipjk,1) - pic_p_s(ijk,1)) /                    &
               (dx(i) + dx(i_of(ipjk)))
         ELSE
            IF(pic_p_s(ijk,1) > zero) THEN
               ps_force_pic(1,ijk) = 2.0d0*                            &
                  (pic_p_s(ipjk,1) - pic_p_s(ijk,1)) /                 &
                  (dx(i) + dx(i_of(ipjk)))
            ELSE
               ps_force_pic(1,ijk) = zero
            ENDIF
         ENDIF

! Calculate the solids pressure graident at the north face.
         IF(fluid_at(ijpk)) THEN
            ps_force_pic(2,ijk) = 2.0d0*                               &
               (pic_p_s(ijpk,1) - pic_p_s(ijk,1)) /                    &
               (dy(j)+dy(j_of(ijpk)))
         ELSE
            IF(pic_p_s(ijk,1) > zero) THEN
               ps_force_pic(2,ijk) = 2.0d0*                            &
                  (pic_p_s(ijpk,1) - pic_p_s(ijk,1))/                  &
                  (dy(j)+dy(j_of(ijpk)))
            ELSE
               ps_force_pic(2,ijk) = zero
            ENDIF
         ENDIF

! Calculate the solids pressure graident at the top face.
         IF(do_k) THEN
            IF(fluid_at(ijkp)) THEN
               ps_force_pic(3,ijk) = 2.0d0*                            &
                  (pic_p_s(ijkp,1) - pic_p_s(ijk,1))/                  &
                  (dz(k)+dz(k_of(ijkp)))
            ELSE
               IF(pic_p_s(ijk,1).GT.zero) then
                  ps_force_pic(3,ijk) = 2.0d0*&
                     (pic_p_s(ijkp,1) - pic_p_s(ijk,1))/               &
                     (dz(k)+dz(k_of(ijkp)))
               ELSE
                  ps_force_pic(3,ijk) = zero
               ENDIF
            ENDIF
         ENDIF
      ENDDO

! Compute the pressure gradients along the west domain bondary which
! is previously skipped.
      i1 = imin2
      IF(i1 == istart2) THEN
         DO k1 = kstart3, kend3
         DO j1 = jstart3, jend3
            IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
            ijk = funijk(i1,j1,k1)
            ipjk = ip_of(ijk)
            IF(pic_p_s(ipjk,1) > zero) THEN
               ps_force_pic(1,ijk) = 2.0d0*                            &
                  (pic_p_s(ipjk,1) - pic_p_s(ijk,1)) /                 &
                  (dx(i1) + dx(i_of(ipjk)))
            ELSE
               ps_force_pic(1,ijk) = zero
            ENDIF
         ENDDO
         ENDDO
      ENDIF

! Compute the pressure gradients along the south domain bondary which
! is previously skipped.
      j1 = jmin2
      IF(j1 == jstart2) THEN
         DO k1 = kstart3, kend3
         DO i1 = istart3, iend3
            IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
            ijk = funijk(i1,j1,k1)
            ijpk = jp_of(ijk)
            IF(pic_p_s(ijpk,1).GT.zero) then
               ps_force_pic(2,ijk) = 2.0d0*                            &
                  (pic_p_s(ijpk,1) - pic_p_s(ijk,1))/                  &
                  (dy(j) + dy(j_of(ijpk)))
            ELSE
               ps_force_pic(2,ijk) = zero
            ENDIF
         ENDDO
         ENDDO
      ENDIF

! Compute the pressure gradients along the south domain bondary which
! is previously skipped.
      IF(do_k) then
         k1 = kmin2
         IF(k1 == kstart2) THEN
            DO j1 = jstart3, jend3
            DO i1 = istart3, iend3
               IF (.NOT.is_on_mype_owns(i1,j1,k1)) cycle
               ijk = funijk(i1,j1,k1)
               ijkp = kp_of(ijk)
               IF(pic_p_s(ijkp,1).GT.zero) THEN
                  ps_force_pic(3,ijk) = 2.0d0*                         &
                     (pic_p_s(ijkp,1) - pic_p_s(ijk,1))/               &
                     (dz(k)+dz(k_of(ijkp)))
               ELSE
                  ps_force_pic(3,ijk) = zero
               ENDIF
            ENDDO
            ENDDO
         ENDIF
      ENDIF

      DO idim = 1, merge(2,3,no_k)
         CALL send_recv(ps_force_pic(idim,:),1)
      ENDDO

      END SUBROUTINE calc_ps_grad_pic0