Optimization of one-parameter family of integration formulae for solving stiff chemical-kinetic ODEs

A fast and robust Jacobian-free time-integration method—called Minimum-error Adaptation of a Chemical-Kinetic ODE Solver (MACKS)—for solving stiff ODEs pertaining to chemical-kinetics is proposed herein. The MACKS formulation is based on optimization of the one-parameter family of integration formulae coupled with a dual time-stepping method to facilitate error minimization. The proposed method demonstrates higher accuracy compared to the method—Extended Robustness-enhanced numerical algorithm (ERENA)—previously proposed by the authors. Additionally, when this method is employed in homogeneous-ignition simulations, it facilitates realization of faster performance compared to CVODE.

www.nature.com/scientificreports/ The proposed study focuses on development of a Jacobian-free implicit method to leverage the advantage associated with reduction in computational time involved in Jacobian generation. This is all-the-more true with regard to cases involving sparse Jacobians. The primary objective of this research concerned the development of a fast and accurate integration method for solving stiff chemical-kinetic ODEs. Formulation of the proposed method has first been explained in this paper. The said formulation was based on optimization of the oneparameter family of integration equations (to facilitate error minimization) coupled with a dual time-stepping method. Lastly, the proposed integration method was incorporated into Zero-RK software 4 and the performance was investigated via comparing results obtained upon its use during ignition simulation against those obtained using the ERENA (proposed previously by the authors) 14 , conventional CVODE 10 , and the efficiency method developed by the LLNL group (hereafter referred to as LLNL) 16 .

Methodology
Chemical-kinetic ODEs and model equation. Chemical-kinetic ODEs under constant-volume adiabatic conditions can be expressed as where y denotes the mass fraction; subscript i denotes the ith chemical species; ω denotes the production rate; ρ denotes density; τ denotes the chemical characteristic time; c denotes the creation rate; and N s is the total number of species. A model equation was considered in this study to simplify the discussion concerning proposition of a new integration method for solving chemical-kinetic ODEs. For construction of the said model equation, the quasi-steady-state (QSS) assumption 11,14 was applied to Eq. (1); i.e., parameters c i and τ i were assumed to have constant values in Eq. (1). Equation (1) in terms of variable u with constants α and β can, thence, be expressed as The following passages describe the development of the proposed integration method based on the modified model equation Eq. (2).

Dual time-stepping Method (DTS).
To investigate characteristics of the above model equation, the exact solution to Eq. (2) can be expressed as wherein u 0 represents the initial value of u. When t → ∞ , u → αβ ; i.e., the product αβ is representative of a converged solution for u.
Next, Eq. (2) can be discretized as follows using the explicit Euler (EE) method.
Here, n denotes the current step number. From the viewpoint of stability analysis and positive mass fraction, h ≤ β . When h = β , Eq. (4) takes the form which is identical to the exact converged solution. This implies that a converged solution can be obtained by means of a single-step calculation when employing the EE technique. Nonlinear chemical-kinetic ODEs are expected to demonstrate the best convergence properties when h ≈ τ . It is to be noted that when adopting the EE approach for solving chemical-kinetic ODEs, the actual time-step size depends on the minimum characteristic time, thereby implying the time-step size to be invariably small. To overcome this limitation with regard to time-step size, the DTS was employed in this study. Application of the DTS method facilitates easy attainment of the converged solution u n+1 . From Eq. (2), the parameter F(u) can be expressed as Using the above result, the converged solution for F(u) = 0 can be expressed as u n+1 . To obtain F(u) = 0 , a steady-state solution to the following ODE must be obtained.
where t ′ denotes the pseudo time. Any value of the pseudo-time-step size h ′ can be selected to facilitate attainment of rapid convergence, whereas that of h influences time accuracy of the integration method. Therefore, different values of h ′ can be chosen corresponding to different chemical species. It can, therefore, be realized u n+1 = u n + β α − u n /β = αβ, www.nature.com/scientificreports/ that use of the DTS approach eliminates the dependency on the minimum characteristic time. Rewriting the expression for F(u) as the corresponding solution could be obtained using h ′ = β ′ in the single-step EE calculation. Consequently, the following subsection focuses exclusively on the development of an efficient solver for F(u).
One-parameter family of integration formulae. Characteristic of one-parameter family of integration formulae for stiff ODEs 24 . Equation (6) discretized using the one-parameter family of integration formulae can be expressed as where θ values of 0, 0.5, and 1 imply use of the implicit Euler (IE), trapezoidal rule (TR), and EE methods, respectively. The one-parameter family of integration formulae is considered to be A-stable if θ ≤ 1/2 , and if θ = 1/2 (i.e., when employing the TR-based method), Eq. (9) becomes one with 2nd order, thereby demonstrating incurrence of the least error. The TR-rule based ( θ = 0.5 ) approach is, therefore, typically employed in this regard. Using eigenvalues of Jacobian matrix, the simple stiff system can be expressed as With reference to Eq. (9), u n+1 can be expressed as where q = h . The exact solution to u = − u satisfies the relation Using Eqs. (11) and (12), the relation can be obtained. Equation (13) holds for cases wherein q ≪ 1 . However, when q ≫ 1 (i.e., when ODEs are very stiff), Eq. (13) is no longer applicable because The TR-based approach is, therefore, inaccurate when handling stiff ODEs ( q ≫ 1 ). To minimize the error incurred when handling stiff ODEs, the value of θ must be optimized. In accordance with Liniger's definition, the error induced based on the value of θ can be expressed as Based on the above equation, Liniger demonstrated that a value of θ = 0.122 is optimum with regard to minimizing the total error incurred. Figure 1 graphically demonstrates the relation between parameters q and E(θ) predicted via application of the IE, TR, and Liniger's (optimized one-parameter integration) approaches. Reference to the said figure demonstrates that the total error incurred when using the TR-based approach increases with increase in ODE stiffness, and that the EI-based method demonstrates comparatively higher accuracy when solving stiff ODEs. Additionally, as can be observed in the figure, the maximum error incurred when employing Liniger's optimized one-parameter integration approach is the smallest. It must, however, be noted that ODEs employed in numerical modeling of combustion simulations, in general, demonstrate stiffness. It can, therefore, be inferred that when attempting to solve stiff chemical-kinetic ODEs, the IE-based approach must be employed to facilitate error reduction. In contrast, the TR-based approach must be exclusively employed under non-stiff conditions.

One-parameter family of integration formulae to solve stiffness chemical-kinetic ODEs. When
Eq. (5) is applied to Eq. (15), the term which is not related to the exponential could not be deleted. Consequently, Eq. (5) was differentiated in this study. The resulting equation can be expressed as In addition, the new error estimation has to be applied to discuss the optimization because Eq. (16) can not be applied to Liniger's error definition. Therefore, we defined the new error estimation as  www.nature.com/scientificreports/ Finally, using Eqs. (16) and (17) along with the assumptions of t n = 0 and t n+1 = h , θ could be expressed using |E| = 0 as The trend concerning variation in the value of θ corresponding to changes in γ is depicted in Fig. 2. As can be seen in the figure, the value of θ can be optimized in accordance with the stiffness of the model equation. This implies that in the proposed approach, the TR-based method is automatically applied when ≪ 1 , and the same is true with regard to the IE method when ≫ 1 . Consequently, the proposed study's objective of developing a numerical method based on the said model equation can be considered accomplished.

Minimum-error adaptation of chemical-kinetic ODEs solver (MACKS). The proposed method
based on the optimized value of θ and coupled with the DTS approach has been named as "Minimum-error adaptation of chemical-kinetic ODEs solver (MACKS)". The corresponding model equations can be expressed as  www.nature.com/scientificreports/ The proposed MACKS approach employs the absolute-and relative-error thresholds (ATOL and RTOL, respectively) to facilitate error evaluation in a manner similar to the CVODE approach. To preserve estimation accuracy, the actual time-step size h was divided by 4 for cases wherein the evaluated error increased within the inner loop.

Results and discussions
Numerical conditions. This study compares the proposed MACKS approach against previously proposed ERENA, CVODE, and LLNL methods in terms of their prediction accuracy and computational cost with regard to solving several zero-dimensional, homogeneous ignition problems involving stoichiometric mixtures of H 2 /air and CH 4 /air -n-C 12 H 26 /air maintained under conditions of p 0 = 0.1 and 1.0 MPa and T 0 = 1300 K . Values of RTOL and ATOL corresponding to the MACKS, CVODE, and LLNL approaches were 1 × 10 −5 and 1 × 10 −13 , respectively. The corresponding threshold value for ERENA was set as 1 × 10 −8 . The tolerances and threshold are the same as in our previous study 14 . Detailed reaction mechanisms considered in this study were generated in accordance with the Knowledge-basing Utilities for Complex Reaction Systems (KUCRS) 25 , and details concerning the number of chemical species and reactions considered have been summarized in Table 1.
It must be noted that the Jacobian matrix concerning the CVODE approach was intentionally initialized during each time step owing to the requirement of an initialization process when performing multidimensional CFD simulations. The corresponding base time-step sizes were h = 1 × 10 −8 and 1 × 10 −7 s. Accuracy and computational cost. In this study, the accuracy with regard to the ignition-delay time estimated when employing the MACKS approach was first compared against that estimated using the combined  where IDT represents the ignition-delay time and subscript TIM corresponds to LLNL, ERENA, or MACKS. As can be realized, if the value of "Error" is small, results obtained using the proposed approach demonstrate good agreement with those obtained using CVODE. Reference to Fig. 3 demonstrates the LLNL aproach was as accurate as the CVODE aproach, and the MACKS approach to be more accurate compared to ERENA. In particular, when the number of chemical species is small, use of the MACKS approach drastically reduces the error incurred when compared against the ERENA approach. This is because when employing the ERENA approach, the mass fraction is divided at every time step to conserve the summation of mass fractions 14 in accordance with the equation Using the ERENA approach's error threshold ε , the resulting summation of mass fractions can be expressed as www.nature.com/scientificreports/ The optimization error associated with the use of the ERENA approach can, therefore, be expressed as Consequently, the maximum error incurred per mass fraction can be represented in terms of ε/N s . Thus, the ERENA approach can be considered less accurate when employed in cases involving a small number of species. Next, the computational cost associated with the MACKS approach was compared against that incurred when employing the CVODE, LLNL, and ERENA methods with h = 1 × 10 −8 and 1 × 10 −7 s . The ratio of the total computational time could be defined as where t CPU denotes the total computational time.
Reference to the above equation demonstrates that when "Efficient ratio" equals less than 1, the corresponding computational cost exceeds that incurred when using the CVODE method; likewise, when "Efficient ratio" exceeds unity, the associated computational cost is lower compared to that incurred using the CVODE approach. Corresponding results are depicted in Fig. 4. As can be observed, the LLNL aproach is drastically faster compared to CVODE, but the MACKS approach is even faster compread to LLNL. However, the MACKS approach is slower compared to ERENA for cases wherein the number of chemical species involved is large. Accordingly, if only a small number of chemical species are involved, performance of the MACKS approach is comparable to that of ERENA. www.nature.com/scientificreports/ To facilitate better understanding of characteristics of the MACKS approach, the following subsection discusses the case involving a CH 4 /air mixture maintained at a pressure of 1 MPa. This is because chemical-kinetic models, typically, are reduced to involve less than 100 species, since the computational cost incurred becomes impractical when considering more than 100 species. Additionally, a pressure of 1 MPa is considered to be appropriate with regard to actual applications.
Accuracy and computational cost of CH 4 /air case. Time histories of temperature variations concerning the CH 4 /air mixture are depicted in Fig. 5 for cases involving h values of = 1 × 10 −8 and 1 × 10 −7 s . As can be observed, results obtained using the MACKS approach demonstrate good agreement with those obtained using CVODE and LLNL. However, results obtained using ERENA demonstrate large deviations from those obtained using CVODE, especially when h = 1 × 10 −7 s. Figure 6 depicts time histories of the minimum mass fraction concerning the the CH 4 /air mixture for cases with h = 1 × 10 −8 and h = 1 × 10 −7 s . As can be observed, use of the ERENA and MACKS approaches demonstrates zero or positive values of the minimum mass fraction, whereas the CVODE and LLNL method demonstrates negative values of corresponding parameters. (Note that the negative values obtained by CVODE are on the order of 10 −19 , which is very small compared to the LLNL results and therefore difficult to read from Fig. 5). This is because both ERENA and MACKS approaches are based on exact solutions of model equations. Additionally, the CVODE approach could produce negative mass fractions 27 .
Lastly, time histories concerning the computational cost incurred per iteration when employing the MACKS, ERENA, CVODE, and LLNL approaches are depicted in Fig. 7 for cases with h = 1 × 10 −8 and h = 1 × 10 −7 s . As can be observed, the peak computational cost corresponding to the MACKS approaches is much smaller compared to those corresponding to the CVODE, LLNL, and ERENA approaches with h = 1 × 10 −8 s . The MPI library is widely employed in high-performance computing, wherein load balancing forms one of the most important considerations. Figure 7 reveals that use of the MACKS approach demonstrates much better load-balancing characteristics compared to the ERENA method. This is because the difference in computational costs between flame zone and the other is small when the MACKS approach is employed during CFD analysis. In addition, the peak computational cost associated with the use of the MACKS approach is nearly identical to that associated with CVODE when h = 1 × 10 −7 s . In other words, the MACKS approach can be expected to demonstrate higher performance compared to the CVODE method even with less than 100 chemical species and a large time interval, and to have much higher performance under this condition than ERENA.

Conclusions
This study proposes a fast and robust Jacobian-free time-integration method called MACKS to facilitate treatment of stiff chemical-kinetic ODEs. The proposed approach was derived using the optimized one-parameter family of integration formulae to facilitate minimization of aggregate error. Results obtained in this study demonstrate that use of the MACKS approach leads to lower computational cost associated with solving zero-dimensional ignition problems whilst achieving higher prediction accuracy compared to the CVODE approach. In addition, the MACKS has been demonstrated to possess higher accuracy compared to the ERENA method previously proposed by the authors.

Appendix A: Accuracy and computational cost of n-C 7 H 16 /air case under engine-like conditions
MACKS is comared against CVODE, LLNL, and ERENA methods in terms of their prediction accuracy and computational cost with regard to solving several zero-dimensional, homogeneous ignition problems involving stoichiometric mixtures of n-C 7 H 16 /air maintained under conditions of p 0 = 2.0 MPa and T 0 = 800 and 1300 K. The corresponding base time-step sizes were h = 1 × 10 −8 s.
The "Efficient ratios" are summarized in Table 2. As can be observed, MACKS is about twice as slow as ERENA, but MACKS is about 8-18 times faster than LLNL.
Time histories of temperature variations concerning the n-C 7 H 16 /air mixture are depicted in Fig. 8. As can be observed, results obtained using the MACKS approach demonstrate good agreement with those obtained using CVODE, LLNL, and ERENA.
Time histories concerning the computational cost incurred per iteration when employing the CVODE, LLNL, ERENA, and MACKS approaches are depicted in Fig. 9. As can be observed, the peak computational cost corresponding to the MACKS approaches is smaller compared to those corresponding to the CVODE, LLNL, and ERENA. www.nature.com/scientificreports/