check_odepack_stiff_chem.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Module name: CHECK_DATA_CHEM                                        !
!  Author: J.Musser                                   Date: 02-Aug-13  !
!                                                                      !
!  Purpose: Check the chemical rxns namelist variables.                !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE check_odepack_stiff_chem

! Global Variables:
!---------------------------------------------------------------------//
! Dimension of IJK arrays.
      use param, only: dimension_3
! Double precision zero.
      use param1, only: zero
! Double precision value for undefined variables.
      use param1, only: undefined
! Constant gas phase density
      use physprop, only : ro_g0
! Runtime logical for solving energy equations.
      use run, only: energy_eq
! Run time logical for solving species equations.
      use run, only: species_eq
! Net rate of gas phase production/consumption
      use rxns, only: sum_r_g
! Net rate of solids phase production/consumption
      use rxns, only: sum_r_s
! Run time logical for using stiff chemistry solver
      use stiff_chem, only: stiff_chemistry
! Run time logicals for identifying cells owned by myPE
      use stiff_chem, only: notowner

! Full access to the following modules:
!---------------------------------------------------------------------//
      use compar
      use geometry
      use indices

      use error_manager
      use functions

      implicit none

! Local Variables:
!---------------------------------------------------------------------//
      INTEGER :: I, J, K, IJK

      CALL init_err_msg('CHECK_ODEPACK_STIFF_CHEM')

! Verify that there is sufficient run complexity to use the stiff solver
      IF(stiff_chemistry) THEN

! Energy equations must be solved.
         IF(.NOT.energy_eq) THEN
            WRITE(err_msg,1004)'ENERGY_EQ = .FALSE.'
            CALL flush_err_msg(abort=.true.)
         ENDIF

! Must be compressible
         IF(ro_g0 /= undefined) THEN
            WRITE(err_msg,1004)'RO_G0 /= UNDEFINED'
            CALL flush_err_msg(abort=.true.)
         ENDIF

         IF(.NOT.species_eq(0)) THEN
            WRITE(err_msg,1004)'SPECIES_EQ(0) = .FALSE.'
            CALL flush_err_msg(abort=.true.)
         ENDIF

! The stiff chemistry solver only needs to loop over physical cells
! owned by a process (e.g., not ghost cells). To avoid having a
! triple do loop, this array is populated to identify the cells that
! are not owned.
         ALLOCATE( notowner(dimension_3) ); notowner = .true.
         do k=kstart, kend
         do j=jstart, jend
         do i=istart, iend
            ijk = funijk(i,j,k)
            notowner(ijk) = .false.
         enddo
         enddo
         enddo

! Initialize ODEPACK operating parameters.
         CALL odepack_init

! Clear the interphase mass transfer terms as the stiff solver
! does no use them.
         IF(allocated(sum_r_g)) sum_r_g = zero
         IF(allocated(sum_r_s)) sum_r_s = zero

      ENDIF

      CALL finl_err_msg



 1004 FORMAT('Error 1004: ',                                           &
         'Invalid parameters for stiff chemistry solver!',//           &
         ' The following criteria must be satisfied:',/                &
         '   > Solving the energy equations.',/                        &
         '   > Compressible gas phase.',/                              &
         '   > Solving gas phase species equations.',//                &
         ' >>> Invalid Parameter: ',a,//                               &
         'Check the user documentation and correct the mfix.dat file.')

      RETURN
      END SUBROUTINE check_odepack_stiff_chem



!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                         C
!     Module name: MCHEM_ODEPACK_INIT                                     C
!     Purpose: controlling values for ODEAPCK(reference to ODEPACK manual)C
!                                                                         C
!     Author: Nan Xie                                   Date: 02-Aug-04   C
!     Reviewer:                                         Date:             C
!                                                                         C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE odepack_init

! Global Variables:
!---------------------------------------------------------------------//
      use physprop, only: mmax
      use physprop, only: nmax

      use run, only: species_eq

      use stiff_chem, only: neq_dimn
      use stiff_chem, only: ode_dimn_g
      use stiff_chem, only: ode_dimn_all

      use stiff_chem, only: ode_itol
      use stiff_chem, only: ode_rtol
      use stiff_chem, only: ode_atol

      use stiff_chem, only: ode_lrw
      use stiff_chem, only: ode_liw
      use stiff_chem, only: ode_jt

      use stiff_chem, only: stiff_chem_max_steps

      use stiff_chem_dbg, only: allocate_stiff_chem_dbg
      use stiff_chem_stats, only: allocate_stiff_chem_stats

! Double precision value for undefined variables.
      use param1, only: undefined_i

      implicit none


! Local Variables:
!---------------------------------------------------------------------//
      INTEGER :: LRN, LRS
      INTEGER :: M

      LOGICAL :: LIMIT_MAX_STEPS

! Number of ODEs (maximum)
      neq_dimn = 2 + mmax

! Dimension of ODEs solved if only the gas phase is present:
! Gas density, temperature, and species
      ode_dimn_g = 2 + nmax(0)
! Solids temperature excluding for all phases.
      ode_dimn_g = ode_dimn_g + mmax

! Calculate the total number of ODEs that are solve.
      ode_dimn_all = ode_dimn_g
      DO m=1, mmax
! Solids bulk density and species.
         IF(species_eq(m)) ode_dimn_all = ode_dimn_all + (1 + nmax(m))
      ENDDO

! Indicates type of Error control.
      ode_itol = 2 ! :: EWT(i) = RTOL * ABS(Y(i)) * ATOL(i)

! Relative error tolerance parameter.
      ode_rtol(1) = 1.0d-5

! Absolute error tolerance parameter.
      IF(.NOT.(allocated(ode_atol))) allocate(ode_atol(ode_dimn_all))
      ode_atol(:) = 1.0d-6  ! All Equations

! Declared length of RWORK.
      lrn = 20 + 16*ode_dimn_all
      lrs = 22 + 9*ode_dimn_all + (ode_dimn_all**2)
      ode_lrw = max(lrn, lrs)

! Declared length of IWORK.
      ode_liw = 20 + ode_dimn_all

! Jacobian type indicator.
      ode_jt = 2 ! Internally generated.


      IF(stiff_chem_max_steps == undefined_i) THEN
         stiff_chem_max_steps = 500000
         limit_max_steps = .false.
      ELSE
         limit_max_steps = .true.
      ENDIF


      CALL allocate_stiff_chem_dbg(ode_dimn_all)
      CALL allocate_stiff_chem_stats

      return
      END SUBROUTINE odepack_init