stiff_chem_mod.f

Go to the documentation of this file.

      MODULE stiff_chem

! External Routines.
!---------------------------------------------------------------------//
! Routine used to calculate the reaction rates and populate the
! fluid variable ODEs for ODEPACK.
      external stiff_chem_rrates
! Routine used to compare to values.
      LOGICAL, external :: compare
! Routine used to calculate species enthalpies.
      DOUBLE PRECISION, external :: calc_h0


! Runtime Flags:
!---------------------------------------------------------------------//
! Flag to invoke stiff chemistry solver.
      LOGICAL :: stiff_chemistry
! Flag to invoke variable solids density.
!      LOGICAL :: VARIABLE_DENSITY
! Flag indicating if cell IJK is own by myPE.
      LOGICAL, dimension(:), allocatable :: notowner

! ODEPACK Controlling parameters:
!---------------------------------------------------------------------//
! Dimension of ODEs solved in stiff solver.
      INTEGER :: ode_dimn_all
! Dimension of ODEs solved in stiff solver for gas phase only.
      INTEGER :: ode_dimn_g

! Dimension of ODEs solved in stiff solver for gas phase only.
      INTEGER :: neq_dimn

! Indicates type of Error control.
      INTEGER :: ode_itol
! Relative error tolerance paramter.
      DOUBLE PRECISION, DIMENSION(1) :: ode_rtol
! Absolute error tolerance parameter. (Dimension (ODE_DIMN))
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: ode_atol
! Declared length of RWORK.
      INTEGER :: ode_lrw
! Declared length of IWORK.
      INTEGER :: ode_liw
! Jacobian type indicator.
      INTEGER :: ode_jt
! The maximum number of steps ODEPACK may use to integrate.
      INTEGER :: stiff_chem_max_steps
! Flag indicating that the max number of steps is unlimited.
      LOGICAL :: unlimited_steps

! Explicit interface for ODEPACK
!---------------------------------------------------------------------//
      INTERFACE
         SUBROUTINE dlsoda (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, &
            itask,istate, iopt, rwork, lrw, iwork, liw, jac, jt)
            external f
            INTEGER :: ITOL, ITASK, ISTATE, IOPT, LRW, LIW, JT
            INTEGER, dimension(2) :: NEQ
            INTEGER, dimension(LIW) :: IWORK
            DOUBLE PRECISION :: T, TOUT
            DOUBLE PRECISION :: JAC
            DOUBLE PRECISION, dimension(1) :: RTOL
            DOUBLE PRECISION, dimension(LRW) :: RWORK
            DOUBLE PRECISION, dimension(NEQ(1)) :: Y, ATOL
         END SUBROUTINE dlsoda
      END INTERFACE

      contains

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                         C
!     Module name: MCHEM_TIME_MARCH                                       C
!     Purpose: Called in time_march.f to do rxns calcs                    C
!                                                                         C
!     Author: Nan Xie                                   Date: 02-Aug-04   C
!     Reviewer:                                         Date:             C
!                                                                         C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE stiff_chem_solver(ODE_DT, iErr)

! External Module Procedures:
!---------------------------------------------------------------------//
      use stiff_chem_maps, only : mapodetomfix
      use stiff_chem_maps, only : mapmfixtoode

! Global Variables:
!---------------------------------------------------------------------//
      use output,   only : full_log
      use param1,   only : zero
      use run,      only : time

      use mpi_utility

      use stiff_chem_dbg
      use stiff_chem_stats
      use functions

      implicit none

! Passed Variables:
!----------------------------------------------------------------------!
! Time integral length.
      DOUBLE PRECISION, intent(IN) :: ODE_DT
! Error Flag
      INTEGER, intent(OUT) :: iErr

! Local Variables:
!----------------------------------------------------------------------!
! Error flag -> Integration failed in one or more fluid cells.
      LOGICAL :: lErr_l   ! local

! Fluid Cell Index
      INTEGER :: IJK
! The maximum number of ODEs to solve.
      DOUBLE PRECISION, dimension(ODE_DIMN_all) :: ODE_VARS

! (1) :: Number of ODEs
! (2) :: Fluid cell index (IJK) passed into ODEPACK
! (:) :: Flag for solving solid phase M density and species equations.
      INTEGER, dimension(NEQ_DIMN) :: lNEQ
! Start time for integration
      DOUBLE PRECISION :: lT
! Stop time for integration
      DOUBLE PRECISION :: lTOUT
! Indicates type of Error control.
      INTEGER :: lITOL
! Relative error tolerance paramter.
      DOUBLE PRECISION :: lRTOL(1)
! Absolute error tolerance parameter. (Dimension (ODE_DIMN))
      DOUBLE PRECISION :: lATOL(ode_dimn_all)
! Index specifying the ODEPACK task.
      INTEGER :: lITASK
! Specifies the state of ODEPACK
      INTEGER :: lISTATE
! Flag indicating optional inputs are used.
      INTEGER :: lIOPT
! Array for REAL* work
      DOUBLE PRECISION :: RWORK(ode_lrw)
! Declared length of RWORK.
      INTEGER :: lLRW
! Array for Integer work
      INTEGER :: IWORK(ode_liw)
! Declared length of IWORK.
      INTEGER :: lLIW
! Jacobain Routine (not used)
      DOUBLE PRECISION :: lJAC
! Jacobian type indicator.
      INTEGER :: lJT

! The number of attempts of a specific fluid cell.
      INTEGER :: lAtps

      LOGICAL :: lReset
      LOGICAL :: lIncpt

      lerr_l = .false.

      CALL init_stiff_chem_stats

      ijk_lp: DO ijk = ijkstart3, ijkend3
         IF(notowner(ijk)) cycle ijk_lp
         IF(fluid_at(ijk)) THEN

            latps = 0
            lreset = .false.
            lincpt = .false.

! Forced restset of tolerance values.
            lrtol = ode_rtol
            latol = ode_atol

! Increment the attempts counter.
  50        latps = latps + 1

! Forced restset to initial values.
            lt    = 0.0d0
            ltout = ode_dt
            litol = ode_itol
            llrw  = ode_lrw
            lliw  = ode_liw
            ljt   = ode_jt

! Fixed parameters
            litask  = 1
            listate = 1
            liopt   = 1

! Calculate the number of ODEs to solve.
            CALL calc_ode_coeff(lneq, ijk)

! Clear the work arrays.
            iwork = 0
            rwork = zero

! The maximum number of internal steps ODEPACK may use to integrate over
! the time interval. The default value is 500.
            iwork(6) = stiff_chem_max_steps

            IF(calc_reactions(ijk)) THEN

! Map MFIX variables to ODE variables.
               CALL mapmfixtoode(neq_dimn, lneq, ode_dimn_all, ode_vars)

! Store a copy of the original field variables. This allows for these
! values to be 'reset' in the event that the stiff solver fails.
               CALL update_ode_old(ode_dimn_all, ode_vars)

! Clear the error flag.
 100           ierr = 0

! Integrate flow field variables to incorporate reactions.
               CALL dlsoda(stiff_chem_rrates, lneq, ode_vars, lt,      &
                  ltout, litol, lrtol, latol, litask, listate, liopt,  &
                  rwork, llrw, iwork, lliw, ljac, ljt)

! Verify that the results are well defined.
               CALL check_ode_data(neq_dimn, lneq, ode_dimn_all,       &
                  ode_vars, unlimited_steps, listate, ierr)

! Successfully Integrated ODEs.
               IF(ierr == 0) THEN
                  lreset = .false.
! Additional integration steps are needed (lT < lTOUT).
               ELSEIF(ierr == -1) THEN
! Reste the state flag and keep integrating.
                  IF(unlimited_steps) THEN
                     listate = 2
                     goto 100
                  ELSE
                     lreset = .false.
                     lincpt = .true.
                  ENDIF

! Too much accuracy was requested.
               ELSEIF(ierr == -2) THEN
                  IF(latps < 3) THEN
! Write to the error log file.
                     IF(ode_debug_level >= 2) CALL write_ode_log(ierr, &
                        neq_dimn, lneq, ode_dimn_all, ode_vars)
! Reset the ODE variable array.
                     CALL reset_ode(ode_dimn_all, ode_vars, latps)
! Reset the field variables.
                     CALL mapodetomfix(neq_dimn, lneq,               &
                        ode_dimn_all, ode_vars)
! Loosen the convergence criteria and try again.
                     lrtol = ode_rtol*10.0d0
                     latol = ode_atol*10.0d0
                     goto 50
                  ELSE
! Write to the error log file.
                     IF(ode_debug_level >= 1) CALL write_ode_log(ierr, &
                        neq_dimn, lneq, ode_dimn_all, ode_vars)
! Set the flag to reset the field variables to the initial values.
                     lreset = .true.
                  ENDIF
! All other errors.
               ELSE
! Tighten the convergence criteria and try again.
                  IF(latps < 3) THEN
! Write to the error log file.
                     IF(ode_debug_level >= 2) CALL write_ode_log(ierr, &
                        neq_dimn, lneq, ode_dimn_all, ode_vars)
! Reset the ODE variable array.
                     CALL reset_ode(ode_dimn_all, ode_vars, latps)
! Rest the filed variables to their original values.
                     CALL mapodetomfix(neq_dimn, lneq,                &
                        ode_dimn_all, ode_vars)
! Reduce the tolerances and try again.
                     lrtol = ode_rtol/(10.0d0**latps)
                     latol = ode_atol/(10.0d0**latps)
                     goto 50
                  ELSE
! Write to the error log file.
                     IF(ode_debug_level >= 1) CALL write_ode_log(ierr, &
                        neq_dimn, lneq, ode_dimn_all, ode_vars)
! Set the flag to reset the field variables to the initial values.
                     lreset = .true.
                  ENDIF

               ENDIF ! IF(iErr == 0)

! Reset the field variables.
               if(lreset) CALL reset_ode(ode_dimn_all, ode_vars, latps)
! Store the results in the field variables.
               CALL mapodetomfix(neq_dimn, lneq, ode_dimn_all, ode_vars)
! Collect solver stats.
               IF(full_log) CALL update_stiff_chem_stats(lneq, &
                  neq_dimn, iwork(11), ode_dimn_all, latps, lincpt)


            ENDIF  ! EndIF CALC_REACTIONS
         ENDIF  ! IF(CALC_REACTIONS(IJK))
      END DO ijk_lp ! End Loop over fluiod Cells, IJK

!      gErr_l = .FALSE.
!      CALL GLOBAL_ALL_OR(lErr_l, gErr_l)
!      IF(gErr_l) CALL WRITE_VTU_FILE


      CALL finalize_stiff_solver()

      IF(full_log) CALL write_stiff_chem_stats()

      ierr = 0

      RETURN
      END SUBROUTINE stiff_chem_solver


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: CALC_REACTIONS                                         !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!  Author: J.Musser                                   Date: 07-Feb-13  !
!                                                                      !
!  Comments:                                                           !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      LOGICAL FUNCTION calc_reactions(IJK)

      use rxns,   only : no_of_rxns

      implicit none

      INTEGER, intent(in) :: IJK

      DOUBLE PRECISION :: RATES(no_of_rxns)

      DOUBLE PRECISION, parameter :: rLimit = 1.0d-8

! Initialize
      rates = 0.0d0

! Calculate user defined reaction rates.
      CALL usr_rates(ijk, rates)

! If there is little to no reaction in the cell, then set the ODE
! Time to zero to avoid calling the stiff solver.
      calc_reactions = .true.
      if(maxval(rates) < rlimit) calc_reactions = .false.


      RETURN
      END FUNCTION calc_reactions


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: CALC_DIMN_ODE                                          !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!  Author: J.Musser                                   Date: 07-Feb-13  !
!                                                                      !
!  Comments:                                                           !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE calc_ode_coeff(lNEQ, IJK)

      use fldvar, only : ep_s
      use physprop, only : mmax
      use run, only : species_eq

      implicit none

      INTEGER, intent(in)  :: IJK
      INTEGER, intent(out) :: lNEQ(neq_dimn)

      INTEGER :: M

      LOGICAL :: USE_SOLIDS_ODEs

! Initialize.
      use_solids_odes = .false.
      lneq(2) = ijk
      lneq(3:) = 0

! If there is little to no solids in the cell, then set the ODE
! dimension to gas phase only.
      DO m=1, mmax
         IF(species_eq(m)) THEN
            IF(ep_s(ijk,m) > 1.0d-6) THEN
               use_solids_odes = .true.
               lneq(2+m) = 1
            ENDIF
         ENDIF
      ENDDO

      IF(use_solids_odes)THEN
         lneq(1) = ode_dimn_all
      ELSE
         lneq(1) = ode_dimn_g
         lneq(3:) = 0
      ENDIF


      RETURN
      END SUBROUTINE calc_ode_coeff


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: FINALIZE_STIFF_SOLVER                                  !
!                                                                      !
!  Purpose:                                                            !
!                                                                      !
!  Author: J.Musser                                   Date: 07-Feb-13  !
!                                                                      !
!  Comments:                                                           !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE finalize_stiff_solver

! Gas phase volume fraction
      use fldvar, only : ep_g
! Gas phase denisty
      use fldvar, only: ro_g, rop_g
! Gas phase pressure
      use fldvar, only: p_g
! Gas phase temperature
      use fldvar, only: t_g
! Gas phase species mass fraction
      use fldvar, only: x_g
! Molecular weight of each gas phase species
      use physprop, only : mw_g
! Gas phase mixture molecular weight.
      use physprop, only : mw_mix_g

! Solids phase bulk density
      use fldvar, only: rop_s
! Solids phase temperature
      use fldvar, only: t_s
! Solids phase species mass fractions
      use fldvar, only: x_s

! Number of solids phases
      use physprop, only: mmax
! Number of species in each phase.
      use physprop, only: nmax

! Double precision:  ONE = 1.0d0
      use param1, only: one
! Universal gas constant
      use constant, only : gas_const

      use compar
      use mpi_utility
      use sendrecv
      use functions
      use utilities

      implicit none

! Local loop indicies.
      INTEGER :: IJK  ! Fluid Cell index.
      INTEGER :: M    ! Solids phase index
      INTEGER :: NN    ! Species index

      CALL send_recv(ep_g,2)
      CALL send_recv(ro_g,2)
      CALL send_recv(rop_g,2)
      CALL send_recv(t_g,2)

      DO nn=1,nmax(0)
         CALL send_recv(x_g(:,nn),2)
         CALL bound_x (x_g(1,nn), ijkmax2)
      ENDDO

      DO m = 1, mmax
! Solids temperature.
         CALL send_recv(t_s(:,m),2)
! Solids volume fraction. (Constant Solids Density)
         CALL send_recv(rop_s(:,m),2)
! Solids phase species mass fractions.
         DO nn=1,nmax(m)
            CALL send_recv(x_s(:,m,nn),2)
            CALL bound_x (x_s(1,m,nn), ijkmax2)
         ENDDO
      ENDDO

      DO ijk = ijkstart3, ijkend3
         IF(.NOT.fluid_at(ijk)) cycle
! Calculate the mixture molecular weight.
         mw_mix_g(ijk) = sum(x_g(ijk,1:nmax(0))/mw_g(1:nmax(0)))
         mw_mix_g(ijk) = one/mw_mix_g(ijk)
! Calculate the gas phase pressure.
         p_g(ijk) = (ro_g(ijk)*gas_const*t_g(ijk))/mw_mix_g(ijk)
      ENDDO

      RETURN
      END SUBROUTINE finalize_stiff_solver

      END MODULE stiff_chem