Automated Time Stepping

A means of “automated time stepping” in RMG was implemented in order to eliminate the need for the user to specify times/conversions for mechanism validity testing. The “AUTO” option implementing this “automated time stepping” was successfully integrated into RMG starting in version 3.00. Several tests have been carried out with the new “AUTO” option. Results were promising, suggesting that the option functions as intended without being unreasonably computationally expensive. The “AUTO” option offers a more rigorous implementation of RMG’s rate-based mechanism generation algorithm while allowing straightforward generation of reaction mechanisms that encompass chemistry at all time-scales of interest.

Motivation

Originally, the implementation of RMG’s reaction mechanism generation algorithm required the user to specify times or conversions at which the program will estimate “edge” reaction fluxes and decide whether to include a new species in the “core” reaction model. The choice of which conversions or times to specify is, to a certain extent, arbitrary, and could have a significant impact on the model generated by the program. This arbitrariness can be effectively eliminated by keeping track of the “edge” reaction flux “continuously” within the ODE solver (analogously to how conversion is handled) and terminating ODE solver “integration” when the threshold reaction flux is reached (whereupon a species would be moved from the “edge” to the “core” as in the current RMG implementation). This new method is here referred to as “automated time stepping”. Implementation of such a method would reduce the number of adjustable parameters that the user must tweak during mechanism generation and provide a more rigorous use of the mechanism generation algorithm. At the same time, such a method would allow straightforward generation of mechanisms that encompass chemistry at all time-scales of interest.

Background

The reaction mechanism is divided into a “core” and an “edge”. The “core” is tracked by integrating the appropriate equations using an ordinary differential equation (ODE) solver. The “edge” includes species that are products of reactions among the “core” species. Periodically during the course of simulating species concentrations using the ODE solver, mechanism validity is tested. Mechanism validity criteria used by RMG are shown in the equations below.

\begin{eqnarray} \frac{{dC_j }}{{dt}} < tol \cdot R_{char} {\rm{ }}\forall {\rm{ }}j \in {\rm{\text{edge species}}} \\ {\rm{\text{with }}}R_{char} = \left[ {\sum\limits_i {\left( {\frac{{dC_i }}{{dt}}} \right)^2 } } \right]^{{\textstyle{1 \over 2}}} {\rm{ }}i \in {\rm{\text{core species}}} \end{eqnarray}

In these criteria, tol is a tolerance specified by the user. If any of these criteria are not valid, the edge species with the highest flux is added to the mechanism “core” and the ODE solver restarts simulation from t = 0. In the original implementation, the times or conversions at which model validity is tested are specified by the user. The ODE solvers used by RMG are FORTRAN programs (DASSL or DASPK), while RMG itself is written in Java.

Implementation

In order to achieve the desired behavior for “automated time stepping”, modifications of the type illustrated in Figure 1 were required.

Figure 1. Modifications to implement “automated time stepping”.

As the figure shows, a significant required change was the modification of the ODE solver to have access to edge species flux information at each ODE solver time step. Thus, much of the modification effort involved implementing a scheme for passing edge reaction information to the ODE solver. The implementation ended up requiring changes to the Java classes ReactionModelGenerator and JDASSL and the FORTRAN file call_dassl.f (the modified version was renamed to call_dasslAUTO.f90 to distinguish from the original version). Also, a minor change to the Java class ReactionTimeTT was made. An example of an ODE solver input file for “automated time stepping” is shown in Figure 2, below. Differences from “normal” operation are highlighted with boxes.

Figure 2. ODE solver input file example for “automated time stepping”. Boxes highlight differences from the “normal” version.

The first box in the above figure highlights an integer flag signalling to the FORTRAN code that “automated time stepping” is desired (as opposed to “normal” operation”). The next box highlights the section containing the information needed by the ODE solver to check mechanism validity at each step. In particular, this includes the user-specified tolerance, the number of edge species/reactions, and a line for each edge reaction containing the number of reactants, the number of products, ID numbers corresponding to reactants and products, and finally the rate coefficient, k, for the reaction. It should be noted that in the current implementation, the rate coefficient is considered to be a constant, as additional information regarding temperature/pressure dependence is not passed to solver. Thus, the “automated time stepping” code (along with other preexisting DASSL code) would need to be modified to handle non-isothermal cases (as well as non-isobaric pressure-dependent cases). In order to reproduce existing validity testing for cases with pressure-dependence, special considerations are necessary. In particular, each pressure-dependent network is considered as a “pseudo-species” with an associated “pseudo-reaction” using the kleak associated with the network as the rate coefficient. Monitoring flux to these “pseudo-species” thus allows us to consider whether a pressure-dependent network needs to be enlarged, analogously to how monitoring of species flux allows us to consider whether a species needs to be added to the “core” reaction mechanism.

Usage notes

The modified code still allows for evaluating flux at user-defined time/conversion points, and the use of “automated time stepping” is invoked by replacing the user-defined time/conversion points with the keyword “AUTO”. Thus, the use of “automated time stepping” will be referred to as the “AUTO” method in subsequent discussion. The modified code has been designed to work for both pressure-dependent and non-pressure-dependent cases; also, the code allows automatic stepping with either a conversion or a reaction time as the user-specified termination criterion.

Results

The new option of “AUTO” time-stepping was tested for several cases to verify that the option was working as desired and to investigate the effects on mechanism generation time and on the resulting mechanism. Note: Throughout this section, “1,3-hexadiene” refers to a system involving methane doped with 1,3-hexadiene (in particular, the system defined by the condition files distributed with RMG 2.00 (March 2007)). A brief summary of findings follows.

Effects on Generated Mechanism

One would expect that mechanisms generated using “AUTO” method would be more rigorously valid than those generated using arbitrary user-specified time/conversion steps. In the case of the 1,3-hexadiene system, mechanisms generated using the “AUTO” method were compared to mechanisms generated using evenly-spaced conversion steps. The results are shown in Table 1, below.

Table 1. Size of mechanisms generated using the “AUTO” method and varying numbers of evenly-spaced conversion steps. In all cases, the goal conversion is 0.9 and the tolerance is 0.05.

# of conversion steps	Core reactions	Core species
1	308	31
2	56	17
3	55	16
4	55	16
6	55	16
9	48	15
Auto	48	15

It should be noted that the mechanism generated with the “AUTO” method and the mechanism generated using nine evenly-spaced conversion steps appear to be essentially identical for this system. Thus, in this case, we approach “AUTO” method behavior as we specify a finer “grid” for conversion steps, as one might intuitively expect. It is also interesting to note that in the test case considered, the “AUTO” method produced a smaller mechanism than cases with a small number of specified conversion steps. Thus, in this case, the mechanism generated by the “AUTO” method (and cases with a sufficiently fine grid of conversion steps) would not be expected to be more accurate than the mechanism generated without any conversion steps (assuming the kinetics and thermochemistry of the extra species and reactions is accurate). However, the time required to generate the “AUTO” mechanism is significantly shorter, due to the smaller size of the mechanism. Furthermore, one can imagine that mechanisms of similar size to the largest mechanism reported above could be generated by using tighter tolerances and the “AUTO” method, and one could be more confident in their validity. In summary, although it produces a larger mechanism, the use of a low number of conversion steps is not “better” than the use of a high number of conversion steps or the “AUTO” method.

Timing Studies

Tests suggest that the time required to generate a mechanism with the “AUTO” method is virtually the same as the time required to generate the same mechanism with conventional use of conversion steps. However, the computational cost associated with tracking edge-species is significant compared to the cost of solving the ODEs, with the “AUTO” cases requiring a longer amount of time by a factor of two to eight (for the same number of ODE time steps). Three phenomena are proposed to explain this “discrepancy”. First, one would expect that in the “AUTO” case, the total time that we must integrate over during mechanism generation will be shorter, in general; that is, the “AUTO” method avoids “wasted” ODE time steps associated with “overshooting” the point where threshold flux has been reached. Second, the time associated with ODE solving is small in comparison to the total time requirement (< 10% of the total time for mechanism generation). Third, and perhaps most significant, there are fewer calls to the ODE solver when the “AUTO” method is used. In the case of the “AUTO” method, each call to the ODE solver is associated with a modification to the reaction mechanism (with the exception of the final integration). However, in the case of the “normal” method, there are extra calls to the ODE solver, as the integration is often stopped despite the mechanism being valid (i.e. no edge species fluxes exceed the threshold). Thus, in order to generate the same mechanism, the “normal” method requires a larger number of calls to the ODE solver, which, in turn is associated with greater overhead such as reading and writing ODE solver input and output files. Timing studies for the 1,3-hexadiene test case suggest that this overhead can be significant. (Note that the use of eighteen rather than nine conversions exacerbates the time cost, in accord with expectations based on the aforementioned overhead considerations.) The contribution of this overhead to FORTRAN time requirements is relatively small, but noticeable, while the contribution from the overhead on the Java side appears to be more significant. In fact, in this case, the total ODE solver time-requirement is noticeably lower in the “AUTO” case (relative to the “normal” cases). Thus, in this case, it appears that the added overhead associated with extra calls to the ODE solver in the “normal” method has outweighed the extra cost associated with tracking the flux at each time step within the ODE solver, using the “AUTO” method.

Summary

The “AUTO” method was implemented and tested during development of RMG 3.00 and provides the option of “automated time stepping”. Even so, timing studies suggest that tracking the edge flux “continuously” in the ODE solver is associated with a significant time penalty. However, various mitigating factors appear to reduce the impact of this time penalty on overall mechanism generation time when the “AUTO” method is used. Thus, the “AUTO” method appears to be an appealing alternative to user-specification of time/conversion steps.