Analytically reduced chemistry via quasisteady state (QSS) assumption in PelePhysics¶
The QSS assumption¶
Given a detailed chemical mechanism, some species can sometimes be assumed to be in a quasisteady state (QSS). Formally, if species \(A\) is assumed to be in a QSS then
where \([A]\) is the molar concentration of species \(A\). If the set of QSS species is judiciously chosen, macroscopic quantities of interest such as ignition delay or laminar flame speed are only mildly affected by the QSS assumption [DRG2005].
The advantage of the QSS assumption¶
Using the elementary reactions of the chemical mechanism, it is possible to write a set of algebraic equations that relate QSS species to nonQSS species. In turn, it is not necessary to transport the QSS species in the flow solver, as their concentration can be directly obtained from the nonQSS species.
In addition, QSS species are typically species that induce stiffness in the chemical mechanism since they evolve on scales different that the nonQSS species.
From nonQSS species to QSS species concentration¶
The set of algebraic equations that result from the QSS assumption does not always easily result in a simple algebraic relation between a given QSS species and all the nonQSS species. To do this, it is necessary to apply some form of linearization to the reactions involved [SYS2006]. In addition, even in the presence of linear relations, QSS species may depend on one another. The following figure shows the relation graph between the QSS species of a reduced \(NC_{12}H_{26}\) mechanism [ND2018]. The arrows indicate that a QSS species concentration depends on another QSS species concentration. It can also be seen that dependency groups exist among the QSS species. In PelePhysics it is possible to deduce invert analytically the linear system of equations given an arbitrary relation graph.
Linearizing the set of equations¶
In general, the QSS assumption results in a set of equations that are nonlinearly coupled making it difficult to invert the system. The nonlinear relations can arise if two or more QSSspecies are on the same side of an elementary reaction, or if the stoichiometric coefficient of a QSSspecies is not equal to one. Below, the relation graph between the QSS species plotted above is expanded with dots that denote reactions that relate QSSspecies. Dots or species are colored in red if they are involved in a quadratic coupling.
From here, it is necessary to eliminate the quadratic coupling to linearize the set of algebraic equations that result from the QSS assumption. Three methods can be envisioned: either one relaxes the QSS assumption by reducing the set of QSS species (Method 1) or one can eliminate the reactions that induce quadratic coupling. In case reactions are reversible, it can happen that either the forward or the backward reaction induces the quadratic coupling. Either one can remove both the forward and the backward reaction (Method 2) or remove either the backward or the forward reaction (Method 3).
All three methods are available in PelePhysics. By default, Method 3 is activated as it is computationally efficient and accurate (as will be shown below). Method 1 is the most accurate and Method 2 is the most computationally efficient in our experience. Given that each method has its advantage, we decided to allow the user to choose either one according to his/her needs.
Validation¶
The three linearization methods are validated against the skeletal \(NC_{12}H_{26}\) [SKEL2017]. Using 343 0D calculation that span the range of applicability of the QSS assumption (\(\phi = [0.5, 2.0], p=[1atm, 50atm], T=[800K, 1600K]\)), the ignition delay is computed using the skeletal mechanism (SK53) and each one of the three linearization methods for the QSS mechanism (RedXX). The left plot shows the correlation between the ignition delay from the skeletal mechanism and the reduced version. The statistics of the relative error between the reduced and the skeletal mechanism are shown in the title of that plot. The right plot shows the ignition delay values at highpressure conditions only.
Analytical Jacobian¶
In several computational experiments, using analytical Jacobians was found to provide better stability or efficiency compared with finite difference approximation or numerical inversion (see also fig:qss_integrator). Compared with nonQSS mechanisms, analytical Jacobians need to reflect the dependence of each QSS species on nonQSS species. However, QSS species may depend on an ensemble of other nonQSS species and therefore ensemble of reactions. Therefore, analytical Jacobian cannot be constructed by sequentially adding the contribution of each reaction. This significantly complicates the analytical jacobian generation. Failure to include the dependence of QSS species with respect to nonQSS species typically results in wrong ignition profiles, unless very small timesteps are used, as seen in figure fig:qss_aj.
To ease the implementation of analytical Jacobian in presence of QSS species, a symbolic approach is used to construct the analytical Jacobian. This strategy has the advantage of not requiring complex logic, and being flexible and readable for future development. For the last row of the Jacobian (partial difference of reaction rate with respect to temperature), finite difference is used since perturbations in temperature are less susceptible to numerical errors than perturbations in species concentrations. During the construction of the reaction rates, the operations printed to file are recorded symbolically using the sympy
and symengine
library [SYMPY]. For computational efficiency during the symbolic differentiation, the chainrules terms are computed and the final expressions are computed and assembled by chainruling using logic internal to CEPTR rather than sympy
. We have found that this speeds up the Jacobian construction cost by a factor 10.
To activate the use of symbolic jacobian, one needs to pass the flag qss_symbolic_jacobian
to CEPTR.
Printing the Jacobian terms one by one is not possible since the expressions that include QSS species are typically very large. Instead, the expressions are reduced via common subexpression precomputing that are later used in each term of the Jacobian. The number of subexpressions may be orders of magnitude larger than the number of Jacobian entries which can be problematic if the computational architecture has limited memory.
Several formatting strategies have been implemented to mitigate the memory footprint of the symbolic Jacobian. They can be adjusted by providing a .toml
file to CEPTR via the qss_format_input
flag. A model input file is provided in the Tutorials section of this documentation. A model execution script for generating a mechanism for dodecane_lu_qss
is available under the Tutorials section of this documentation.
The formatting options are the following
hformat
(string)cpu
will print intermediate variables used for chainruling. This gives a “readable” version of the Jacobian entries, albeit memory consuming.gpu
will not print intermediate variables used for chainruling, and instead will replace them directly in the Jacobian entries. This gives a less readable version of the Jacobian, but more memory efficient.
remove_1
(boolean) will replace expressions of the type1.0 * xxx
intoxxx
.remove_pow
(boolean) will convert expressions of the typepow(xxx,n)
into multiplications or division. The conversion occurs forn<=3
andn>=3
consistent with optimCudaremove_pow10
(boolean) will convert expressions of the typepow(10,xxx)
intoexp(ln(10)*xxx)
, consistent with optimCudamin_op_count
(integer) counts number operations used to construct each common subexpression and replace the common subexpression if the number of operations is less or equal ton
min_op_count_all
(integer) is similar tomin_op_count
but also counts how many times that common subexpression is used later. The meaning ofn
is different than formin_op_count
as it refers to how many more operations will be done if the common subexpression is eliminated. This option should be prefered tomin_op_count
as it tends to only marginally increase the file size (therefore compile time), while still being memory efficient.gradual_op_count
(boolean) is useful ifmin_op_count
ormin_op_count_all
are active. It loops from 1 ton
and gradually eliminate the common subexpressions. This has the advantage of ensuring that the memory footprint is strictly monotonically decreasing as n is increased.store_in_jacobian
(boolean) will use the Jacobian array as a temporary space to store intermediate variables. In particular, the last row of the Jacobian (dependence with respect to temperature) is done by finite difference which requires storing intermediate variables (production rate, forward and backward reactions). When the option is active, theproductionRate
function used to compute the finite difference is replaced with aproductionRate_light
functions where references to different parts of the Jacobian are used in place of allocating new arrays.round_decimals
(boolean) will round floats printed bysympy
when possible to minimize character count in themechanism.H
file.recycle_cse
(boolean) will reuse subexpressions that are not used later to avoid declaring new temporary reals.remove_single_symbols_cse
(boolean) will remove common subexpressions that are made of 1 operation and 1 symbol. Those common subexpressions are typicallyxxx
and may not appear as worth replacing because they save 1 operations and are reused multiple times. However, when replaced in the later expressions, the
operations typically disappear or is merged into another operations which actually does not increase the total number of operations.
The analytical Jacobian for QSS mechanisms is typically more accurate and stable than GMRES, and is on par with the finite difference Jacobian of CVODE as seen in fig:qss_integrator
In terms of speed, the analytical Jacobian 0D reactor is faster on CPU than finite difference Jacobian and GMRES. For the piston bowl challenge problem, the analytical Jacobian relative speed depends on the prevalence of chemical reactions. At some points, the AJ is slower than GMRES with PeleC, at others, AJ is faster. In PeleLM cases, AJ was found to be faster than GMRES. Further optimization and tests are still ongoing.
Using the Analytical Jacobian¶
To use the analytic Jacobian QSS mechanisms in a Pele run, one must specify the correct solve_type
. The choice of the solve_type
will determine whether or not the analytic Jacobian is used. With CUDA, both magma_direct
and sparse_direct
will use the analytic Jacobian, while on HIP only magma_direct
will trigger the use of the analytic Jacobian parts.
[DRG2005] 

[SYS2006] 

[ND2018]  (1, 2, 3)

[SKEL2017] 

[SYMPY] 
