solve_species_eq.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SOLVE_SPECIES_EQ                                        C
!  Purpose: Solve species mass balance equations in matrix equation    C
!     form Ax=b. The center coefficient (ap) and source vector (b)     C
!     are negative.  The off-diagonal coefficients are positive.       C
!                                                                      C
!  Author: M. Syamlal                                 Date: 11-FEB-98  C
!  Reviewer:                                          Date:            C
!                                                                      C
!  Literature/Document References:                                     C
!                                                                      C
!  Variables referenced:                                               C
!  Variables modified:                                                 C
!  Local variables:                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE solve_species_eq(IER)

! Modules
!---------------------------------------------------------------------//
      use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
      use bc, only: bc_x_g, bc_xw_g, bc_hw_x_g, bc_c_x_g
      use bc, only: bc_x_s, bc_xw_s, bc_hw_x_s, bc_c_x_s
      use chischeme, only: set_chi, unset_chi
      use compar, only: ijkstart3, ijkend3, mype, numpes
      use fldvar, only: ep_g, u_g, v_g, w_g, x_g, x_go, rop_go
      use fldvar, only: ep_s, u_s, v_s, w_s, x_s, x_so, rop_so
      use functions, only: fluid_at, zmax
      use geometry, only: ijkmax2, vol
      use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
      use mflux, only: flux_ge, flux_gse, flux_gn, flux_gsn
      use mflux, only: flux_se, flux_sse, flux_sn, flux_ssn
      use mflux, only: flux_gt, flux_gst, flux_st, flux_sst
      use mpi_utility, only: global_all_sum
      use param, only: dimension_3, dimension_m
      use param, only: dimension_n_s, dimension_n_g
      use param1, only: zero
      use physprop, only: nmax, dif_g, dif_s
      use physprop, only: smax
      use ps, only: point_source
      use ps, only: ps_x_g, ps_massflow_g
      use ps, only: ps_x_s, ps_massflow_s
      use residual, only: resid, max_resid, ijk_resid
      use residual, only: num_resid, den_resid
      use residual, only: resid_x
      use run, only: species_eq, discretize, odt, added_mass, m_am
      use run, only: chi_scheme
      use rxns, only: sum_r_g, rox_gc, r_gp
      use rxns, only: sum_r_s, rox_sc, r_sp
      use toleranc, only: zero_x_gs
      use ur_facs, only: ur_fac
      use usr_src, only: call_usr_source, calc_usr_source
      use usr_src, only: gas_species, solids_species
      use utilities, only: bound_x
      IMPLICIT NONE

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

! Local variables
!---------------------------------------------------------------------//
! phase index
      INTEGER :: M
! species index
      INTEGER :: LN
! previous time step term
      DOUBLE PRECISION :: apo
! Indices
      INTEGER :: IJK
! linear equation solver method and iterations
      INTEGER :: LEQM, LEQI
! tmp array to pass to set_chi
      DOUBLE PRECISION :: X_s_temp(dimension_3, dimension_n_s)

! Arrays for storing errors:
! 130 - Gas phase species equation diverged
! 131 - Solids phase species equation diverged
! 13x - Unclassified
      INTEGER :: Err_l(0:numpes-1)  ! local
      INTEGER :: Err_g(0:numpes-1)  ! global

! temporary use of global arrays:
! array1 (locally s_p)
! source lhs: coefficient of dependent variable
! becomes part of a_m matrix; must be positive
      DOUBLE PRECISION :: S_P(dimension_3)
! array2 (locally s_c)
! source rhs vector: constant part becomes part of b_m vector
      DOUBLE PRECISION :: S_C(dimension_3)
! array3 (locally eps)
! alias for solids volume fraction
      DOUBLE PRECISION :: eps(dimension_3)
! array4 (locally vxgama)
      DOUBLE PRECISION :: vxgama(dimension_3)
! Septadiagonal matrix A_m, vector b_m
!      DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
!      DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)

! External functions
!---------------------------------------------------------------------//
      DOUBLE PRECISION , EXTERNAL :: Check_conservation
!---------------------------------------------------------------------//

      call lock_ambm       ! locks arrys a_m and b_m

! Initialize error flag.
      err_l = 0

! Fluid phase species mass balance equations
! ---------------------------------------------------------------->>>
      IF (species_eq(0)) THEN
         IF(chi_scheme) call set_chi(discretize(7), x_g, nmax(0), &
                                     u_g, v_g, w_g)

! looping over species
         DO ln = 1, nmax(0)
            CALL init_ab_m (a_m, b_m, ijkmax2, 0)
!!$omp    parallel do private(IJK, APO)
            DO ijk = ijkstart3, ijkend3
               IF (fluid_at(ijk)) THEN
! calculate the source terms to be used in the a matrix and b vector
                   apo = rop_go(ijk)*vol(ijk)*odt
                   s_p(ijk) = apo+&
                      (zmax(sum_r_g(ijk))+rox_gc(ijk,ln))*vol(ijk)
                   s_c(ijk) = apo*x_go(ijk,ln) + &
                      x_g(ijk,ln)*zmax((-sum_r_g(ijk)))*vol(ijk) +&
                      r_gp(ijk,ln)*vol(ijk)
               ELSE
                  s_p(ijk) = zero
                  s_c(ijk) = zero
               ENDIF
            ENDDO

! calculate the convection-diffusion terms
            IF(.NOT.added_mass) THEN
               CALL conv_dif_phi (x_g(1,ln), dif_g(1,ln), &
                  discretize(7), u_g, v_g, w_g, &
                  flux_ge, flux_gn, flux_gt, 0, a_m, b_m)
            ELSE
               CALL conv_dif_phi (x_g(1,ln), dif_g(1,ln),&
                  discretize(7), u_g, v_g, w_g, &
                  flux_gse, flux_gsn, flux_gst, 0, a_m, b_m)
            ENDIF

! calculate standard bc
            CALL bc_phi (x_g(1,ln), bc_x_g(1,ln), bc_xw_g(1,ln), &
               bc_hw_x_g(1,ln), bc_c_x_g(1,ln), 0, a_m, b_m)

! set the source terms in a and b matrix form
            CALL source_phi (s_p, s_c, ep_g, x_g(1,ln), 0, a_m, b_m)

! Add point sources.
            IF(point_source) CALL point_source_phi (x_g(1,ln), &
               ps_x_g(:,ln), ps_massflow_g, 0, a_m, b_m)

! usr sources
            IF(call_usr_source(8)) CALL calc_usr_source(gas_species, &
                                 a_m, b_m, lm=0, ln=ln)

            CALL calc_resid_s (x_g(1,ln), a_m, b_m, 0, &
               num_resid(resid_x+(ln-1),0), &
               den_resid(resid_x+(ln-1),0), resid(resid_x+(ln-1),0), &
               max_resid(resid_x+(ln-1),0), ijk_resid(resid_x+(ln-1),0), &
               zero_x_gs)

            CALL under_relax_s (x_g(1,ln), a_m, b_m, 0, ur_fac(7))

            CALL adjust_leq (resid(resid_x+(ln-1),0), leq_it(7), &
               leq_method(7), leqi, leqm)

! solve the phi equation
            CALL solve_lin_eq ('X_g', 7, x_g(1,ln), a_m, b_m, 0, &
               leqi, leqm, leq_sweep(7), leq_tol(7), leq_pc(7), ier)

! Check for linear solver divergence.
               IF(ier == -2) err_l(mype) = 130

            CALL bound_x (x_g(1,ln), ijkmax2)

         ENDDO    ! end do loop (ln = 1, nmax(0)
         IF(chi_scheme) call unset_chi()
      ENDIF
! end fluid phase species equations
! ----------------------------------------------------------------<<<

! Solids phase species balance equations
! ---------------------------------------------------------------->>>
      DO m = 1, smax
         IF (species_eq(m)) THEN
            IF(chi_scheme) THEN
               DO ln = 1, nmax(m)
                  DO ijk = ijkstart3, ijkend3
                     x_s_temp(ijk, ln) = x_s(ijk,m,ln)
                  ENDDO
                ENDDO
              call set_chi(discretize(7), x_s_temp, nmax(m), &
                           u_s(1,m), v_s(1,m), w_s(1,m))
            ENDIF ! for chi_scheme

            DO ln = 1, nmax(m)
               CALL init_ab_m (a_m, b_m, ijkmax2, m)

!!$omp    parallel do private(IJK, APO)
               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))+rox_sc(ijk,m,ln))*vol(ijk)
                    s_c(ijk) = apo*x_so(ijk,m,ln) + &
                       x_s(ijk,m,ln)*zmax((-sum_r_s(ijk,m)))*vol(ijk) + &
                       r_sp(ijk,m,ln)*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 (x_s(1,m,ln), dif_s(1,m,ln), &
                    discretize(7), 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
                  CALL conv_dif_phi (x_s(1,m,ln), dif_s(1,m,ln), &
                    discretize(7), 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 (x_s(1,m,ln), bc_x_s(1,m,ln), &
                  bc_xw_s(1,m,ln), bc_hw_x_s(1,m,ln), &
                  bc_c_x_s(1,m,ln), m, a_m, b_m)

               CALL source_phi (s_p, s_c, eps, x_s(1,m,ln), m, a_m, b_m)

! Add point sources.
               IF(point_source) CALL point_source_phi (x_s(1,m,ln), &
                  ps_x_s(:,m,ln), ps_massflow_s(:,m), m, a_m, b_m)

! usr sources
            IF(call_usr_source(8)) CALL calc_usr_source(solids_species,&
                                 a_m, b_m, lm=m, ln=ln)

               CALL calc_resid_s (x_s(1,m,ln), a_m, b_m, m, &
                  num_resid(resid_x+(ln-1),m), &
                  den_resid(resid_x+(ln-1),m), resid(resid_x+(ln-1),&
                  m), max_resid(resid_x+(ln-1),m), ijk_resid(resid_x+(ln-1),m), &
                  zero_x_gs)

               CALL under_relax_s (x_s(1,m,ln), a_m, b_m, m, ur_fac(7))

!               call check_ab_m(a_m, b_m, m, .false., ier)
!               write(*,*) resid(resid_x+(LN-1), m), &
!                  max_resid(resid_x+(LN-1), m), &
!                  ijk_resid(resid_x+(LN-1), m)
!               call write_ab_m(a_m, b_m, ijkmax2, m, ier)
!
!               call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1,-3,M),&
!                  1, DO_K, ier)

               CALL adjust_leq (resid(resid_x+(ln-1),m), leq_it(7), &
                  leq_method(7), leqi, leqm)

               CALL solve_lin_eq ('X_s', 7, x_s(1,m,ln), a_m, b_m, m,&
                  leqi, leqm, leq_sweep(7), leq_tol(7), leq_pc(7), ier)

! Check for linear solver divergence.
               IF(ier == -2) err_l(mype) = 131

               CALL bound_x (x_s(1,m,ln), ijkmax2)
!               call out_array(X_s(1,m,LN), 'X_s')

            END DO

            if(chi_scheme) call unset_chi()
         ENDIF ! check for any species in phase m
      END DO ! for m = 1, mmax
! end solids phases species equations
! ----------------------------------------------------------------<<<

      call unlock_ambm

! If the linear solver diverged, species mass fractions may take on
! unphysical values. To prevent them from propogating through the domain
! or causing failure in other routines, force an exit from iterate and
! reduce the time step.
      CALL global_all_sum(err_l, err_g)
      ier = maxval(err_g)

      RETURN
      END SUBROUTINE solve_species_eq