source_phi.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOURCE_phi                                              C
!  Purpose: Determine source terms for phi eq. The terms               C
!     appear in the center coefficient and RHS vector. The center      C
!     coefficient and source vector are negative. The off-diagonal     C
!     coefficients are positive. S_p must be positive.                 C
!                                                                      C
!  Note: This routine is now restricted to Non-Negative scalers when   C
!     using deferred correction.                                       C
!                                                                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 30-APR-97  C
!  Reviewer:                                          Date:            C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE source_phi(S_P, S_C, EP, PHI, M, A_M, B_M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE scales
      USE physprop
      USE fldvar
      USE visc_s
      USE rxns
      USE run
      USE toleranc
      USE geometry
      USE indices
      USE is
      USE tau_s
      USE compar
      USE fun_avg
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Dummy Arguments
!-----------------------------------------------
! Source term on LHS.  Must be positive.
      DOUBLE PRECISION, INTENT(IN) :: S_p(dimension_3)
! Source term on RHS
      DOUBLE PRECISION, INTENT(IN) :: S_C(dimension_3)
! Phase volume fraction
      DOUBLE PRECISION, INTENT(IN) :: EP(dimension_3)
! Phi
      DOUBLE PRECISION, INTENT(IN) :: Phi(dimension_3)
! phase index
      INTEGER, INTENT(IN) :: M
! 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)
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: IJK, IMJK, IJMK, IJKM
!-----------------------------------------------

      DO ijk = ijkstart3, ijkend3

         IF (fluid_at(ijk)) THEN

! dilute flow
            IF (ep(ijk) <= dil_ep_s) THEN
               a_m(ijk,east,m) = zero
               a_m(ijk,west,m) = zero
               a_m(ijk,north,m) = zero
               a_m(ijk,south,m) = zero
               a_m(ijk,top,m) = zero
               a_m(ijk,bottom,m) = zero
               a_m(ijk,0,m) = -one
               b_m(ijk,m) = zero

! use the average boundary cell values to compute phi (sof, Aug 23 2005)
               imjk = im_of(ijk)
               ijmk = jm_of(ijk)
               ijkm = km_of(ijk)
               IF (ep(west_of(ijk)) > dil_ep_s .AND. &
                   .NOT.is_at_e(imjk)) a_m(ijk,west,m) = one
               IF (ep(east_of(ijk)) > dil_ep_s .AND. &
                   .NOT.is_at_e(ijk)) a_m(ijk,east,m) = one
               IF (ep(south_of(ijk)) > dil_ep_s .AND. &
                   .NOT.is_at_n(ijmk)) a_m(ijk,south,m) = one
               IF (ep(north_of(ijk)) > dil_ep_s .AND. &
                   .NOT.is_at_n(ijk)) a_m(ijk,north,m) = one
               IF(.NOT. no_k) THEN
                  IF (ep(bottom_of(ijk)) > dil_ep_s .AND. &
                     .NOT.is_at_t(ijkm)) a_m(ijk,bottom,m) = one
                  IF (ep(top_of(ijk)) > dil_ep_s .AND. &
                     .NOT.is_at_t(ijk)) a_m(ijk,top,m) = one
               ENDIF

               IF((a_m(ijk,west,m)+a_m(ijk,east,m)+&
                  a_m(ijk,south,m)+a_m(ijk,north,m)+ &
                  a_m(ijk,bottom,m)+a_m(ijk,top,m)) == zero) THEN
                  b_m(ijk,m) = -phi(ijk)
               ELSE
                  a_m(ijk,0,m) = -(a_m(ijk,east,m) + a_m(ijk,west,m) +&
                     a_m(ijk,north,m) + a_m(ijk,south,m) + &
                     a_m(ijk,top,m)+a_m(ijk,bottom,m))
               ENDIF

! Normal case
            ELSE

! Collect the terms
               a_m(ijk,0,m) = -(a_m(ijk,east,m)+a_m(ijk,west,m)+&
                                a_m(ijk,north,m)+a_m(ijk,south,m)+&
                                a_m(ijk,top,m)+a_m(ijk,bottom,m)+s_p(ijk))

! B_m and A_m are corrected in case deferred corrections computes B_m > S_c
! see CONV_DIF_PHI_DC.
               IF(b_m(ijk,m) < s_c(ijk) .OR. phi(ijk) == zero) THEN
                  b_m(ijk,m) = -s_c(ijk)+b_m(ijk,m)
               ELSE ! disable ELSE statememt if PHI can be negative
                  a_m(ijk,0,m) = a_m(ijk,0,m) - b_m(ijk,m)/phi(ijk)
                  b_m(ijk,m) = -s_c(ijk)
               ENDIF

            ENDIF   ! end if/else (ep_g(ijk)<=dil_ep_s)
         ENDIF    ! end if (fluid_at(ijk))
      ENDDO   ! end do (ijk=ijkstart3,ijkend3)


      RETURN
      END SUBROUTINE source_phi



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: POINT_SOURCE_PHI                                        C
!                                                                      C
!  Purpose: Adds point sources to the scalar equations.                C
!                                                                      C
!  Author: J. Musser                                  Date: 10-JUN-13  C
!  Reviewer:                                          Date:            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE point_source_phi(PHI, PS_PHI, PS_FLOW,  &
         m, a_m, b_m)

      use compar
      use geometry
      use indices
      use physprop
      use ps
      use run

! To be removed upon complete integration of point source routines.
      use bc
      use usr
      use functions
      use param, only: dimension_3, dimension_m
      use param1, only: one, small_number

      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------

! Vector b_m
      DOUBLE PRECISION, INTENT(IN) :: PHI(dimension_3)


! maximum term in b_m expression
      DOUBLE PRECISION, INTENT(IN) :: PS_PHI(dimension_bc)

! maximum term in b_m expression
      DOUBLE PRECISION, INTENT(IN) :: PS_FLOW(dimension_bc)

      INTEGER, intent(in) :: M

! 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)

!-----------------------------------------------
! Local Variables
!-----------------------------------------------

! Indices
      INTEGER :: IJK, I, J, K
      INTEGER :: PSV

! terms of bm expression
      DOUBLE PRECISION pSource

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

      psv_lp: do psv = 1, dimension_ps

         if(.NOT.ps_defined(psv)) cycle psv_lp
         if(abs(ps_flow(psv)) < small_number) cycle psv_lp

         do k = ps_k_b(psv), ps_k_t(psv)
         do j = ps_j_s(psv), ps_j_n(psv)
         do i = ps_i_w(psv), ps_i_e(psv)

            if(.NOT.is_on_mype_plus2layers(i,j,k)) cycle
            IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells

            ijk = funijk(i,j,k)
            if(.NOT.fluid_at(ijk)) cycle

            if(a_m(ijk,0,m) == -one .AND. b_m(ijk,m) == -phi(ijk)) then
               b_m(ijk,m) = -ps_phi(psv)
            else

! Calculate the mass flow rate for this cell. (mass/time)
               psource = ps_flow(psv) * (vol(ijk)/ps_volume(psv))

               a_m(ijk,0,m) = a_m(ijk,0,m) - psource
               b_m(ijk,m) = b_m(ijk,m) - ps_phi(psv) * psource

            endif

         enddo
         enddo
         enddo

      enddo psv_lp

      RETURN
      END SUBROUTINE point_source_phi