solve_scalar_eq.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOLVE_SCALAR_EQ                                         C
!  Purpose: Solve scalar transport equations                           C
!                                                                      C
!                                                                      C
!  Author: M. Syamlal                                 Date: 4-12-99    C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE solve_scalar_eq(IER)

! Modules
!---------------------------------------------------------------------//
      use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
      use bc, only: bc_scalar, bc_hw_scalar, bc_c_scalar, bc_scalarw
      use compar, only: ijkstart3, ijkend3
      use fldvar, only: rop_g, rop_go, ep_g, u_g, v_g, w_g
      use fldvar, only: rop_s, rop_so, ep_s, u_s, v_s, w_s
      use fldvar, only: scalar, scalaro
      use functions, only: fluid_at, zmax
      use geometry, only: vol, ijkmax2
      use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
      use mflux, only: flux_ge, flux_gn, flux_gt
      use mflux, only: flux_se, flux_sn, flux_st
      use mflux, only: flux_gse, flux_gsn, flux_gst
      use mflux, only: flux_sse, flux_ssn, flux_sst
      use param, only: dimension_3
      use param1, only: zero
      use residual, only: resid, max_resid, ijk_resid
      use residual, only: num_resid, den_resid
      use residual, only: resid_sc
      use run, only: discretize, added_mass, m_am, odt
      use rxns, only: sum_r_g, sum_r_s
      use scalars, only: nscalar, phase4scalar
      use scalars, only: scalar_c, scalar_p, dif_scalar
      use toleranc, only: zero_ep_s
      use ur_facs, only: ur_fac
      use usr_src, only: call_usr_source, calc_usr_source
      use usr_src, only: usr_scalar
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Error index
      INTEGER, INTENT(INOUT) :: IER

! Local variables
!---------------------------------------------------------------------//
! phase index
      INTEGER :: M
! species index
      INTEGER :: NN
! Indices
      INTEGER :: IJK
!
      DOUBLE PRECISION :: APO
!
! linear equation solver method and iterations
      INTEGER :: LEQM, LEQI
! temporary variables in residual computation
      DOUBLE PRECISION :: res1, mres1, num_res, den_res
      INTEGER :: ires1
! source vector: coefficient of dependent variable
! becomes part of a_m matrix; must be positive
      DOUBLE PRECISION :: S_P(dimension_3)
! source vector: constant part becomes part of b_m vector
      DOUBLE PRECISION :: S_C(dimension_3)
!
      DOUBLE PRECISION :: EPS(dimension_3)
!
      character(LEN=8) :: Vname
!---------------------------------------------------------------------//

      call lock_ambm

      resid(resid_sc,0) = zero
      num_resid(resid_sc,0) = zero
      den_resid(resid_sc,0) = zero
      max_resid(resid_sc,0) = zero
      ijk_resid(resid_sc,0) = 0

! Loop over number of scalar equations
      DO nn = 1, nscalar

! Setting M to the phase assigned to convecting the indicated scalar
! equation
         m = phase4scalar(nn)
         CALL init_ab_m (a_m, b_m, ijkmax2, m)

! Gas phase
! ---------------------------------------------------------------->>>
         IF(m == 0) THEN

            DO ijk = ijkstart3, ijkend3
               IF (fluid_at(ijk)) THEN
                  apo = rop_go(ijk)*vol(ijk)*odt
                  s_p(ijk) = apo + (zmax(sum_r_g(ijk)) + &
                                         scalar_p(ijk,nn))*vol(ijk)
                  s_c(ijk) = apo*scalaro(ijk,nn) + &
                             scalar(ijk,nn)*zmax((-sum_r_g(ijk)))*vol(ijk) + &
                             scalar_c(ijk, nn)*vol(ijk)
               ELSE
                  s_p(ijk) = zero
                  s_c(ijk) = zero
               ENDIF
            ENDDO

            IF(.NOT.added_mass) THEN
               CALL conv_dif_phi (scalar(1,nn), dif_scalar(1,nn), &
                                  discretize(9), u_g, v_g, w_g, &
                                  flux_ge, flux_gn, flux_gt, m, a_m, b_m)
            ELSE
               CALL conv_dif_phi (scalar(1,nn), dif_scalar(1,nn), &
                                  discretize(9), u_g, v_g, w_g, &
                                  flux_gse, flux_gsn, flux_gst, m, a_m, b_m)
            ENDIF

            CALL bc_phi (scalar(1,nn), bc_scalar(1,nn), &
                         bc_scalarw(1,nn), bc_hw_scalar(1,nn), &
                         bc_c_scalar(1,nn), m, a_m, b_m)

            CALL source_phi (s_p, s_c, ep_g, scalar(1,nn), m, a_m, b_m)

! calculate any usr source terms
            IF (call_usr_source(9)) CALL calc_usr_source(usr_scalar,&
                                  a_m, b_m, lm=m, ln=nn)

            CALL calc_resid_s (scalar(1,nn), a_m, b_m, m, num_res, &
                               den_res, res1, mres1, ires1, zero)
            resid(resid_sc,0) = resid(resid_sc,0)+res1
            num_resid(resid_sc,0) = num_resid(resid_sc,0)+num_res
            den_resid(resid_sc,0) = den_resid(resid_sc,0)+den_res
            if(mres1 .gt. max_resid(resid_sc,0))then
               max_resid(resid_sc,0) = mres1
               ijk_resid(resid_sc,0) = ires1
            endif

            CALL under_relax_s (scalar(1,nn), a_m, b_m, m, ur_fac(9))

!           call check_ab_m(a_m, b_m, m, .false., ier)
!           call write_ab_m(a_m, b_m, ijkmax2, m, ier)
!           write(*,*) res1, mres1, ires1
!           call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, &
!                            DO_K, ier)

            CALL adjust_leq (res1, leq_it(9), leq_method(9), &
               leqi, leqm)

            write(vname, '(A,I2)')'Scalar',nn
            CALL solve_lin_eq (vname, 9, scalar(1,nn), a_m, b_m, m,&
                               leqi, leqm, leq_sweep(9), leq_tol(9),&
                               leq_pc(9), ier)
!           call out_array(Scalar(1, N), Vname)

! Solids phase
! ---------------------------------------------------------------->>>
         ELSE   ! (m/=0) , i.e., solids phase

            DO ijk = ijkstart3, ijkend3
               IF (fluid_at(ijk)) THEN
                  apo = rop_so(ijk, m)*vol(ijk)*odt
                  s_p(ijk) = apo + (zmax(sum_r_s(ijk, m)) + &
                                         scalar_p(ijk, nn))*vol(ijk)
                  s_c(ijk) = apo*scalaro(ijk,nn) + &
                             scalar(ijk,nn)*zmax((-sum_r_s(ijk, m)))*vol(ijk)+&
                             scalar_c(ijk, nn)*vol(ijk)
                  eps(ijk) = ep_s(ijk, m)
               ELSE
                  s_p(ijk) = zero
                  s_c(ijk) = zero
                  eps(ijk) = zero
               ENDIF
            ENDDO

            IF(.NOT.added_mass .OR. m /= m_am) THEN
               CALL conv_dif_phi (scalar(1,nn), dif_scalar(1,nn), &
                                  discretize(9), u_s(1,m), v_s(1,m), &
                                  w_s(1,m), flux_se(1,m), flux_sn(1,m),&
                                  flux_st(1,m), m, a_m, b_m)
            ELSE ! virtual mass term added for M = M_AM ONLY!!!!
               CALL conv_dif_phi (scalar(1,nn), dif_scalar(1,nn), &
                                  discretize(9), u_s(1,m), v_s(1,m), &
                                  w_s(1,m), flux_sse, flux_ssn, flux_sst, &
                                  m, a_m, b_m)
            ENDIF

            CALL bc_phi (scalar(1,nn), bc_scalar(1,nn), bc_scalarw(1,nn), &
                         bc_hw_scalar(1,nn), bc_c_scalar(1,nn), m, a_m, b_m)

            CALL source_phi (s_p, s_c, eps, scalar(1,nn), m, a_m, b_m)

! calculate any usr source terms
            IF (call_usr_source(9)) CALL calc_usr_source(usr_scalar,&
                                  a_m, b_m, lm=m, ln=nn)


            CALL calc_resid_s (scalar(1,nn), a_m, b_m, m, num_res, &
                              den_res, res1, mres1, ires1, zero)
            resid(resid_sc,0) = resid(resid_sc,0)+res1
            num_resid(resid_sc,0) = num_resid(resid_sc,0)+num_res
            den_resid(resid_sc,0) = den_resid(resid_sc,0)+den_res
            if(mres1 .gt. max_resid(resid_sc,0))then
               max_resid(resid_sc,0) = mres1
               ijk_resid(resid_sc,0) = ires1
            endif

            CALL under_relax_s (scalar(1,nn), a_m, b_m, m, ur_fac(9))

!            call check_ab_m(a_m, b_m, m, .false., ier)
!            call write_ab_m(a_m, b_m, ijkmax2, m, ier)
!            write(*,*) res1, mres1, ires1
!            call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, &
!                             DO_K, ier)

            CALL adjust_leq (res1, leq_it(9), leq_method(9), &
                             leqi, leqm)

            write(vname, '(A,I2)')'Scalar',nn
            CALL solve_lin_eq (vname, 9, scalar(1,nn), a_m, b_m, m, &
                               leqi, leqm, leq_sweep(9), leq_tol(9), &
                               leq_pc(9), ier)
!            call out_array(Scalar(1, N), Vname)

         ENDIF
      ENDDO

      call unlock_ambm

      RETURN
      END SUBROUTINE solve_scalar_eq