source_rop_s.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOURCE_ROP_s                                            C
!  Purpose: Determine source terms for solids continuity equation.     C
!                                                                      C
!  Notes: The off-diagonal coefficients are positive. The center       C
!         coefficient and the source vector are negative.              C
!         See conv_rop_s0 for additional details.                      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 source_rop_s(A_M, B_M, M)

!-----------------------------------------------
! Modules
!-----------------------------------------------
      USE param
      USE param1
      USE parallel
      USE fldvar
      USE rxns
      USE run
      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)
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
! solids phase index
      INTEGER, INTENT(IN) :: M
!-----------------------------------------------
! Local variables
!-----------------------------------------------
! DEL dot V
      DOUBLE PRECISION :: DEL_V
! Mass source
      DOUBLE PRECISION :: Src
! Indices
      INTEGER :: I, J, K, IJK, IMJK, IJMK, IJKM
! error message
      CHARACTER(LEN=80) :: LINE
!-----------------------------------------------

!!$omp  parallel do private( I, J, K, IJK, IMJK, IJMK, IJKM,  DEL_V, &
!!$omp&  Src, LINE) &
!!$omp&  schedule(static)
      DO ijk = ijkstart3, ijkend3
         IF (fluid_at(ijk) .AND. 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)
            imjk = im_of(ijk)
            ijmk = jm_of(ijk)
            ijkm = km_of(ijk)

            del_v = u_s(ijk,m)*ayz(ijk) - u_s(imjk,m)*ayz(imjk) +&
                    v_s(ijk,m)*axz(ijk) - v_s(ijmk,m)*axz(ijmk) + &
                    w_s(ijk,m)*axy(ijk) - w_s(ijkm,m)*axy(ijkm)

            IF (rop_s(ijk,m) > zero) THEN
               src = vol(ijk)*zmax((-sum_r_s(ijk,m)))/rop_s(ijk,m)
            ELSE
               src = zero
            ENDIF

            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)+&
                             vol(ijk)*odt+zmax(del_v)+src)
            b_m(ijk,m) = -(rop_so(ijk,m)*vol(ijk)*odt+&
                           zmax((-del_v))*rop_s(ijk,m)+&
                           zmax(sum_r_s(ijk,m))*vol(ijk))

            IF (abs(a_m(ijk,0,m)) < small_number) THEN
               IF (abs(b_m(ijk,m)) < small_number) THEN
                  a_m(ijk,0,m) = -one            ! Equation is undefined.
                  b_m(ijk,m) = -rop_s(ijk,m)     ! Use existing value
               ELSE
!!$omp             critical
                  WRITE (line, '(A,I6,A,I1,A,G12.5)') 'Error: At IJK = ', ijk, &
                     ' M = ', m, ' A = 0 and b = ', b_m(ijk,m)
                  CALL write_error ('SOURCE_ROP_s', line, 1)
!!$omp             end critical
               ENDIF
            ENDIF
         ELSE
! set the value of rop_s in all wall and flow boundary cells to what is
! known for that cell
            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) = -rop_s(ijk,m)
         ENDIF
      ENDDO

      RETURN
      END SUBROUTINE source_rop_s

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: POINT_SOURCE_ROP_S                                      C
!  Purpose: Adds point sources to the solids continuity equation.      C
!                                                                      C
!                                                                      C
!  Author: J. Musser                                  Date: 10-JUN-13  C
!  Reviewer:                                          Date:            C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE point_source_rop_s(B_M, M)

      use compar
      use constant
      use geometry
      use indices
      use param, only: dimension_m, dimension_3
      use param1, only: small_number
      use physprop
      use ps
      use run
      use functions

      IMPLICIT NONE
!-----------------------------------------------
! Dummy arguments
!-----------------------------------------------
! Vector b_m
      DOUBLE PRECISION, INTENT(INOUT) :: B_m(dimension_3, 0:dimension_m)
! Solids phase index.
      INTEGER, INTENT(IN) :: M

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

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

! terms of bm expression
      DOUBLE PRECISION :: pSource

!-----------------------------------------------
      ps_lp: do psv = 1, dimension_ps

         if(.NOT.ps_defined(psv)) cycle ps_lp
         if(ps_massflow_s(psv,m) < small_number) cycle ps_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(fluid_at(ijk)) then

               psource =  ps_massflow_s(psv,m) * &
                  (vol(ijk)/ps_volume(psv))

               b_m(ijk,m) = b_m(ijk,m) - psource
            endif
         enddo
         enddo
         enddo

      enddo ps_lp

      RETURN
      END SUBROUTINE point_source_rop_s