FDS Source, Verification, Manual, Utilities: Allow 2-3 steps mechanisms with CVODE

Manuals/FDS_Technical_Reference_Guide/Combustion_Chapter.tex

@@ -342,16 +342,16 @@ \subsection{Finite-Rate Chemistry (Detailed Chemical Mechanism)}
 \label{detailed_chemistry_using_mechanism}
 A chemical mechanism represents different chemical pathways using multiple species and reactions. Typically, these mechanisms are available in Chemkin or Cantera YAML formats. Such a mechanism can be represented as:
 \begin{equation}\label{eq:chemistry_mechanism}
-\sum_{j=1}^{N_{sp}}{\nu}_{ji}^{'} \ X_j \Leftrightarrow \sum_{j=1}^{N_{sp}}{\nu}_{ji}^{''} \ X_j, i=1,2,3,...,N_{reac}
+\sum_{j=1}^{N_{\text{s}}}{\nu}_{ji}^{'} \ X_j \Leftrightarrow \sum_{j=1}^{N_{\text{s}}}{\nu}_{ji}^{''} \ X_j, i=1,2,3,...,N_\text{r}
 \end{equation}
-Here, $X_j$ represents the chemical symbol of the $j$th species (i.e. $\mathrm{CH_4, H_2, O_2}$); ${\nu}_{ji}^{'}$ and ${\nu}_{ji}^{''}$ are the stoichiometric coefficients of the $j$th species in the $i$th reaction; ${N_{sp}}$ and ${N_{reac}}$ are the total number of species and total number of reactions, respectively.
+Here, $X_j$ represents the chemical symbol of the $j$th species (i.e. $\mathrm{CH_4, H_2, O_2}$); ${\nu}_{ji}^{'}$ and ${\nu}_{ji}^{''}$ are the stoichiometric coefficients of the $j$th species in the $i$th reaction; ${N_{\text{s}}}$ and ${N_\text{r}}$ are the total number of species and total number of reactions, respectively.
 The rate of change of the molar concentration \si{(kmol/m^3/s)} of each species can be represented using the following system of ordinary differential equations (ODEs):
 \begin{equation}\label{eq:chemistry_ode_system}
-\frac{\d C_k}{\d t} =  \dot{\omega_k} = \sum_{i=1}^{N_{reac}} b_{i} \ \nu_{ki} \ r_i
+\frac{\d C_k}{\d t} =  \dot{\omega_k} = \sum_{i=1}^{N_\text{r}} b_{i} \ \nu_{ki} \ r_i
 \end{equation}
 The right-hand side of the above equation represents contributions from each reaction to the $j$th species. Here, $C_k$ is the molar concentration \si {(kmol/m^3)} of the $k$th species; $b_i$ is the reaction rate modification coefficient of the $i$th reaction due to third-body effects and pressure; $\nu_{ki} = {\nu}_{ki}^{''} - {\nu}_{ki}^{'}$; and $r_i$ is the reaction progress rate of the $i$th reaction.
 \begin{equation}\label{eq:reaction_progress_rate}
-r_i =  k_{f,i} \prod_{j=1}^{N_{sp}} (C_j)^{{\nu}_{ji}^{'}}  -  k_{r,i} \prod_{j=1}^{N_{sp}} (C_j)^{{\nu}_{ji}^{''}}
+r_i =  k_{f,i} \prod_{j=1}^{N_{\text{s}}} (C_j)^{{\nu}_{ji}^{'}}  -  k_{r,i} \prod_{j=1}^{N_{\text{s}}} (C_j)^{{\nu}_{ji}^{''}}
 \end{equation}
 Here, $k_{f,i}$ and $k_{r,i}$ are the forward and reverse reaction rate coefficients, respectively. The $k_{f,i}$ can be calculated using Eq.~(\ref{eq:rate_cons}), and the $k_{r,i}$ can be calculated by obtaining the equilibrium constant, as shown in Eq.~(\ref{eq:equilibrium_const}).
 
@@ -365,7 +365,7 @@ \subsection{Finite-Rate Chemistry (Detailed Chemical Mechanism)}
 \end{equation}
 For third-body reactions, presence of (or collission with) a third body (any other molecule) modify the reaction rate. If the third-body efficiency of species $j$ in reaction $i$ (denoted as $\alpha_{ij}$) is provided, then
 \be
-\eta_{\mathrm{tb},i} = \sum_{j=1}^{N_{sp}} \alpha_{ij}C_{j}
+\eta_{\mathrm{tb},i} = \sum_{j=1}^{N_{\text{s}}} \alpha_{ij}C_{j}
 \ee
 If the third-body efficiency is not provided, the default value of $\alpha_{ij}$ is 1.0.
 
@@ -398,7 +398,7 @@ \subsection{Finite-Rate Chemistry (Detailed Chemical Mechanism)}
 
 To solve the reactive system, we assume a constant pressure reactor. For that, the concentration ODEs in the form of Eq.~(\ref{eq:chemistry_ode_system}) need to be solved along with a ODE of temperature given by:
 \begin{equation}\label{eq:TemperatureDerivative}
-\frac{\d T}{\d t} = \dot{T} = -\frac{1}{\rho c_p} \sum_{j=1}^{N_{sp}}h_j W_j \dot{\omega_j} 
+\frac{\d T}{\d t} = \dot{T} = -\frac{1}{\rho c_p} \sum_{j=1}^{N_{\text{s}}}h_j W_j \dot{\omega_j} 
 \end{equation} 
 Here,  $\rho$  is the density $\mathrm{(kg/m^3)}$, $c_p$ is the specific heat of the mixture (J/kg/K), $h_k$ is the absolute enthalpy that includes enthalpy of formation (J/kg), $W_j$ is the molecular weight (kg/kmol) of species $j$, and $\dot{\omega}$ is the species production rate given by Eq.~(\ref{eq:chemistry_ode_system}).
 To solve the system of ODEs, we first convert the mass fraction of species to molar concentrations using $C_j=Y_j\rho/W_j$. Then, CVODE from Sundials \cite{cvodeDoc:2024} is used to solve the system of ODEs by supplying an analytical Jacobian. The details of analytical Jacobian formulation is provided in the Appendix \ref{chemistry_analytical_jacobian}.

Source/chem.f90

@@ -17,6 +17,7 @@ MODULE CVODE_INTERFACE
 
 LOGICAL :: DEBUG=.FALSE.
 REAL(EB), PARAMETER :: MIN_CHEM_TMP=200._EB
+REAL(EB) :: CUR_CFD_TIME
 
 PUBLIC CVODE_SERIAL
 
@@ -27,8 +28,8 @@ MODULE CVODE_INTERFACE
 !> \param TN_C is the current time
 !> \param SUNVEC_Y is the current array of molar concentrations, temperature and pressure.
 !> \param SUNVEC_F is the array of derivatives returned
-!> \param USER_DATA ??
-
+!> \param USER_DATA is the user data array. Not yet used in FDS.
+!> \details The right hand side function of the ode d[c]/dt = wdot (=f). Provides the Derivative function to CVODE.
 INTEGER(C_INT) FUNCTION RHSFN(TN_C, SUNVEC_Y, SUNVEC_F, USER_DATA) &
     RESULT(IERR) BIND(C,NAME='RHSFN')
 
@@ -65,7 +66,7 @@ END FUNCTION RHSFN
 
 
 
-!> \brief Calculate derivative (fvec) for a given cvec(n_tracked_species+2)
+!> \brief Calculate derivative (fvec) for a given concentration vector cvec(n_tracked_species+2)
 !> \param CVEC is the current array of molar concentrations, temperature and pressure.
 !> \param FVEC is the array of derivatives returned
 
@@ -77,7 +78,7 @@ SUBROUTINE DERIVATIVE(CVEC,FVEC)
 REAL(EB) :: R_F,MIN_SPEC(N_TRACKED_SPECIES), KG,  TMP, RHO, &
           K_0, K_INF, P_RI, FCENT, B_I, RRTMP, THIRD_BODY_ENHANCEMENT, PR
 INTEGER :: I,NS, ITMP
-REAL(EB) :: ZZ(N_TRACKED_SPECIES), CP, HS_I, DG 
+REAL(EB) :: ZZ(N_TRACKED_SPECIES), CP, HS_I, DG
 TYPE(REACTION_TYPE), POINTER :: RN
 
 TMP = MAX(CVEC(N_TRACKED_SPECIES+1), MIN_CHEM_TMP)
@@ -161,7 +162,7 @@ END SUBROUTINE DERIVATIVE
 !> \brief Calculate fall-off function 
 !> \param TMP is the current temperature.
 !> \param P_RI is the reduced pressure
-!> \param RN is the reaction
+!> \param RNI is the index of reaction
 
 REAL(EB) FUNCTION CALCFCENT(TMP, P_RI, RNI)
 REAL(EB), INTENT(IN) :: TMP, P_RI
@@ -194,10 +195,16 @@ END FUNCTION CALCFCENT
 
 
 
-!> \the jacobian of the ode right hand side function j = df/dy
-!> \param TN_C is the current time
+!> \brief The jacobian of the ode right hand side function j = df/dy
+!> \param TN_C is the current time provided by CVODE during callback, not the actual CFD time.
 !> \param SUNVEC_Y is the current array of molar concentrations, temperature and pressure.
 !> \param SUNVEC_F is the array of derivatives returned
+!> \param SUNMAT_J is the Jacobian array returned to CVODE
+!> \param USER_DATA is the user data array. Not yet used in FDS.
+!> \param TMP1 is not yet used in FDS.
+!> \param TMP2 is not yet used in FDS.
+!> \param TMP3 is not yet used in FDS.
+!> \details The jacobian of the ode right hand side function j = df/dy. Provides the Jacobian function to CVODE.
 
 INTEGER(C_INT) FUNCTION JACFN(TN_C, SUNVEC_Y, SUNVEC_F, SUNMAT_J, &
  USER_DATA, TMP1, TMP2, TMP3) &
@@ -255,7 +262,7 @@ INTEGER(C_INT) FUNCTION JACFN(TN_C, SUNVEC_Y, SUNVEC_F, SUNMAT_J, &
 END FUNCTION JACFN
 
 
-!> \brief Calculate the jacobian matrix (jmat[n_tracked_species+2,
+!> \brief Calculate the analytical jacobian jmat[n_tracked_species+2,n_tracked_species+2]
 !> \param CVEC is the current array of molar concentrations, temperature and pressure.
 !> \param FVEC is the array of derivatives passed as input
 !> \param JMAT is the jacobian matrix returned
@@ -352,7 +359,7 @@ SUBROUTINE JACOBIAN(CVEC,FVEC, JMAT)
    !Contribution of qi
    DO NS1 = 1, RN%N_SPEC
       DO NS=1,RN%N_SMIX_FR
-         DCVEC1 = R_F*RN%NU_NN(RN%NU_INDEX(NS))/CVEC(YP2ZZ(RN%N_S_INDEX(NS1)))
+         DCVEC1 = R_F*RN%NU_NN(RN%NU_INDEX(NS))*RN%N_S(NS1)/CVEC(YP2ZZ(RN%N_S_INDEX(NS1)))
          JMAT((YP2ZZ(RN%N_S_INDEX(NS1))),RN%NU_INDEX(NS)) = &
          JMAT((YP2ZZ(RN%N_S_INDEX(NS1))),RN%NU_INDEX(NS))+ DCVEC1 
       ENDDO   
@@ -429,6 +436,7 @@ END SUBROUTINE JACOBIAN
 
 
 !> \brief Print the component of jacobian matrix
+!> \param JMAT is the Jacobian matrix.
 SUBROUTINE PRINT_JMAT(JMAT)
 REAL(EB), INTENT(IN) :: JMAT(N_TRACKED_SPECIES+2, N_TRACKED_SPECIES+2)
 INTEGER :: NS, NS2
@@ -477,10 +485,13 @@ END SUBROUTINE PRINT_JMAT
 
 !> \brief Calculate DBIDC of reactions 
 !> \param TMP is the current temperature.
-!> \param RN is the reaction
-!> \param CVEC is the array of concentrations
-!> \param PR is the PRESSURE RATIO.
-!> \param F is the FALLOFF FUNCTION VALUE.
+!> \param RNI is the reaction index.
+!> \param K0 is the low pressure rate coeff.
+!> \param KINF is the high pressure rate coeff.
+!> \param PR is the pressure ratio.
+!> \param F is the falloff function value.
+!> \param DBIDC is the derivative of modification factor w.r.t concentration (out).
+!> \param DBIDT is the derivative of modification factor w.r.t temperature (out).
 
 SUBROUTINE CALC_FALLOFF_DBIDC_AND_DBIDT(TMP, RNI, K0, KINF, PR, F, DBIDC, DBIDT)
 REAL(EB), INTENT(IN) :: TMP, PR, K0, KINF, F
@@ -515,6 +526,11 @@ SUBROUTINE CALC_FALLOFF_DBIDC_AND_DBIDT(TMP, RNI, K0, KINF, PR, F, DBIDC, DBIDT)
 END SUBROUTINE CALC_FALLOFF_DBIDC_AND_DBIDT
 
 !> \brief Calculate derivative of TROE function w.r.t concentration  
+!> \param PR is the pressure ratio.
+!> \param F is the falloff function value.
+!> \param DPRDC is the derivative of TROE function w.r.t concentration.
+!> \param TMP is the current temperature.
+!> \param RNI is the reaction index
 REAL(EB) FUNCTION DDC_TROE(PR, F, DPRDC, TMP, RNI) 
 REAL(EB), INTENT(IN) :: PR, F, DPRDC, TMP
 INTEGER, INTENT(IN) :: RNI
@@ -547,6 +563,11 @@ END FUNCTION DDC_TROE
 
 
 !> \brief Calculate derivative of TROE function w.r.t temperature  
+!> \param PR is the pressure ratio.
+!> \param F is the falloff function value.
+!> \param DPRDT is the derivative of TROE function w.r.t temperature.
+!> \param TMP is the current temperature.
+!> \param RNI is the reaction index
 REAL(EB) FUNCTION DDTMP_TROE(PR, F, DPRDT, TMP, RNI) 
 REAL(EB), INTENT(IN) :: PR, F, DPRDT, TMP
 INTEGER, INTENT(IN) :: RNI
@@ -586,15 +607,15 @@ END FUNCTION DDTMP_TROE
 
 
 
-!> \cvode interface for ODE integrator 
-!> \Call sundials cvode in serial mode.
+!> \brief cvode interface for ODE integrator. Call sundials cvode in serial mode.
 !> \param ZZ species mass fraction array
 !> \param TMP_IN is the temperature
 !> \param PR_IN is the pressure
 !> \param TCUR is the start time in seconds
 !> \param TEND is the end time in seconds
 !> \param RTOL is the relative error for all the species (REAL_EB)
 !> \param ATOL is the absolute error tolerance array for the species (REAL_EB)
+!> \details This is the interface subroutine to the other modules.
 
 SUBROUTINE CVODE_SERIAL(CC,TMP_IN, PR_IN, TCUR,TEND, RTOL, ATOL)
 USE PHYSICAL_FUNCTIONS, ONLY : MOLAR_CONC_TO_MASS_FRAC, CALC_EQUIV_RATIO
@@ -740,12 +761,13 @@ SUBROUTINE CVODE_SERIAL(CC,TMP_IN, PR_IN, TCUR,TEND, RTOL, ATOL)
 
    IF (IERR_C .NE. CV_SUCCESS) THEN
       IF (IERR_C == CV_TOO_MUCH_WORK) THEN
-         WRITE(LU_ERR,'(A, 3E12.2, A)')" WARN: CVODE took all internal substeps. TSTART, TEND, DT=", TCUR, TEND, (TEND-TCUR), &
+         WRITE(LU_ERR,'(A, 2E18.8, A)')" WARN: CVODE took all internal substeps. CUR_CFD_TIME, DT=", CUR_CFD_TIME, (TEND-TCUR), &
                      ". If the warning persists, reduce the timestep."
       ELSE
-         WRITE(LU_ERR,'(A, I4, A, 3E12.2, A)')" WARN: CVODE didn't finish ODE solution with message code:", IERR_C, &
-                        " and TSTART, TEND, DT=", TCUR, TEND, (TEND-TCUR), ". If the warning persists, reduce the timestep."
+         WRITE(LU_ERR,'(A, I4, A, 2E18.8, A)')" WARN: CVODE didn't finish ODE solution with message code:", IERR_C, &
+                        " and CUR_CFD_TIME, DT=", CUR_CFD_TIME, (TEND-TCUR), ". If the warning persists, reduce the timestep."
       ENDIF   
+
       IF (DEBUG) THEN
          CALL MOLAR_CONC_TO_MASS_FRAC(CC(1:N_TRACKED_SPECIES), ZZ(1:N_TRACKED_SPECIES))
          CALL CALC_EQUIV_RATIO(ZZ(1:N_TRACKED_SPECIES), EQUIV)
@@ -772,7 +794,12 @@ SUBROUTINE CVODE_SERIAL(CC,TMP_IN, PR_IN, TCUR,TEND, RTOL, ATOL)
 END SUBROUTINE CVODE_SERIAL
 
 
-!> \brief error handler function, such that CVODE() doesn't output the error directly to stderr. 
+!> \brief CVODE error handler callback function, such that CVODE() doesn't output the error directly to stderr.
+!> \param ERR_CODE The error code send by CVODE
+!> \param MOD_NAME The module name where error occured send by CVODE
+!> \param FUNC_NAME The functio name where error occured send by CVODE
+!> \param MESSAGE The error message
+!> \param USER_DATA User data, not used in FDS.
 SUBROUTINE FDS_CVODE_ERR_HANDLER( ERR_CODE, MOD_NAME, FUNC_NAME, MESSAGE, USER_DATA) &
    BIND(C,NAME='FDS_CVODE_ERR_HANDLER')
 INTEGER(C_INT), VALUE :: ERR_CODE
@@ -837,6 +864,8 @@ END SUBROUTINE FDS_CVODE_ERR_HANDLER
 
 
 
+!> \brief print CVODE statistics from a CVODE memory object.
+!> \param CVODE_MEM CVODE memory object.
 SUBROUTINE CVODESTATS(CVODE_MEM)
 
 !======= INCLUSIONS ===========

Source/fire.f90

@@ -8,7 +8,7 @@ MODULE FIRE
 USE COMP_FUNCTIONS, ONLY: CURRENT_TIME
 USE SOOT_ROUTINES, ONLY: SOOT_SURFACE_OXIDATION
 #ifdef WITH_SUNDIALS
-USE CVODE_INTERFACE
+USE CVODE_INTERFACE, ONLY: CUR_CFD_TIME, CVODE_SERIAL
 #endif
 
 IMPLICIT NONE (TYPE,EXTERNAL)
@@ -821,7 +821,6 @@ SUBROUTINE COMBUSTION_MODEL(T,DT,ZZ_GET,Q_OUT,MIX_TIME_OUT,CHI_R_OUT,CHEM_SUBIT_
          ! May be used with N_FIXED_CHEMISTRY_SUBSTEPS, but default mode is to use error estimator and variable DT_SUB
 
          RICH_EX_LOOP: DO RICH_ITER = 1,RICH_ITER_MAX
-
             DT_SUB = MIN(DT_SUB_NEW,DT-DT_ITER)
             ! FDS Tech Guide (E.3), (E.4), (E.5)
             CALL FIRE_RK2(A1,ZZ_MIXED,ZZ_0,ZETA_1,ZETA_0,DT_SUB,1,TMP_IN,RHO_HAT,CELL_MASS,TAU_MIX,&
@@ -869,6 +868,9 @@ SUBROUTINE COMBUSTION_MODEL(T,DT,ZZ_GET,Q_OUT,MIX_TIME_OUT,CHI_R_OUT,CHEM_SUBIT_
          DO NS =1,N_TRACKED_SPECIES
             ATOL(NS) = DBLE(SPECIES_MIXTURE(NS)%ODE_ABS_ERROR)
          ENDDO
+#ifdef WITH_SUNDIALS
+         CUR_CFD_TIME = T ! Set current cfd time in cvode, for logging purpose.
+#endif
          CALL  CVODE(ZZ_MIXED,TMP_IN,PRES_IN, T1,T2, GLOBAL_ODE_REL_ERROR, ATOL)
          Q_REAC_SUB = 0._EB
    END SELECT INTEGRATOR_SELECT
@@ -995,7 +997,7 @@ SUBROUTINE CVODE(ZZ, TMP_IN, PRES_IN,  TCUR,TEND, GLOBAL_ODE_REL_ERROR, ATOL)
 ! Convert to concentration
 CC = 0._EB
 DO NS =1,N_TRACKED_SPECIES
-  CC(NS) = RHO_IN*ZZ(NS)/SPECIES_MIXTURE(NS)%MW  ! [RHO]= kg/m3, [MW] = gm/mol = kg/kmol, [CC] = kmol/m3
+  CC(NS) = RHO_IN*ZZ(NS)/SPECIES_MIXTURE(NS)%MW  ! [RHO]= kg/m3, [MW] = g/mol = kg/kmol, [CC] = kmol/m3
 ENDDO
 WHERE(CC<0._EB) CC=0._EB
 

Source/read.f90

@@ -4512,7 +4512,7 @@ SUBROUTINE READ_REAC
             CRITICAL_FLAME_TEMPERATURE,HEAT_OF_COMBUSTION,HOC_COMPLETE,E,C,H,N,O, &
             N_T,NU(MAX_SPECIES),N_S(MAX_SPECIES),THIRD_EFF(MAX_SPECIES),&
             FUEL_C_TO_CO_FRACTION,FUEL_H_TO_H2_FRACTION,FUEL_N_TO_HCN_FRACTION,AUTO_IGNITION_TEMPERATURE, A_LOW_PR, &
-            E_LOW_PR,N_T_LOW_PR, A_TROE, T1_TROE, T2_TROE, T3_TROE, SUMNU
+            E_LOW_PR,N_T_LOW_PR, A_TROE, T1_TROE, T2_TROE, T3_TROE
 REAL(EB) :: E_TMP=0._EB,S_TMP=0._EB,ATOM_COUNTS(118),MW_FUEL=0._EB,H_F=0._EB,PR_TMP=0._EB
 LOGICAL :: L_TMP,CHECK_ATOM_BALANCE,REVERSE,THIRD_BODY
 TYPE(REACTION_TYPE), POINTER, DIMENSION(:) :: REAC_TEMP
@@ -4826,23 +4826,7 @@ SUBROUTINE READ_REAC
       ALLOCATE(RN%SPEC_ID_NU_READ(RN%N_SMIX))
       RN%SPEC_ID_NU_READ = 'null'
       RN%SPEC_ID_NU_READ(1:RN%N_SMIX)=SPEC_ID_NU(1:RN%N_SMIX)
-
-      SUMNU = 0._EB
-      DO NS = 1, RN%N_SMIX
-         IF (NU(NS) .LT. 0._EB) THEN
-            SUMNU = SUMNU - NU(NS)
-         ENDIF
-      ENDDO
-
-      IF (RN%REACTYPE==THREE_BODY_ARRHENIUS_TYPE) THEN
-         RN%A_SI = RN%A_IN*(1000._EB)**(-SUMNU)  ![kmol/m3]^(-sum(nu)). Here 1000 is for conversion from mol/cm3 to kmol/m3
-      ELSE
-         RN%A_SI = RN%A_IN*(1000._EB)**(1-SUMNU) ![kmol/m3]^(1-sum(nu))
-      ENDIF
-
-      IF (RN%REACTYPE == FALLOFF_LINDEMANN_TYPE .OR. RN%REACTYPE == FALLOFF_TROE_TYPE) THEN
-         RN%A_LOW_PR = RN%A_LOW_PR*(1000._EB)**(-SUMNU)
-      ENDIF
+
    ELSE SIMPLE_IF
       RN%N_SMIX = 3
       RN%N_SPEC_READ = 0
@@ -5254,6 +5238,21 @@ SUBROUTINE PROC_REAC_1
       RN%RHO_EXPONENT = RN%RHO_EXPONENT + RN%N_S(NS)
    ENDDO
 
+   IF(.NOT. RN%REVERSE) THEN
+      IF (RN%REACTYPE==THREE_BODY_ARRHENIUS_TYPE) THEN
+         RN%A_SI = RN%A_IN*(1000._EB)**(-RN%RHO_EXPONENT)  ![kmol/m3]^(-sum(nu)). Here 1000 is for conversion from mol/cm3 to kmol/m3
+      ELSE
+         RN%A_SI = RN%A_IN*(1000._EB)**(1-RN%RHO_EXPONENT) ![kmol/m3]^(1-sum(nu))
+      ENDIF
+
+      IF (RN%REACTYPE == FALLOFF_LINDEMANN_TYPE .OR. RN%REACTYPE == FALLOFF_TROE_TYPE) THEN
+         RN%A_LOW_PR = RN%A_LOW_PR*(1000._EB)**(-RN%RHO_EXPONENT)
+      ENDIF
+   ELSE
+      RN%A_SI = REACTION(RN%REVERSE_INDEX)%A_SI
+      RN%A_LOW_PR = REACTION(RN%REVERSE_INDEX)%A_LOW_PR
+   ENDIF
+
    RN%RHO_EXPONENT = RN%RHO_EXPONENT - 1._EB ! subtracting 1 accounts for division by rho in Eq. (5.40)
    RN%A_PRIME = RN%A_PRIME * 1000._EB*SPECIES_MIXTURE(RN%FUEL_SMIX_INDEX)%MW ! conversion terms in Eq. (5.37)