source_pp_g.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOURCE_Pp_g                                             C
!  Purpose: Determine source terms for pressure correction equation.   C
!                                                                      C
!  Notes: The off-diagonal coefficients are positive. The center       C
!         coefficient and the source vector are negative. See          C
!         conv_Pp_g                                                    C
!                                                                      C
!  Author: M. Syamlal                                 Date: 21-JUN-96  C
!                                                                      C
!  Revision Number: 1                                                  C
!  Purpose: To incorporate Cartesian grid modifications                C
!  Author: Jeff Dietiker                              Date: 01-Jul-09  C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE source_pp_g(A_M, B_M, B_MMAX)

! Modules
!---------------------------------------------------------------------//
      USE bc, ONLY: ijk_p_g
      USE compar, ONLY: ijkstart3, ijkend3
      USE cutcell, ONLY: cartesian_grid, a_upg_e, a_vpg_n, a_wpg_t
      USE fldvar, ONLY: u_g, v_g, w_g, u_s, v_s, w_s, rop_g, rop_s
      USE fldvar, only: rop_go, rop_so, ro_g
      USE geometry, ONLY: vol
      USE param, only: dimension_3, dimension_m
      USE param, only: east, west, south, north, top, bottom
      USE param1, ONLY: small_number, one, zero, undefined
      USE pgcor, ONLY: d_e, d_n, d_t
      USE physprop, ONLY: close_packed, mmax, ro_g0
      USE run, ONLY: shear, odt, undefined_i
      USE rxns, ONLY: sum_r_g, sum_r_s
      USE vshear, ONLY: vsh
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
! maximum term in b_m expression
      DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(dimension_3, 0:dimension_m)

! Local Variables
!---------------------------------------------------------------------//
! solids phase index
      INTEGER :: M
! Indices
      INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
! terms of bm expression
      DOUBLE PRECISION bma, bme, bmw, bmn, bms, bmt, bmb, bmr
! error message
      CHARACTER(LEN=80) :: LINE(1)
!---------------------------------------------------------------------//
! loezos
! update to true velocity
      IF (shear) THEN
!!$omp parallel do private(IJK)
         DO ijk = ijkstart3, ijkend3
            IF (fluid_at(ijk)) THEN
               v_g(ijk)=v_g(ijk)+vsh(ijk)
            ENDIF
         ENDDO
      ENDIF

! Calculate convection-diffusion fluxes through each of the faces

!$omp parallel default(none) &
!$omp          private(IJK, IMJK, IJMK, IJKM, M, bma,bme,bmw,bmn,&
!$omp                  bms,bmt,bmb,bmr,line)  &
!$omp          shared(ijkstart3,ijkend3,cartesian_grid,rop_g,rop_go,&
!$omp                 rop_s,rop_so,vol,odt,u_g,v_g,w_g,u_s,v_s,w_s,b_m, &
!$omp                 b_mmax,d_e,d_n,d_t,a_m,a_upg_e,a_vpg_n,a_wpg_t,&
!$omp                 mmax,close_packed,sum_r_s,sum_r_g,ro_g0)
!$omp do
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ijkm = km_of(ijk)

            bma = (rop_g(ijk)-rop_go(ijk))*vol(ijk)*odt
            bme = a_m(ijk,east,0)*u_g(ijk)
            bmw = a_m(ijk,west,0)*u_g(imjk)
            bmn = a_m(ijk,north,0)*v_g(ijk)
            bms = a_m(ijk,south,0)*v_g(ijmk)
            bmt = a_m(ijk,top,0)*w_g(ijk)
            bmb = a_m(ijk,bottom,0)*w_g(ijkm)
            bmr = sum_r_g(ijk)*vol(ijk)
            b_m(ijk,0) = -((-(bma + bme - bmw + &
                              bmn - bms + bmt - bmb ))+ bmr )
            b_mmax(ijk,0) = max(abs(bma),abs(bme),abs(bmw),&
               abs(bmn),abs(bms),abs(bmt),abs(bmb),abs(bmr) )

            a_m(ijk,east,0) = a_m(ijk,east,0)*d_e(ijk,0)
            a_m(ijk,west,0) = a_m(ijk,west,0)*d_e(imjk,0)
            a_m(ijk,north,0) = a_m(ijk,north,0)*d_n(ijk,0)
            a_m(ijk,south,0) = a_m(ijk,south,0)*d_n(ijmk,0)
            a_m(ijk,top,0) = a_m(ijk,top,0)*d_t(ijk,0)
            a_m(ijk,bottom,0) = a_m(ijk,bottom,0)*d_t(ijkm,0)

            IF(cartesian_grid) THEN
               a_m(ijk,east,0) = a_m(ijk,east,0) * a_upg_e(ijk)
               a_m(ijk,west,0) = a_m(ijk,west,0) * a_upg_e(imjk)
               a_m(ijk,north,0) = a_m(ijk,north,0) * a_vpg_n(ijk)
               a_m(ijk,south,0) = a_m(ijk,south,0) * a_vpg_n(ijmk)
               a_m(ijk,top,0) = a_m(ijk,top,0) * a_wpg_t(ijk)
               a_m(ijk,bottom,0) = a_m(ijk,bottom,0) * a_wpg_t(ijkm)
            ENDIF

! solids phase terms are needed for those solids phases that do not
! become close packed. these terms are used to essentially make the
! matrix equation for gas pressure correction a mixture pressure
! correction equation by adding solids phases that do not become
! close packed
            DO m = 1, mmax
               IF (.NOT.close_packed(m)) THEN
                  b_m(ijk,0) = b_m(ijk,0) - &
                     ((-((rop_s(ijk,m)-rop_so(ijk,m))*vol(ijk)*odt+&
                     a_m(ijk,east,m)*u_s(ijk,m)-a_m(ijk,west,m)*u_s(imjk,m)+&
                     a_m(ijk,north,m)*v_s(ijk,m)-a_m(ijk,south,m)*v_s(ijmk,m)+&
                     a_m(ijk,top,m)*w_s(ijk,m)-a_m(ijk,bottom,m)*w_s(ijkm,m)))+&
                     sum_r_s(ijk,m)*vol(ijk))

                  IF(.NOT.cartesian_grid) THEN
                     a_m(ijk,east,0) = a_m(ijk,east,0) + a_m(ijk,east,m)*d_e(ijk,m)
                     a_m(ijk,west,0) = a_m(ijk,west,0) + a_m(ijk,west,m)*d_e(imjk,m)
                     a_m(ijk,north,0) = a_m(ijk,north,0) + a_m(ijk,north,m)*d_n(ijk,m)
                     a_m(ijk,south,0) = a_m(ijk,south,0) + a_m(ijk,south,m)*d_n(ijmk,m)
                     a_m(ijk,top,0) = a_m(ijk,top,0) + a_m(ijk,top,m)*d_t(ijk,m)
                     a_m(ijk,bottom,0) = a_m(ijk,bottom,0) + a_m(ijk,bottom,m)*d_t(ijkm,m)
                  ELSE
                     a_m(ijk,east,0) = a_m(ijk,east,0) + a_m(ijk,east,m)*d_e(ijk,m)  * a_upg_e(ijk)
                     a_m(ijk,west,0) = a_m(ijk,west,0) + a_m(ijk,west,m)*d_e(imjk,m) * a_upg_e(imjk)
                     a_m(ijk,north,0) = a_m(ijk,north,0) + a_m(ijk,north,m)*d_n(ijk,m)  * a_vpg_n(ijk)
                     a_m(ijk,south,0) = a_m(ijk,south,0) + a_m(ijk,south,m)*d_n(ijmk,m) * a_vpg_n(ijmk)
                     a_m(ijk,top,0) = a_m(ijk,top,0) + a_m(ijk,top,m)*d_t(ijk,m)  * a_wpg_t(ijk)
                     a_m(ijk,bottom,0) = a_m(ijk,bottom,0) + a_m(ijk,bottom,m)*d_t(ijkm,m) * a_wpg_t(ijkm)
                  ENDIF
               ENDIF
            ENDDO

            a_m(ijk,0,0) = -(a_m(ijk,east,0)+a_m(ijk,west,0)+&
                             a_m(ijk,north,0)+a_m(ijk,south,0)+&
                             a_m(ijk,top,0)+a_m(ijk,bottom,0))

            IF (abs(a_m(ijk,0,0)) < small_number) THEN
               IF (abs(b_m(ijk,0)) < small_number) THEN
                  a_m(ijk,0,0) = -one
                  b_m(ijk,0) = zero
               ELSEIF (ro_g0 .NE. undefined) THEN !This is an error only in incompressible flow
!!$omp             critical
                  WRITE (line, '(A,I6,A,I1,A,G12.5)') 'Error: At IJK = ', ijk, &
                     ' M = ', 0, ' A = 0 and b = ', b_m(ijk,0)
                  CALL write_error ('SOURCE_Pp_g', line, 1)
!!$omp             end critical
               ENDIF
            ENDIF


         ELSE   ! if/else branch .not.fluid_at(ijk)
! set the value (correction) in all wall and flow boundary cells to zero
! note the matrix coefficients and source vector should already be zero
! from the initialization of A and B but the following ensures the zero
! value.
            a_m(ijk,east,0) = zero
            a_m(ijk,west,0) = zero
            a_m(ijk,north,0) = zero
            a_m(ijk,south,0) = zero
            a_m(ijk,top,0) = zero
            a_m(ijk,bottom,0) = zero
            a_m(ijk,0,0) = -one
            b_m(ijk,0) = zero
         ENDIF   ! end if/else branch fluid_at(ijk)
      ENDDO    ! end do loop (ijk=ijkstart3,ijkend3)
!$omp end parallel


! loezos
      IF (shear) THEN
!!$omp parallel do private(IJK)
         DO ijk = ijkstart3, ijkend3
            IF (fluid_at(ijk)) THEN
               v_g(ijk)=v_g(ijk)-vsh(ijk)
            ENDIF
         ENDDO
      ENDIF

! Modification for compressibility 
      IF (ro_g0 == undefined) CALL compressible_pp_g(a_m)

! Remove the asymmetry in matrix caused by the pressure outlet or inlet
! boundaries.  Because the P' at such boundaries is zero we may set the
! coefficient in the neighboring fluid cell to zero without affecting
! the linear equation set.
!!$omp    parallel do                                                     &
!!$omp&   private(IJK,IMJK, IPJK, IJMK, IJPK, IJKM, IJKP)
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN
            imjk = im_of(ijk)
            ipjk = ip_of(ijk)
            ijmk = jm_of(ijk)
            ijpk = jp_of(ijk)
            ijkm = km_of(ijk)
            ijkp = kp_of(ijk)
! Cutting the neighbor link between fluid cell and adjacent p_flow_at cell
            if(p_flow_at(imjk)) a_m(ijk, west, 0) = zero
            if(p_flow_at(ipjk)) a_m(ijk, east, 0) = zero
            if(p_flow_at(ijmk)) a_m(ijk, south, 0) = zero
            if(p_flow_at(ijpk)) a_m(ijk, north, 0) = zero
            if(p_flow_at(ijkm)) a_m(ijk, bottom, 0) = zero
            if(p_flow_at(ijkp)) a_m(ijk, top, 0) = zero
         ENDIF
      ENDDO

! Specify P' to zero for incompressible flows. Check set_bc0
! for details on selection of IJK_P_g.
      IF (ijk_p_g /= undefined_i) THEN
         b_m(ijk_p_g,0) = zero
         a_m(ijk_p_g,:,0) = zero
         a_m(ijk_p_g,0,0) = -one
      ENDIF

      RETURN

CONTAINS

      include 'functions.inc'

END SUBROUTINE source_pp_g


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: COMPRESSIBLE_PP_G                                       C
!  Purpose: Correction for incompressible flows.                       C
!                                                                      C
!  Notes: A HS_CORRECT option is available as a better approximation   C
!  for high speed flows because it considers density changes in the    C
!  neighboring C cells. However, the code runs faster without it for   C
!  low speed flows. The gas phase mass balance cannot be maintained    C
!  to machine precision with this approximation.                       C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE compressible_pp_g(A_M)

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      use eos, only: droodp_g
      use fldvar, only: ep_g, u_g, v_g, w_g, rop_g, ro_g, p_g
      use functions, only: fluid_at, im_of, jm_of, km_of
      use functions, only: east_of, west_of, north_of, south_of
      use functions, only: top_of, bottom_of
      use geometry, only: vol, ayz, axz, axy, do_k
      use param, only: dimension_3, dimension_m
      use param, only: east, west, south, north, top, bottom
      use param1, only: one
      use run, only: discretize, odt
      use ur_facs, only: ur_fac
      use xsi
      IMPLICIT NONE

! Parameters
!---------------------------------------------------------------------//
      LOGICAL, PARAMETER :: HS_CORRECT = .false.

! Dummy arguments
!---------------------------------------------------------------------//
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)

! Local variables
!---------------------------------------------------------------------//
! loezos: used for including shearing
      integer :: incr
! temporary use of global arrays:
! xsi_array: convection weighting factors
      DOUBLE PRECISION :: XSI_e(dimension_3)
      DOUBLE PRECISION :: XSI_n(dimension_3)
      DOUBLE PRECISION :: XSI_t(dimension_3)
! under relaxation factor for pressure
      double precision :: fac
! Indices
      INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
      INTEGER :: IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
!---------------------------------------------------------------------//

! since p_g = p_g* + ur_fac * pp_g
      fac = ur_fac(1)  

     IF (.NOT.hs_correct) THEN
!!$omp    parallel do  &
!!$omp&   private(IJK)
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN
! account for change in density w.r.t. pressure
            a_m(ijk,0,0) = a_m(ijk,0,0) - ur_fac(1)*&
               droodp_g(ro_g(ijk),p_g(ijk))*ep_g(ijk)*vol(ijk)*odt
         ENDIF
      ENDDO
      ENDIF

      IF (hs_correct) THEN
! since p_g = p_g* + ur_fac * pp_g
! loezos
      incr=0
      CALL calc_xsi(discretize(1),rop_g,u_g,v_g,w_g,xsi_e,xsi_n,xsi_t,incr)

!!$omp    parallel do  &
!!$omp&   private(IJK,I,J,K, &
!!$omp&            IMJK,IJMK,IJKM,IJKE,IJKW,IJKN,IJKS,IJKT,IJKB)
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk)) THEN
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ijkm = km_of(ijk)
            ijke = east_of(ijk)
            ijkw = west_of(ijk)
            ijkn = north_of(ijk)
            ijks = south_of(ijk)
            ijkt = top_of(ijk)
            ijkb = bottom_of(ijk)
            a_m(ijk,0,0) = a_m(ijk,0,0) - fac*droodp_g(ro_g(ijk),p_g(ijk))*&
               ep_g(ijk)*((one - xsi_e(ijk))*u_g(ijk)*ayz(ijk)-&
                                 xsi_e(imjk)*u_g(imjk)*ayz(imjk)+&
                             (one-xsi_n(ijk))*v_g(ijk)*axz(ijk)-&
                                  xsi_n(ijmk)*v_g(ijmk)*axz(ijmk))

            a_m(ijk,east,0) = a_m(ijk,east,0) - ep_g(ijke)*fac*&
               droodp_g(ro_g(ijke),p_g(ijke))*xsi_e(ijk)*u_g(ijk)*ayz(ijk)
            a_m(ijk,west,0) = a_m(ijk,west,0) + ep_g(ijkw)*fac*&
               droodp_g(ro_g(ijkw),p_g(ijkw))*(one - xsi_e(imjk))*u_g(imjk)*ayz(imjk)
            a_m(ijk,north,0) = a_m(ijk,north,0) - ep_g(ijkn)*fac*&
               droodp_g(ro_g(ijkn),p_g(ijkn))*xsi_n(ijk)*v_g(ijk)*axz(ijk)
            a_m(ijk,south,0) = a_m(ijk,south,0) + ep_g(ijks)*fac*&
               droodp_g(ro_g(ijks),p_g(ijks))*(one - xsi_n(ijmk))*v_g(ijmk)*axz(ijmk)
            IF (do_k) THEN
               a_m(ijk,0,0) = a_m(ijk,0,0) - fac*droodp_g(ro_g(ijk),p_g(ijk))*&
                  ep_g(ijk)*((one - xsi_t(ijk))*w_g(ijk)*axy(ijk)-&
                                    xsi_t(ijkm)*w_g(ijkm)*axy(ijkm))
               a_m(ijk,top,0) = a_m(ijk,top,0) - ep_g(ijkt)*fac*&
                  droodp_g(ro_g(ijkt),p_g(ijkt))*xsi_t(ijk)*w_g(ijk)*axy(ijk)
               a_m(ijk,bottom,0) = a_m(ijk,bottom,0) + ep_g(ijkb)*fac*&
                  droodp_g(ro_g(ijkb),p_g(ijkb))*(one - xsi_t(ijkm))*w_g(ijkm)*axy(ijkm)
            ENDIF

         ENDIF   !end if (fluid_at(ijk))
      ENDDO    ! end do (ijk=ijkstart3,ijkend3)
      ENDIF   ! end if (hs_correct)

      RETURN
      END SUBROUTINE compressible_pp_g