usr_src_mod.f

Go to the documentation of this file.

! or put the actual routine calc_usr_source here? move logic variable
!  usr_source into run or something...
      module usr_src
      use param, only: dim_eqs

! If .TRUE. call user-defined source subroutines
      LOGICAL :: call_usr_source(dim_eqs)

      enum, bind(c)
         enumerator :: pressure_correction, solids_correction
         enumerator :: gas_continuity, solids_continuity
         enumerator :: gas_u_mom, solids_u_mom, gas_v_mom, solids_v_mom
         enumerator :: gas_w_mom, solids_w_mom
         enumerator :: gas_energy, solids_energy
         enumerator :: gas_species, solids_species
         enumerator :: gran_energy
         enumerator :: usr_scalar, k_epsilon_k, k_epsilon_e
         enumerator :: blank
         enumerator :: dummy
      END enum

      CONTAINS
!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: CALC_USR_SOURCE                                         C
!  Purpose: Driver routine to calculate user defined source terms      C
!  for each of the various equations                                   C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE calc_usr_source(lEQ_NO, A_M, B_M, lB_MMAX, lM, lN)

! Modules
!-----------------------------------------------
      use compar, only: ijkstart3, ijkend3
      use fldvar, only: ep_s, ep_g
      use fun_avg, only: avg_x, avg_y, avg_z
      use functions, only: fluid_at, wall_at
      use functions, only: ip_at_e, ip_at_n, ip_at_t
      use functions, only: sip_at_e, sip_at_n, sip_at_t
      use functions, only: east_of, north_of, top_of
      use indices, only: i_of, j_of, k_of
      use param, only: dimension_3, dimension_m
      use param1, only: undefined_i, zero
      use pgcor, only: phase_4_p_g
      use physprop, only: smax, mmax
      use pscor, only: phase_4_p_s
      use run, only: kt_type_enum, ghd_2007
      use run, only: momentum_x_eq, momentum_y_eq, momentum_z_eq
      use toleranc, only: dil_ep_s
      use error_manager
      IMPLICIT NONE
! Dummy arguments
!-----------------------------------------------
! Equation number
! Table relating the eq_no set in the mfix.dat to the corresponding
! pde/variable of interest
!     1, PP_g (pressure correction eq)
!     2, EPP (solids correction eq)
!     2, rop_g (gas continuity eq)
!     2, rop_s (solids continuity eq)
!     3, u_g (gas u-momentum eq)
!     3, u_s (solids u-momentum eq)
!     4, v_g (gas v-momentum eq)
!     4, v_s (solids v-momentum eq)
!     5, w_g (gas w-momentum eq)
!     5, w_s (solids w-momentum eq)
!     6, T_g (gas energy eq)
!     6, T_s (solids energy eq)
!     7, x_g (gas species eq)
!     7, x_s (solids species eq)
!     8, theta_m (granular energy eq)
!     9, scalar_eq (scalar eq)
!     9, k_epsilon k (k_epsilon eqs)
!     9, k_epsilon e (k_epsilon_eqs)
      INTEGER, INTENT(IN) :: lEQ_NO

! 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)
! Vector b_mmax
      DOUBLE PRECISION, OPTIONAL, INTENT(INOUT) :: lB_mmax(dimension_3, 0:dimension_m)
! Phase index
      INTEGER, OPTIONAL, INTENT(IN) :: lM
! Species index OR scalar equation number
      INTEGER, OPTIONAL, INTENT(IN) :: lN

! Local variables
!-----------------------------------------------
! source terms which appear appear in the
! center coefficient (lhs) - part of a_m matrix
! source vector (rhs) - part of b_m vector
      DOUBLE PRECISION :: sourcelhs, sourcerhs
! indices
      INTEGER :: IJK, I, J, K, IJKE, IJKN, IJKT
! tmp indices
      INTEGER :: L, ll, M, N
! volume fractions
      DOUBLE PRECISION :: EPGA, EPSA, EPtmp
!-----------------------------------------------

! set local values for phase index and species or scalar index
      IF (.NOT.present(lm)) THEN
         m = undefined_i
      ELSE
         m = lm
      ENDIF
      IF (.NOT.present(ln)) THEN
         n = undefined_i
      ELSE
         n = ln
      ENDIF

      SELECT CASE(leq_no)

! gas pressure correction equation
      CASE(pressure_correction)
         DO ijk=ijkstart3,ijkend3
            IF (fluid_at(ijk)) THEN
               CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, m, n)
               a_m(ijk,0,m) = a_m(ijk,0,m) - sourcelhs
               b_m(ijk,m) = b_m(ijk,m) - sourcerhs
               lb_mmax(ijk,m) = max(abs(lb_mmax(ijk,m)), abs(b_m(ijk,m)))
            ENDIF
         ENDDO

! solids correction currently only called when smax==1 and
! mcp/=undefined_i
      CASE(solids_correction)
         DO ijk=ijkstart3,ijkend3
            IF (fluid_at(ijk)) THEN
               CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, m, n)
               a_m(ijk,0,0) = a_m(ijk,0,0) - sourcelhs
               b_m(ijk,0) = b_m(ijk,0) - sourcerhs
               lb_mmax(ijk,0) = max(abs(lb_mmax(ijk,0)), abs(b_m(ijk,0)))
            ENDIF
         ENDDO

! gas continuity
      CASE(gas_continuity)
         DO ijk=ijkstart3,ijkend3
            IF (fluid_at(ijk)) THEN
!                            .AND. PHASE_4_P_G(IJK)/=M) THEN
! have not included phase_4_p_g logic...which would require unique
! structure here
               CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, m, n)
               a_m(ijk,0,m) = a_m(ijk,0,m) - sourcelhs
               b_m(ijk,m) = b_m(ijk,m) - sourcerhs
            ENDIF
         ENDDO

! soldis continuity
      CASE(solids_continuity)
! have not included phase_4_p_g/p_s logic...which would require unique
! structure here
         DO ijk=ijkstart3,ijkend3
            IF (fluid_at(ijk)) THEN
!               .AND. PHASE_4_P_G(IJK)/=M .AND. PHASE_4_P_S(IJK)/=M) THEN
               CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, m, n)
               a_m(ijk,0,m) = a_m(ijk,0,m) - sourcelhs
               b_m(ijk,m) = b_m(ijk,m) - sourcerhs
            ENDIF
         ENDDO

! u-momentum equations
      CASE(gas_u_mom)
         m=0
         IF(.NOT.momentum_x_eq(m)) RETURN
         DO ijk=ijkstart3,ijkend3
            IF (.NOT.fluid_at(ijk)) cycle
            i = i_of(ijk)
            ijke = east_of(ijk)
            epga = avg_x(ep_g(ijk),ep_g(ijke),i)
            IF (ip_at_e(ijk) .OR. epga <= dil_ep_s) cycle
            CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, m, n)
            a_m(ijk,0,m) = a_m(ijk,0,m) - sourcelhs
            b_m(ijk,m) = b_m(ijk,m) - sourcerhs
         ENDDO

! u-momentum equations
      CASE(solids_u_mom)
         DO l=1,mmax
            IF (kt_type_enum /= ghd_2007 .OR. &
               (kt_type_enum == ghd_2007 .AND. l==mmax)) THEN
               IF(.NOT.momentum_x_eq(l)) RETURN
               DO ijk=ijkstart3,ijkend3
                  IF (.NOT.fluid_at(ijk)) cycle
                  i = i_of(ijk)
                  ijke = east_of(ijk)
                  eptmp = zero
                  IF (kt_type_enum == ghd_2007) THEN
                     DO ll = 1, smax
                        eptmp = eptmp + avg_x(ep_s(ijk,ll),ep_s(ijke,ll),i)
                     ENDDO
                     epsa = eptmp
                  ELSE
                     epsa = avg_x(ep_s(ijk,l),ep_s(ijke,l),i)
                  ENDIF
                  IF (ip_at_e(ijk) .OR. sip_at_e(ijk) .OR. &
                      epsa <= dil_ep_s) cycle
                  CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, l, n)
                  a_m(ijk,0,l) = a_m(ijk,0,l) - sourcelhs
                  b_m(ijk,l) = b_m(ijk,l) - sourcerhs
               ENDDO    ! enddo ijk
            ENDIF
         ENDDO   ! enddo mmax

! v-momentum equations
      CASE(gas_v_mom)
         m=0
         IF(.NOT.momentum_y_eq(m)) RETURN
         DO ijk=ijkstart3,ijkend3
            IF (.NOT.fluid_at(ijk)) cycle
            j = j_of(ijk)
            ijkn = north_of(ijk)
            epga = avg_y(ep_g(ijk),ep_g(ijkn),j)
            IF (ip_at_n(ijk) .OR. epga <= dil_ep_s) cycle
            CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, m, n)
            a_m(ijk,0,m) = a_m(ijk,0,m) - sourcelhs
            b_m(ijk,m) = b_m(ijk,m) - sourcerhs
         ENDDO

! v-momentum equations
      CASE(solids_v_mom)
         DO l=1,mmax
            IF (kt_type_enum /= ghd_2007 .OR. &
                (kt_type_enum == ghd_2007 .AND. l==mmax)) THEN
               IF(.NOT.momentum_y_eq(l)) RETURN
               DO ijk=ijkstart3,ijkend3
                  IF (.NOT.fluid_at(ijk)) cycle
                  j = j_of(ijk)
                  ijkn = north_of(ijk)
                  eptmp = zero
                  IF (kt_type_enum == ghd_2007) THEN
                     DO ll = 1, smax
                        eptmp = eptmp + avg_y(ep_s(ijk,ll),ep_s(ijkn,ll),j)
                     ENDDO
                     epsa = eptmp
                  ELSE
                     epsa = avg_y(ep_s(ijk,l),ep_s(ijkn,l),j)
                  ENDIF
                  IF (ip_at_n(ijk) .OR. sip_at_n(ijk) .OR. &
                      epsa <= dil_ep_s) cycle
                  CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, l, n)
                  a_m(ijk,0,l) = a_m(ijk,0,l) - sourcelhs
                  b_m(ijk,l) = b_m(ijk,l) - sourcerhs
               ENDDO   ! enddo ijk
            ENDIF
         ENDDO   ! enddo mmax

! w-momentum equations
      CASE (gas_w_mom)
         m=0
         IF(.NOT.momentum_z_eq(m)) RETURN
         DO ijk=ijkstart3,ijkend3
            IF (.NOT.fluid_at(ijk)) cycle
            k = k_of(ijk)
            ijkt = top_of(ijk)
            epga = avg_z(ep_g(ijk),ep_g(ijkt),k)
            IF (ip_at_t(ijk) .OR. epga <= dil_ep_s) cycle
            CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, m, n)
            a_m(ijk,0,m) = a_m(ijk,0,m) - sourcelhs
            b_m(ijk,m) = b_m(ijk,m) - sourcerhs
         ENDDO

! w-momentum equations
      CASE (solids_w_mom)
         DO l=1,mmax
            IF (kt_type_enum /= ghd_2007 .OR. &
                (kt_type_enum == ghd_2007 .AND. l==mmax)) THEN
               IF(.NOT.momentum_z_eq(l)) RETURN
               DO ijk=ijkstart3,ijkend3
                  IF (.NOT.fluid_at(ijk)) cycle
                  k = k_of(ijk)
                  ijkt = top_of(ijk)
                  eptmp = zero
                  IF (kt_type_enum == ghd_2007) THEN
                     DO ll = 1, smax
                        eptmp = eptmp + avg_z(ep_s(ijk,ll),ep_s(ijkt,ll),k)
                     ENDDO
                     epsa = eptmp
                  ELSE
                     epsa = avg_z(ep_s(ijk,l),ep_s(ijkt,l),k)
                  ENDIF
                  IF (ip_at_t(ijk) .OR. sip_at_t(ijk) .OR. &
                      epsa <= dil_ep_s) cycle
                  CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, l, n)
                  a_m(ijk,0,l) = a_m(ijk,0,l) - sourcelhs
                  b_m(ijk,l) = b_m(ijk,l) - sourcerhs
               ENDDO   ! enddo ijk
            ENDIF
         ENDDO   ! enddo mmax


! gas and solids energy equations (6)
! gas and solids species mass fractions (7)
! scalars, k_epsilon k and e (9)
      CASE (gas_energy, solids_energy, gas_species, solids_species,&
            usr_scalar, k_epsilon_k, k_epsilon_e )
          DO ijk = ijkstart3, ijkend3
            IF (fluid_at(ijk)) THEN
               IF (m==0) THEN
                  eptmp=ep_g(ijk)
               ELSE
                  eptmp=ep_s(ijk,m)
               ENDIF
               IF (eptmp <= dil_ep_s) cycle
               CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, m, n)
               a_m(ijk,0,m) = a_m(ijk,0,m) - sourcelhs
               b_m(ijk,m) = b_m(ijk,m) - sourcerhs
            ENDIF
         ENDDO


      CASE (gran_energy)   ! unique due to ghd
! granular temperature
         DO ijk = ijkstart3, ijkend3
            IF (fluid_at(ijk)) THEN
               IF (kt_type_enum == ghd_2007) THEN
                  eptmp = zero
                  DO ll=1,smax
                     eptmp=eptmp+ep_s(ijk,ll)
                  ENDDO
               ELSE
                  eptmp = ep_s(ijk,m)
               ENDIF
            ENDIF
            IF (eptmp<=dil_ep_s) cycle
            CALL usr_sources(leq_no, ijk, sourcelhs, sourcerhs, m, n)
            a_m(ijk,0,m) = a_m(ijk,0,m) - sourcelhs
            b_m(ijk,m) = b_m(ijk,m) - sourcerhs
         ENDDO


! error out
      CASE DEFAULT
! should never hit this
! Initialize the error manager.
         CALL init_err_msg("CALC_USR_SOURCE")
         WRITE(err_msg, 1001) ival(leq_no)
         CALL flush_err_msg(abort=.true.)
 1001 FORMAT('Error 1101: Unknown Equation= ', a)
         CALL finl_err_msg
      END SELECT   ! end selection of usr_source equation

      RETURN
      END SUBROUTINE calc_usr_source

      end module usr_src