conv_rop_s.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CONV_ROP_s                                              C
!  Purpose: Determine convection terms for solids continuity           C
!           equation.  Master routine.                                 C
!                                                                      C
!  Author: M. Syamlal                                 Date: 18-MAR-97  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE conv_rop_s(A_M, B_M, M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE fldvar
      USE run
      USE compar
      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)
! Solids phase index
      INTEGER, INTENT(IN) :: M
!-----------------------------------------------

      IF (discretize(2) == 0) THEN               ! 0 & 1 => first order upwinding
         CALL conv_rop_s0 (a_m, m)
      ELSE
         CALL conv_rop_s1 (a_m, m)
      ENDIF

      RETURN
      END SUBROUTINE conv_rop_s


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CONV_ROP_s0                                             C
!  Purpose: Determine convection terms for solids continuity equation. C
!                                                                      C
!  Notes: The off-diagonal coefficients calculated here must be        C
!         positive. The center coefficient and the source vector       C
!         are negative. -- First order upwinding                       C
!                                                                      C
!  Author: M. Syamlal                                 Date: 3-JUL-96   C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE conv_rop_s0(A_M, M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE fldvar
      USE run
      USE parallel
      USE physprop
      USE geometry
      USE indices
      USE pgcor
      USE pscor
      USE compar
      USE functions
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
! Solids phase index
      INTEGER, INTENT(IN) :: M
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP
      INTEGER :: IMJK, IJMK, IJKM
!-----------------------------------------------

! Calculate convection-diffusion fluxes through each of the faces

!!$omp  parallel do private( I, J, K, IJK, IPJK, IJPK, IJKP,  &
!!$omp&  IMJK, IJMK, IJKM) &
!!$omp&  schedule(static)
      DO ijk = ijkstart3, ijkend3
         IF (phase_4_p_g(ijk)/=m .AND. phase_4_p_s(ijk)/=m) THEN
! if either phase_4_p_g or phase_4_p_s are assigned the index of the
! current solids phase then this section is skipped. currently,
! phase_4_p_g and phase_4_p_s are always assigned to 0 and undefined,
! respectively.  Hence this section should always be entered.

! it was previously possible for phase_4_p_g to be assigned the index of
! the gas phase in some cells, and the index of a solids phase in other
! cells if that solids phase has close_packed=F and was in higher
! concentrations than the gas. in such a case this branch would become
! skipped, while the corresponding gas continuity would become
! activated. moreover, that solids phase's volume fraction would be
! corrected rather than the gas phases void fraction in calc_vol_fr.

! if phase_4_p_s were to be assigned the index of a solids phase that
! can close pack if the current cell exhibited close packing, then
! this branch would become skipped.  moreover, that solids phase's
! volume fraction would then be corected based on the value of
! maximum packing in that cell and the sum of all other solids that
! can close pack in calc_vol_fr.
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)
            ipjk = ip_of(ijk)
            ijpk = jp_of(ijk)
            ijkp = kp_of(ijk)
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ijkm = km_of(ijk)

! East face (i+1/2, j, k)
            a_m(ijk,east,m) = zmax((-u_s(ijk,m)))*ayz(ijk)
            a_m(ipjk,west,m) = zmax(u_s(ijk,m))*ayz(ijk)

! North face (i, j+1/2, k)
            a_m(ijk,north,m) = zmax((-v_s(ijk,m)))*axz(ijk)
            a_m(ijpk,south,m) = zmax(v_s(ijk,m))*axz(ijk)

! Top face (i, j, k+1/2)
            IF (do_k) THEN
               a_m(ijk,top,m) = zmax((-w_s(ijk,m)))*axy(ijk)
               a_m(ijkp,bottom,m) = zmax(w_s(ijk,m))*axy(ijk)
            ENDIF

! Modify west (i-1/2,j,k), south (i j-1/2,k) and bottom (i,j,k-1/2)
! faces if the neighboring west, south, bottom cells have
! phase_4_p_g of m or phase_4_p_s of m.
            IF (phase_4_p_g(imjk)==m .OR. phase_4_p_s(imjk)==m) &
               a_m(ijk,west,m) = zmax(u_s(imjk,m))*ayz(imjk)
            IF (phase_4_p_g(ijmk)==m .OR. phase_4_p_s(ijmk)==m) &
               a_m(ijk,south,m) = zmax(v_s(ijmk,m))*axz(ijmk)

            IF (do_k) THEN
               IF (phase_4_p_g(ijkm)==m .OR. phase_4_p_s(ijkm)==m) &
                  a_m(ijk,bottom,m) = zmax(w_s(ijkm,m))*axy(ijkm)
            ENDIF
         ENDIF
      END DO

      RETURN
      END SUBROUTINE conv_rop_s0


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CONV_ROP_s1                                             C
!  Purpose: Determine convection terms for solids continuity equation. C
!  Notes: The off-diagonal coefficients calculated here must be        C
!         positive. The center coefficient and the source vector       C
!         are negative. -- Higher order methods                        C
!         See conv_rop_s0 for additional details.                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 18-MAR-97  C
!  Reviewer:                                          Date:            C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!                                                                      C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE conv_rop_s1(A_M, M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE compar
      USE fldvar
      USE functions
      USE geometry
      USE indices
      USE parallel
      USE param
      USE param1
      USE pgcor
      USE physprop
      USE pscor
      USE run
      USE xsi
      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! Septadiagonal matrix A_m
      DOUBLE PRECISION, INTENT(INOUT) :: A_m(dimension_3, -3:3, 0:dimension_m)
! Solids phase index
      INTEGER, INTENT(IN) :: M
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! Indices
      INTEGER :: I, J, K, IJK, IPJK, IJPK, IJKP
      INTEGER :: IMJK, IJMK, IJKM
! loezos
      INTEGER :: incr
! temporary use of global arrays:
! xsi_array: convection weighting factors
      DOUBLE PRECISION :: XSI_e(dimension_3), XSI_n(dimension_3),&
                          XSI_t(dimension_3)
!-----------------------------------------------

! Loezos:
      incr=0

! Calculate convection factors
      CALL calc_xsi(discretize(2), rop_s(1,m), u_s(1,m), v_s(1,m),&
         w_s(1,m), xsi_e, xsi_n, xsi_t, incr)

! Calculate convection-diffusion fluxes through each of the faces

!!$omp  parallel do private( I, J, K, IJK, IPJK, IJPK, IJKP,  &
!!$omp&  IMJK, IJMK, IJKM) &
!!$omp&  schedule(static)
      DO ijk = ijkstart3, ijkend3
         IF (phase_4_p_g(ijk)/=m .AND. phase_4_p_s(ijk)/=m) THEN
            i = i_of(ijk)
            j = j_of(ijk)
            k = k_of(ijk)
            ipjk = ip_of(ijk)
            ijpk = jp_of(ijk)
            ijkp = kp_of(ijk)
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ijkm = km_of(ijk)

! East face (i+1/2, j, k)
            a_m(ijk,east,m) = -xsi_e(ijk)*u_s(ijk,m)*ayz(ijk)
            a_m(ipjk,west,m) = (one - xsi_e(ijk))*u_s(ijk,m)*ayz(ijk)

! North face (i, j+1/2, k)
            a_m(ijk,north,m) = -xsi_n(ijk)*v_s(ijk,m)*axz(ijk)
            a_m(ijpk,south,m) = (one - xsi_n(ijk))*v_s(ijk,m)*axz(ijk)

! Top face (i, j, k+1/2)
            IF (do_k) THEN
               a_m(ijk,top,m) = -xsi_t(ijk)*w_s(ijk,m)*axy(ijk)
               a_m(ijkp,bottom,m) = (one - xsi_t(ijk))*w_s(ijk,m)*axy(ijk)
            ENDIF

            IF (phase_4_p_g(imjk)==m .OR. phase_4_p_s(imjk)==m) a_m(ijk,west,m) = &
               (one - xsi_e(imjk))*u_s(imjk,m)*ayz(imjk)
            IF (phase_4_p_g(ijmk)==m .OR. phase_4_p_s(ijmk)==m) a_m(ijk,south,m) = &
               (one - xsi_n(ijmk))*v_s(ijmk,m)*axz(ijmk)
            IF (do_k) THEN
               IF (phase_4_p_g(ijkm)==m .OR. phase_4_p_s(ijkm)==m) a_m(ijk,bottom,m)&
                   = (one - xsi_t(ijkm))*w_s(ijkm,m)*axy(ijkm)
            ENDIF
         ENDIF

      ENDDO

      RETURN
      END SUBROUTINE conv_rop_s1