Learning physics-based reduced-order models for a single-injector combustion process

TL;DR

This paper introduces a physics-informed machine learning approach to develop reduced-order models for a single-injector combustion process, achieving high accuracy and significant computational speedups without needing access to the high-fidelity model implementation.

Contribution

It combines physics-based model reduction with machine learning to identify transformed variables and learn quadratic ROMs directly from data, enhancing predictive capabilities.

Findings

ROM predicts key combustion variables accurately

Speedups of over five orders of magnitude compared to high-fidelity models

ROM remains predictive 200% beyond training interval

Abstract

This paper presents a physics-based data-driven method to learn predictive reduced-order models (ROMs) from high-fidelity simulations, and illustrates it in the challenging context of a single-injector combustion process. The method combines the perspectives of model reduction and machine learning. Model reduction brings in the physics of the problem, constraining the ROM predictions to lie on a subspace defined by the governing equations. This is achieved by defining the ROM in proper orthogonal decomposition (POD) coordinates, which embed the rich physics information contained in solution snapshots of a high-fidelity computational fluid dynamics (CFD) model. The machine learning perspective brings the flexibility to use transformed physical variables to define the POD basis. This is in contrast to traditional model reduction approaches that are constrained to use the physical variables of the high-fidelity code. Combining the two perspectives, the approach identifies a set of transformed physical variables that expose quadratic structure in the combustion governing equations and learns a quadratic ROM from transformed snapshot data. This learning does not require access to the high-fidelity model implementation. Numerical experiments show that the ROM accurately predicts temperature, pressure, velocity, species concentrations, and the limit-cycle amplitude, with speedups of more than five orders of magnitude over high-fidelity models. Our ROM simulation is shown to be predictive 200% past the training interval. Moreover, ROM-predicted pressure traces accurately match the phase of the pressure signal and yield good approximations of the limit-cycle amplitude.

Tables

$𝐀$  =	System matrix for linear part
$𝐁$  =	Input matrix
$𝐇$  =	Matricized quadratic tensor
$𝐪\left(t\right)$  =	State vector in finite dimensions
$𝐮\left(t\right)$  =	External input vector
$𝐐$  =	Snapshot matrix
$𝐕$  =	Matrix of POD basis vectors
$x,y$  =	Spatial coordinates
$d$  =	Number of physical variables
$n_x$  =	Spatial discretization dimension
$r$  =	Reduced model dimension
$t$  =	Time
${\overrightarrow{q}}_p(t,x,y),{\overrightarrow{q}}_c(t,x,y),{\overrightarrow{q}}_L(t,x,y)$  =	State vector in primitive, conservative, and learning variables
$p(t,x,y),𝐩\left(t\right)$  =	Pressure, continuous and discretized
$T(t,x,y),𝐓\left(t\right)$  =	Temperature, continuous and discretized
$\rho (t,x,y),𝝆\left(t\right)$  =	Density, continuous and discretized
$\xi (t,x,y),𝝃\left(t\right)$  =	Specific volume, continuous and discretized
$v_x(t,x,y),{𝐯}_x\left(t\right)$  =	Velocity in $x$ direction, continuous and discretized
$Y_l(t,x,y)$  =	Species mass fraction, $l=1,2,\dots ,n_{\text{sp}}$
$c_l(t,x,y)$  =	Species molar concentrations, $l=1,2,\dots ,n_{\text{sp}}$, also denoted as $\left[S\right]$
$\otimes$  =	Kronecker product
$\hat{\cdot }$  =	Notation for ROM quantities
Abbreviations
CFD  =	Computational Fluid Dynamics
GEMS  =	General equation and mesh solver; a CFD code
PDE  =	Partial differential equation
POD  =	Proper orthogonal decomposition
ROM  =	Reduced-order model

Table 2. Table 2 : Range of variable values for GEMS data.

State variable	Minimum	Mean	Maximum
Pressure $p$ in Pa	$9.226\times {10}^5$	$1.142\times {10}^6$	$1.433\times {10}^6$
Velocity $v_x$ in $\frac{\text{m}}{\text{s}}$	-222.930	69.637	307.147
Velocity $v_y$ in $\frac{\text{m}}{\text{s}}$	-206.990	1.304	186.548
Specific volume $\xi ={\rho }^{-1}$ in $\frac{{\text{m}}^3}{\text{kg}}$	$0.106$	$0.333$	$1.021$
Molar concentration $\left[C{H}_4\right]$	$0.000$	$0.035$	$0.586$
Molar concentration $\left[O_2\right]$	$0.000$	$0.038$	$0.066$
Molar concentration $\left[C{O}_2\right]$	$0.000$	$0.002$	$0.008$
Molar concentration $\left[H_{2}O\right]$	$0.000$	$0.104$	$0.161$

Table 3. Table 3 : CPU times for two ROMs with time step size Δ ​ t = 1 × 10 − 7 ​ s Δ 𝑡 1 superscript 10 7 s \Delta t=1\times 10^{-7}\text{s} and 30000 time steps.

ROM order	LS solve of (16)	ROM simulation	Reconstructing high-dim. field
$r=24$	2.80s	6.06s	0.04s
$r=29$	6.05s	6.22s	0.05s

Equations

$$\frac{\partial {\vec{q}}_c}{\partial t}+\nabla \cdot (\vec{K}-{\vec{K}}_v)=\vec{S}$$

$${\vec{q}}_c=[\rho \rho v_x\rho v_y\rho e\rho Y_1\dots \rho Y_{n_{\text{sp}}}{]}^{\top },$$

$$e=\sum _{l=1}^{n_{\text{sp}}}{h}_l{Y}_l+\frac{1}{2}\left(v_x^2+v_y^2\right)-\frac{p}{\rho }=h^0-\frac{p}{\rho },$$

$$\vec{K}=\left[\begin{array}{c}\rho v_x\\ \rho v_x^2+p\\ \rho v_x{v}_y\\ \rho v_{x}e+p{v}_x\\ \rho v_x{Y}_1\\ ⋮\\ \rho v_x{Y}_{n_{\text{sp}}}\end{array}\right]\vec{i}+\left[\begin{array}{c}\rho v_y\\ \rho v_x{v}_y\\ \rho v_y^2+p\\ \rho v_{y}e+p{v}_y\\ \rho v_y{Y}_1\\ ⋮\\ \rho v_y{Y}_{n_{\text{sp}}}\end{array}\right]\vec{j},{\vec{K}}_v=\left[\begin{array}{c}0\\ {\tau }_{xx}\\ {\tau }_{yx}\\ {\tau }_{xx}{v}_x+{\tau }_{yx}{v}_y-j_x^q\\ -j_{1,x}^m\\ ⋮\\ -j_{n_{\text{sp}},x}^m\end{array}\right]\vec{i}+\left[\begin{array}{c}0\\ {\tau }_{xy}\\ {\tau }_{yy}\\ {\tau }_{xy}{v}_x+{\tau }_{yy}{v}_y-j_y^q\\ -j_{1,y}^m\\ ⋮\\ -j_{n_{\text{sp}},y}^m\end{array}\right]\vec{j}.$$

$$\tau =\left[\begin{array}{cc}{\tau }_{xx}& {\tau }_{xy}\\ {\tau }_{xy}& {\tau }_{yy}\end{array}\right]=\hat{\mu }\left[\begin{array}{cc}\frac{1}{3}\frac{\partial v_x}{\partial x}& \frac{\partial v_x}{\partial y}+\frac{\partial v_y}{\partial x}\\ \frac{\partial v_y}{\partial x}+\frac{\partial v_x}{\partial y}& \frac{1}{3}\frac{\partial v_y}{\partial y}\end{array}\right],$$

$${\vec{j}}^q={\left[j_x^q{j}_y^q\right]}^{\top }=-\kappa \nabla T+\rho \sum _{l=1}^{n_{\text{sp}}}{D}_l{h}_l\nabla Y_l,$$

$${\vec{j}}_l^m={\left[\begin{array}{c}{j}_{l,x}^m{j}_{l,y}^m\end{array}\right]}^{\top }={\left[\begin{array}{c}\rho D_l\frac{\partial Y_l}{\partial x}\rho D_l\frac{\partial Y_l}{\partial y}\end{array}\right]}^{\top }.$$

$$\vec{S}={\left[0000{\dot{\omega }}_1\dots {\dot{\omega }}_{n_{\text{sp}}}\right]}^{\top }$$

$${\text{CH}}_4+2{\text{O}}_2\to {\text{CO}}_2+2{\text{H}}_2\text{O},$$

$$Y_l=\frac{c_l{M}_l}{\rho },$$

$${\dot{\omega }}_l=M_l\frac{\partial c_l^{\text{reaction}}}{\partial t}={\nu }_l{\Gamma }_r,$$

$${\Gamma }_r=k\prod _{l=1}^{n_{\text{reactant}}}{c}_l^{o_l},$$

$$k=A\exp \left(\frac{-E_a}{R_{\text{u}}T}\right),$$

$$\rho =\frac{p}{RT},$$

$$p_{\text{back}}(t)=p_{\text{back,ref}}[1+A\sin (2\pi ft)],$$

$$\frac{\text{d}\mathbf{q}}{\text{d}t}=\mathbf{G}(\mathbf{q},\mathbf{u}(t)),\mathbf{q}(0)={\mathbf{q}}_0,$$

$$\frac{\text{d}\mathbf{q}}{\text{d}t}=\mathbf{A}\mathbf{q}+\mathbf{H}(\mathbf{q}\otimes \mathbf{q})+\mathbf{C}(\mathbf{q}\otimes \mathbf{q}\otimes \mathbf{q})+\mathbf{B}\mathbf{u}+\mathbf{c}+\text{HOT}.$$

$$\frac{\text{d}\hat{\mathbf{q}}}{\text{d}t}=\hat{\mathbf{A}}\hat{\mathbf{q}}+\hat{\mathbf{H}}(\hat{\mathbf{q}}\otimes \hat{\mathbf{q}})+\hat{\mathbf{C}}(\hat{\mathbf{q}}\otimes \hat{\mathbf{q}}\otimes \hat{\mathbf{q}})+\hat{\mathbf{B}}\mathbf{u}+\hat{\mathbf{c}}+\text{HOT},$$

$$\mathbf{Q}=[{\mathbf{q}}_0\dots {\mathbf{q}}_K]\in {\mathbb{R}}^{d{n}_x\times K},\mathbf{U}=[{\mathbf{u}}_0,\dots ,{\mathbf{u}}_K]\in {\mathbb{R}}^{m\times K},$$

$$\mathbf{Q}=\mathbf{V}\mathbf{\Sigma }{\mathbf{W}}^{\top },$$

$$\hat{\mathbf{Q}}={\mathbf{V}}_r^{\top }\mathbf{Q}=[{\hat{\mathbf{q}}}_0\dots {\hat{\mathbf{q}}}_K]\in {\mathbb{R}}^{r\times K},\dot{\hat{\mathbf{Q}}}=[{\dot{\hat{\mathbf{q}}}}_0{\dot{\hat{\mathbf{q}}}}_1\dots {\dot{\hat{\mathbf{q}}}}_K]\in {\mathbb{R}}^{r\times K},$$

$$\dot{\hat{\mathbf{q}}}=\hat{\mathbf{A}}\hat{\mathbf{q}}+\hat{\mathbf{H}}(\hat{\mathbf{q}}\otimes \hat{\mathbf{q}})+\hat{\mathbf{B}}\mathbf{u}+\hat{\mathbf{c}},$$

$$\underset{\hat{\mathbf{A}}\in {\mathbb{R}}^{r\times r},\hat{\mathbf{H}}\in {\mathbb{R}}^{r\times r^2},\hat{\mathbf{B}}\in {\mathbb{R}}^{r\times m},\hat{\mathbf{c}}\in {\mathbb{R}}^r}{\min }{\Vert {\hat{\mathbf{Q}}}^{\top }{\hat{\mathbf{A}}}^{\top }+(\hat{\mathbf{Q}}\otimes \hat{\mathbf{Q}}{)}^{\top }{\hat{\mathbf{H}}}^{\top }+{\mathbf{U}}^{\top }{\hat{\mathbf{B}}}^{\top }+1_K{\hat{\mathbf{c}}}^{\top }-{\dot{\hat{\mathbf{Q}}}}^{\top }\Vert }_2^2,$$

$$\mathbf{O}=[\hat{\mathbf{A}}\hat{\mathbf{H}}\hat{\mathbf{B}}\hat{\mathbf{c}}]\in {\mathbb{R}}^{r\times (r+r^2+m+1)},$$

$$\mathbf{D}=\left[\begin{array}{c}{\hat{\mathbf{Q}}}^{\top }(\hat{\mathbf{Q}}\otimes \hat{\mathbf{Q}}{)}^{\top }{\mathbf{U}}^{\top }{1}_K\end{array}\right]\in {\mathbb{R}}^{K\times (r+r^2+m+1)},$$

$$\underset{\mathbf{O}\in {\mathbb{R}}^{r\times (r+r^2+m+1)}}{\min }{\Vert \mathbf{D}{\mathbf{O}}^{\top }-{\dot{\hat{\mathbf{Q}}}}^{\top }\Vert }_2^2.$$

$$\underset{{\mathbf{o}}_i\in {\mathbf{R}}^{r+r^2+m+1}}{\min }{\Vert \mathbf{D}{\mathbf{o}}_i-{\mathbf{r}}_i\Vert }_2^2+\lambda {\Vert {\mathbf{P}}_i{\mathbf{o}}_i\Vert }_2^2\text{for }i=1,\dots ,r,$$

$${\vec{q}}_L$$

$$\frac{\partial \xi }{\partial x}=\frac{\partial }{\partial x}\frac{1}{\rho }=-\frac{1}{{\rho }^2}\frac{\partial \rho }{\partial x}=-{\xi }^2\frac{\partial \rho }{\partial x},$$

$$\frac{\partial \xi }{\partial t}=\frac{\partial }{\partial t}\frac{1}{\rho }=-\frac{1}{{\rho }^2}\dot{\rho }={\xi }^2\nabla \cdot \left(\rho v_x\vec{i}+\rho v_y\vec{j}\right)={\xi }^2\left[\frac{\partial \rho }{\partial x}{v}_x+\rho \frac{\partial v_x}{\partial x}\right]+{\xi }^2\left[v_y\frac{\partial \rho }{\partial y}+\rho \frac{\partial v_y}{\partial y}\right].$$

Full text

Nomenclature

1 Introduction

This paper presents an approach to learning low-dimensional surrogate models for a complex, nonlinear, multi-physics, multi-scale dynamical system in form of a multi-species combustion process. The need for repeated model evaluations in optimization, design, uncertainty quantification and control of aerospace systems has driven the development of reduced-order models (ROMs) for applications in aerodynamics [1, 2, 3, 4, 5, 6, 7], reacting flows [8, 9, 10] and combustion [11, 12, 13]. ROMs combine the rich information embedded in high-fidelity simulations with the efficiency of low-dimensional surrogate models; yet, effective and robust ROM methods for nonlinear, multi-scale applications such as combustion have remained an open challenge.

Most existing nonlinear model reduction methods are intrusive—that is, they derive the ROM by projecting the high-fidelity model operators onto a low-dimensional subspace. In doing so, the physics of the problem is embedded in the reduced-order representation. The proper orthogonal decomposition (POD) [14, 15] is the most common way to define the low-dimensional subspace, using the singular value decomposition to identify low-dimensional structure based on training data. For some problems, the projection approach is amenable to rigorous error analysis and structure-preservation guarantees [16, 17, 18, 19], but these rigorous guarantees do not apply to nonlinear, multi-physics, multi-scale models, for which projection-based ROMs remain challenging to implement (due to the need for access to the high-dimensional operators). The compressible flow setting of the combustion process poses numerous problems with respect to stability of the projection-based ROMs, see [20, 21, 22, 23, 24, 25] for several approaches to address this stability problem. Furthermore, ROMs for these problems typically require relatively high dimensionality (and thus high cost) to avoid problems with robustness and stability [12, 26]. For instance, Huang et al. [11] construct two separate ROMs for the same single-injector combustion simulation as presented herein, one that uses POD and the other uses least-squares Petrov-Galerkin projection. That work finds that well over 100 modes are necessary to obtain stable ROMs and sufficient accuracy.

There is increasing attention to non-intrusive model reduction methods (sometimes called black-box or data-driven methods) that learn a model based on training data, without requiring explicit access to the high-fidelity model operators. The non-intrusive philosophy aligns directly with the field of machine learning, where representations such as neural networks have been shown to induce nonlinear model forms that can approximate many physical processes [27]. However, neural networks require a large amount of training data, limiting their utility when the data comes from expensive large-scale partial differential equation (PDE) simulations [28]. Moreover, the parametrization of the learned low-dimensional model is critical to the predictive accuracy and success of the learned model—in particular, it is critical to determining whether the ROM can issue reliable predictions in regimes outside of the training data.

For large-scale PDE models, an important class of non-intrusive learning approaches tackle this challenge of model parametrization by embedding the structure of the problem into the learning formulation. Some approaches use sparse learning techniques to identify PDE model terms that explain the data [29, 30, 31]. Dynamic mode decomposition [32, 33] extracts spectral information of the infinite dimensional linear Koopman operator from observed data of the nonlinear system. This spectral information can then be used to build data-driven predictive models. When the model can be expressed in the form of a dynamical system with polynomial terms, then the learning problem can be formulated as a parameter estimation problem, as in the operator inference approach of [34]. An important advantage of non-intrusive learning approaches is that the user has the flexibility to choose the variables that drive the learning. This opens the way for variable transformations that expose system structure and, in doing so, transform the ROM learning task into a structured form. In some cases, the governing PDEs naturally admit variable transformations that reveal polynomial form, such as the specific volume representation of the Euler equations [35]. More generally, one can introduce new auxiliary variables to the problem—known as lifting—to produce a system that is polynomial in its expanded set of state variables [36, 37, 38]. This allows for a much broader class of nonlinear systems to be learned using the operator inference framework.

In this work, we build on the operator inference framework to learn structured, polynomial ROMs from simulated snapshot data of a single-injector combustion model. While in our case the polynomial model parametrization is a model approximation, we show that the predictive capabilities of the learned ROM are excellent beyond the training data. Our proposed approach follows the steps below, which we develop in detail in the ensuing sections:

We obtain high-dimensional simulation snapshot data for a spatially two-dimensional combustion process from the General Equation and Mesh Solver (GEMS) CFD code [39] developed at Purdue University. The governing equations and combustion problem setup are described in Section 2. 2. 2.

We identify a set of state variables in which many of the terms in the governing equations have quadratic form. We transform the snapshot data to these new state variables, as described in Section 3.3. 3. 3.

We use operator inference to learn a ROM that evolves the combustion dynamics in a low-dimensional subspace. Details of the model learning are given in Section 3 and Section 4.2.

We present numerical results comparing our learned ROMs with GEMS test data in Section 4 and conclude the paper in Section 5.

2 Combustion model

Section 2.1 defines the computational domain under consideration, Section 2.2 presents the governing equations for the combustion model, and Section 2.3 briefly summarizes the numerical implementation. The combustion model follows the implementation of the General Equation and Mesh Solver (GEMS) CFD code [39] and more details can be found in [40]. The GEMS code has been successfully used for rocket engine simulations [11] and in high-pressure gas turbines [41].

2.1 Computational domain

A single-injector combustor as in [42] is shown in Figure 1(a), with the computational domain outlined in red dashed lines. Our domain is a simplified two-dimensional version of the computational domain, shown in Figure 1(b), which also shows the four locations where we monitor the state variables.

2.2 Governing equations

The dynamics of the combustor are governed by the conservation equations for mass, momentum, energy and species mass fractions. For this two-dimensional problem, the conservation equations are

[TABLE]

and they describe the evolution of the conservative variables

[TABLE]

where $\rho$ is the density ($\frac{\text{kg}}{\text{m}^{3}}$), $v_{x}$ and $v_{y}$ are the $x$ and $y$ velocity ($\frac{\text{m}}{\text{s}}$), $e$ is the total energy ($\frac{\text{J}}{\text{m}^{3}}$), and $Y_{l}$ is the $l$th species mass fraction with $l=1,2,\ldots,n_{\textrm{sp}}$ and $n_{\textrm{sp}}$ is defined as the number of chemical species that are included in the model.

The total energy is defined as

[TABLE]

where pressure $p$ is given in (Pa), $h_{l}=h_{l}(T)$ is the enthalpy corresponding to the $l$th species and is a highly nonlinear function of temperature, $T$, and $h^{0}$ is the stagnation enthalpy. The inviscid flux $\vec{K}$ and viscous flux $\vec{K}_{v}$ in Eq. (1) are

[TABLE]

The two-dimensional viscous shear tensor is defined as

[TABLE]

where $\hat{\mu}$ is the mixture viscosity coefficient. The diffusive heat flux vector is defined as

[TABLE]

where $\kappa$ quantifies thermal conductivity and $D_{l}$ is the diffusion coefficient for the $l$th species into the mixture, which is an approximation used to model the multi-component diffusion as the binary diffusion of each species into a mixture. The two terms in the definition of the heat flux (Eq. (3)) represent heat transfer due to conductivity and species diffusion. The diffusive mass flux vector of species $l$ is modeled as

[TABLE]

The source term, $\vec{S}$ in Eq. (1) is

[TABLE]

and is defined by considering a 1-step combustion reaction governed by

[TABLE]

as presented in [43], with $n_{\textrm{sp}}=4$. The corresponding general stoichiometric equation is defined as $0=\sum_{l=1}^{n_{\textrm{sp}}}\nu_{l}\chi_{l},$ where $\chi_{1}=\textrm{CH}_{4},\chi_{2}=\textrm{O}_{2},\chi_{3}=\textrm{CO}_{2}$, $\chi_{4}=$ H2O and $\nu_{l}$ is the net stoichiometric coefficients of each species with $\nu_{1}=-1,\nu_{2}=-2,\nu_{3}=1$ and $\nu_{4}=2$. The molar concentration of the $l$th species is denoted by $c_{l}$. In our case, $l\in\{1,2,3,4\}$, so $c_{1}=[\textrm{CH}_{4}],\ c_{2}=[\textrm{O}_{2}],\ c_{3}=[\textrm{CO}_{2}],$ and $c_{4}=[\textrm{H}_{2}\textrm{O}]$ are the molar concentrations. Here, we use the standard bracket notation $[\cdot]$ to indicate molar concentration of a species. The general relationship between a species molar concentration, $c_{l}$, and a species mass fraction, $Y_{l}$, is

[TABLE]

where $Y_{1}$ is the mass fraction of CH4, $Y_{2}$ is the mass fraction of O2, $Y_{3}$ is the mass fraction of CO2 and $Y_{4}$ is the mass fraction of H2O. The production rate of the $l$th species in the source term $\vec{S}$ in Eq. 4 is modeled as

[TABLE]

where $c_{l}^{\text{reaction}}$ are chemical reaction source terms whose dynamics are described below and $\Gamma_{r}$ is the reaction rate. The molar mass of CH4 is $M_{1}=16.04\frac{\textrm{g}}{\textrm{mol}}$, the molar mass of O2 is $M_{2}=32.0\frac{\textrm{g}}{\textrm{mol}}$, the molar mass of CO2 is $M_{3}=44.01\frac{\textrm{g}}{\textrm{mol}}$, and the molar mass of H2O is $M_{4}=18.0\frac{\textrm{g}}{\textrm{mol}}$.

The reaction rate is approximated by

[TABLE]

where $n_{\textrm{reactant}}=2$ is the number of reactants, $k$ is the rate coefficient and $o_{l}$ is the reaction order of the $l$th reactant. In our case $o_{1}=0.2$ and $o_{2}=1.3$. The rate coefficient, $k$, is described by the Arrhenius equation as

[TABLE]

where $R_{\textrm{u}}=8.314\frac{\text{J}}{\text{mol K}}$ is the universal gas constant, $A=2\times 10^{10}$ is the pre-exponential constant and $E_{a}=2.025\times 10^{5}$ is the energy required to reach a chemical reaction, measured in Joules and referred to as the activation energy. In this work, we use the ideal gas state equation that relates density and pressure to temperature

[TABLE]

where $R=\frac{R_{u}}{M}$ and $M=\left(\sum^{n_{sp}}_{l=1}\left(\frac{Y_{l}}{M_{l}}\right)\right)^{-1}$ is the average molar mass of the mixtures. Thus, we can obtain temperature via $T=\frac{p}{\rho\ R(Y_{l})}$ from the states $\rho,p,Y_{l}$.

At the downstream end of the combustor, we impose a non-reflecting boundary condition while maintaining the chamber pressure via

[TABLE]

where $p_{\textrm{back,ref}}=1.0\times 10^{6}$ Pa, $A=0.1$ and $f=5000$Hz. The top and bottom wall boundary conditions are no-slip conditions, and for the upstream boundary we impose constant mass flow at the inlets.

2.3 Numerical model

GEMS uses the finite volume method to discretize the conservation equations (1). The primitive variables $\vec{q}_{p}=[p\ v_{x}\ v_{y}\ T\ Y_{1}\ldots Y_{n_{\textrm{sp}}}]^{\top}$ are chosen as solution variables in GEMS, since they allow for easier computation of thermal properties and provide more flexibility when extending to complex fluid problems like liquid and supercritical fluids. For a spatial discretization with $n_{x}$ cells, this results in a $dn_{x}$-dimensional system of nonlinear ordinary differential equations (ODEs)

[TABLE]

for $0<t\leq T$, where $d$ is the number of unknowns in the PDE governing equations and here $d=8$ (four flow variables and four species concentrations). In Eq. (10), $\mathbf{q}(t)\in\mathbb{R}^{dn_{x}}$ is the discretized state vector at time $t$ (for GEMS, it is the discretization of the primitive variables $\vec{q}_{p}=[p\ v_{x}\ v_{y}\ T\ Y_{1}\ldots Y_{n_{\textrm{sp}}}]^{\top}$), $\mathbf{q}_{0}$ are the specified initial conditions, and $\frac{\text{d}\mathbf{q}}{\text{d}t}$ is the time derivative of the state vector at time $t$. The $m$ inputs $\mathbf{u}(t)\in\mathbb{R}^{m}$ arise from the time-dependent boundary condition, defined in Eq. (9), applied at the combustor downstream end. The nonlinear function $\mathbf{G}:\mathbb{R}^{dn_{x}}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{dn_{x}}$ maps the discretized states $\mathbf{q}$ and the input $\mathbf{u}$ to the state time derivatives, representing the spatial discretization of the governing equations described in Section 2.2.

Solving these high-dimensional nonlinear ODEs is expensive, motivating the derivation of a ROM that can yield approximate solutions at reduced cost. The nonlinear multi-scale dynamics represented by these equations makes this a challenging task. To maintain computational efficiency in the ROM, state-of-the-art nonlinear model reduction methods combine POD with a sparse interpolation method (often called hyperreduction) by evaluating the nonlinear functions only at a select number of points. For instance, POD together with the discrete empirical interpolation method (DEIM) has had some success, but also encountered problems in combustion applications [12] . Of particular challenge is the need to include a large number of interpolation points in the POD-DEIM approximation, which means that the ROM loses its computational efficiency. Robustness and stability of the POD-DEIM models is also a challenge [26]. In the next section, we present a different approach that uses non-intrusive ROM learning to enable variable transformations that expose system structure. This structure is then exploited in the derivation of the ROM and removes the need for the DEIM approximation.

3 Non-intrusive learning of a combustion reduced model

This section presents our approach to learn ROMs for the unsteady combustion dynamics simulation from GEMS. Section 3.1 writes a general nonlinear system in a form that exposes the underlying structure of the governing equations and shows how projection preserves that structure. Section 3.2 presents the operator inference approach from [34], which learns structured ROM operators from simulation data. Section 3.3 describes variable transformations that lead to the desired polynomial structure for the combustion governing equations presented in Section 2. These transformations yield the structure needed to apply the operator inference approach.

3.1 Projection preserves polynomial structure in the governing equations

Consider a large-scale sytem of nonlinear ODEs written in polynomial form

[TABLE]

Relating this equation to the general nonlinear system in Eq. (10), we see that ${\mathbf{A}}{\mathbf{q}}$ are the terms in $\mathbf{G(\cdot)}$ that are linear in the state $\mathbf{q}$, with ${\mathbf{A}}\in\mathbb{R}^{dn_{x}\times dn_{x}}$; ${\mathbf{H}}({\mathbf{q}}\otimes{\mathbf{q}})$ are the terms in $\mathbf{G(\cdot)}$ that are quadratic in $\mathbf{q}$, with ${\mathbf{H}}\in\mathbb{R}^{dn_{x}\times(dn_{x})^{2}}$; ${\mathbf{C}}({\mathbf{q}}\otimes{\mathbf{q}}\otimes{\mathbf{q}})$ are the terms in $\mathbf{G(\cdot)}$ that are cubic in $\mathbf{q}$, with ${\mathbf{C}}\in\mathbb{R}^{dn_{x}\times(dn_{x})^{3}}$; ${\mathbf{B}}\mathbf{u}$ are the terms in $\mathbf{G(\cdot)}$ that are linear in the input $\mathbf{u}$, with ${\mathbf{B}}\in\mathbb{R}^{dn_{x}\times m}$; and $\mathbf{c}\in\mathbb{R}^{dn_{x}}$ are constant terms in $\mathbf{G(\cdot)}$ that do not depend on state or input. The abbreviation “HOT” in Eq. (11) denotes higher-order terms, and represents terms that are quartic and higher order, as well as any other nonlinear terms that cannot be represented in polynomial form.

We emphasize that we are not (yet) introducing approximations—rather, we are explicitly writing out the discretized equations in the form (11) to expose the system structure that arises from the form of the terms in the governing PDEs. For example, a term such as $\frac{\partial}{\partial x}\rho v_{x}$ in Eq. (1) is linear in the state $\rho v_{x}$, while a term such as $\frac{\partial}{\partial x}\rho v_{x}Y_{1}$ is quadratic in the states $\rho v_{x}$ and $Y_{1}$. Also note that the term $\frac{\partial}{\partial x}\rho v_{x}$ is quadratic in the states $\rho$ and $v_{x}$, highlighting the important point that the structure of the nonlinear model depends on the particular choice of state variables.

A projection-based ROM of Eq. (11) preserves the polynomial structure. Approximating the high-dimensional state $\mathbf{q}$ in a low-dimensional basis $\mathbf{V}\in\mathbb{R}^{dn_{x}\times r}$, with $r\ll dn_{x}$, we write $\mathbf{q}\approx\mathbf{V}\widehat{\mathbf{q}}$. Using a Galerkin projection, this yields the ROM of Eq. (11) as

[TABLE]

where $\widehat{\mathbf{A}}=\mathbf{V}^{\top}\mathbf{A}\mathbf{V}\in\mathbb{R}^{r\times r}$, $\widehat{\mathbf{H}}=\mathbf{V}^{\top}\mathbf{H}(\mathbf{V}\otimes\mathbf{V})\in\mathbb{R}^{r\times r^{2}}$, $\widehat{\mathbf{C}}=\mathbf{V}^{\top}\mathbf{C}(\mathbf{V}\otimes\mathbf{V}\otimes\mathbf{V})\in\mathbb{R}^{r\times r^{3}}$, and $\widehat{\mathbf{B}}=\mathbf{V}^{\top}\mathbf{B}\in\mathbb{R}^{r\times m}$ are the ROM operators corresponding respectively to $\mathbf{A}$, $\mathbf{H}$, $\mathbf{C}$, and $\mathbf{B}$, and $\widehat{\mathbf{c}}=\mathbf{V}^{\top}\mathbf{c}\in\mathbb{R}^{r}$ is a constant vector. We note again that projection preserves polynomial structure, that is, (12) has the same polynomial form as (11), but in the reduced subspace defined by $\mathbf{V}$.

In what follows, we will work with a quadratic system in order to simplify notation. We note that the least squares learning approach described below applies directly to cubic, quartic and all higher-order polynomial terms (although it should be noted that the number of elements in the ROM operators scales with $r^{4}$ for the cubic operator, $r^{5}$ for the quartic operator, etc.). Higher-order terms often exhibit significant block-sparsity that can be exploited in numerical implementations, which limits the growth of computational cost to solve the ROM. For terms in the governing equations that are not in polynomial form (such as terms involving $\frac{1}{\rho}$, and the Arrhenius reaction terms) we discuss in Section 3.3 the introduction of variable transformations and auxiliary variables via the process of lifting [36, 37] to convert these terms to polynomial form.

3.2 Operator inference for learning reduced models

Here we summarize the steps of the operator inference approach from [34]. First, we collect $K$ snapshots of the state by solving the high-fidelity model. We store the snapshots and the inputs used to generate them in the matrices:

[TABLE]

where $\mathbf{u}_{i}\equiv\mathbf{u}(t_{i})$ and $\mathbf{q}_{i}\equiv\mathbf{q}(t_{i})$ with $0=t_{0}<t_{1}<\dots<t_{K}=T$. In general, $dn_{x}\gg K$, so the matrix $\mathbf{Q}$ is tall and skinny. Second, we identify the low-dimensional subspace in which we will learn the ROM. In this work, we use the POD to define the low-dimensional subspace, by computing the singular value decomposition of the snapshot matrix

[TABLE]

where $\mathbf{V}\in\mathbb{R}^{dn_{x}\times K}$, $\boldsymbol{\Sigma}\in\mathbb{R}^{K\times K}$ and $\mathbf{W}\in\mathbb{R}^{K\times K}$. The $r\ll dn_{x}$ dimensional POD basis, $\mathbf{V}_{r}=[\mathbf{v}_{1},...,\mathbf{v}_{r}]$, is given by the first $r$ columns of $\mathbf{V}$. Third, we project the state snapshot data onto the POD subspace spanned by the columns of $\mathbf{V}_{r}$ and obtain the reduced snapshot matrices

[TABLE]

where the columns of $\dot{\widehat{\mathbf{Q}}}$ are computed from $\widehat{\mathbf{Q}}$ using any time derivative approximation (see, e.g., [44, 45, 46]), or can be obtained—if available—by collecting and projecting snapshots of $\mathbf{G}(\mathbf{q}_{i},\mathbf{u}_{i})$.

Operator inference solves a least squares problem to find the reduced operators that yield the ROM that best matches the projected snapshot data in a minimum residual sense. For the quadratic ROM

[TABLE]

operator inference solves the least squares problem

[TABLE]

where $\mathbf{1}_{K}\in\mathbb{R}^{K}$ is the length $K$ column vector with all entries set to unity. Note that this least squares problem is linear in the coefficients of the unknown ROM operators $\widehat{\mathbf{A}}$, $\widehat{\mathbf{H}}$, $\widehat{\mathbf{B}}$ and $\widehat{\mathbf{c}}$. Also note that the operator inference approach permits us to compute the ROM operators $\widehat{\mathbf{A}}$, $\widehat{\mathbf{H}}$, $\widehat{\mathbf{B}}$ and $\widehat{\mathbf{c}}$ without needing explicit access to the original high-dimensional operators ${\mathbf{A}}$, ${\mathbf{H}}$, ${\mathbf{B}}$ and ${\mathbf{c}}$.

We combine the unknown operators of Eq. (13) in the matrix

[TABLE]

and the known low-dimensional data in the data matrix

[TABLE]

and then solve the minimization problem

[TABLE]

For $K>r+r^{2}+m+1$ this overdetermined linear least-squares problem has a unique solution [47, Sec. 5.3]. It was proven in [34] that Eq. (15) can be written as $r$ independent least squares problems of the form $\min_{\mathbf{o}_{i}\in\mathbb{R}^{r+r^{2}+m+1}}\left\|\mathbf{D}\mathbf{o}_{i}-\mathbf{r}_{i}\right\|_{2}^{2}$, for $i=1,\dots,r$, where $\mathbf{o}_{i}$ is a column of $\mathbf{O}^{\top}$ (row of $\mathbf{O}$) and $\mathbf{r}_{i}$ is a column of $\dot{\widehat{\mathbf{Q}}}^{\top}$. This makes the operator inference approach efficient and scalable.

Regularization becomes necessary to avoid overfitting and to infer operators that produce a stable ROM. In this work, we use an L2 regularization penalty on the off-diagonal elements of the operator $\widehat{\mathbf{A}}$ and on all elements of the remaining operators. With this regularization, our least squares problem becomes

[TABLE]

where $\lambda$ is the regularization parameter and $\mathbf{P}_{i}$ is the $r+r^{2}+m+1$ identity matrix with the $i$th diagonal set to zero so that we avoid regularizing the diagonal elements of $\mathbf{A}$. It should be noted that the regularization parameter, $\lambda$, is problem specific and should be chosen accordingly. In Section 4.2, we discuss details of the operator inference implementation, a method for selecting $\lambda$, and the removal of redundant terms in the least squares problem in Eq. (16).

3.3 A structure-exploiting ROM learning formulation for GEMS

A key contribution of this work is to recognize that the non-intrusive operator inference approach gives us complete flexibility in the set of physical variables we work with to define the ROM. We can identify choices of physical variables that expose the desired polynomial structure in the governing equations, and then extract snapshots for those variables by applying transformations to the snapshot data—we do not need to make any modifications to the high-fidelity CFD simulation model itself. In theory, a classical intrusive ROM approach could work with transformed variables (e.g., in the work of [24]); however, this would involve rewriting the high-fidelity simulator, a task that would be not only time-consuming but also fraught with mathematical pitfalls, especially for unusual choices of variables. This is where the data-driven perspective of machine learning becomes extremely valuable.

The Euler equations admit a quadratic representation in the specific volume variables; in that case, a transformation of the snapshots from conservative (or primitive) variables to specific volume variables can be exploited to create quadratic ROMs [35]. Other PDEs may not admit polynomial structure via such straightforward transformations, but the process of lifting the equations via the introduction of new auxiliary variables can produce a set of coordinates in which the governing equations become polynomial in the lifted state [36, 37, 38]. For example, the tubular reactor example of [37] includes Arrhenius-type reaction terms similar to those in Eq. (7). The introduction of auxiliary variables permits the governing equations to be written equivalently with quartic nonlinearity in the lifted variables.555The Arrhenius reaction terms can be lifted further to quadratic form, but then require the inclusion of algebraic constraints, which makes the model reduction task more difficult, see [37].

Lifting to polynomial form for the GEMS equations described in Section 2 is made difficult by several of the terms, in particular through some of the gas thermal properties such as the nonlinear dependence of enthalpy on temperature. A complete lifting that converts all equations to a polynomial form is possible, but would require the introduction of a large number of auxiliary variables and would also result in the introduction of some algebraic equations. However, as the analysis below shows, the GEMS governing equations admit a transformation for which many terms in the governing equations take polynomial form when we use the variables

[TABLE]

Here $\xi=\frac{1}{\rho}$ is the specific volume, and recall that $c_{1}=[\textrm{CH}_{4}],\ c_{2}=[\textrm{O}_{2}],\ c_{3}=[\textrm{CO}_{2}],$ and $c_{4}=[\textrm{H}_{2}\textrm{O}]$ are the molar concentrations with $c_{l}=\frac{\rho Y_{l}}{M_{l}}$.

Below, we derive the governing PDEs for specific volume $\xi$ and velocities $v_{x},v_{y}$. These three governing PDEs all turn out to be quadratic in the learning variables $\vec{q}_{L}$. In Appendix A we present the lifting transformations for the source term dynamics $c_{l}^{\text{reaction}}$ in the vector $\vec{S}$ in Eq. (4). In Appendix B we derive the equations governing the pressure $p$ and the species molar concentrations $c_{i}$. These equations have some terms that are not polynomial in the chosen learning variables $\vec{q}_{L}$.

To keep notation clean in application of the chain rule, let the conservative variables be denoted as $g_{1}=\rho,\ g_{2}=\rho v_{x},\ g_{3}=\rho v_{y},\ g_{4}=\rho e,\ g_{5}=\rho Y_{1},\ g_{6}=\rho Y_{2},\ g_{7}=\rho Y_{3},\ g_{8}=\rho Y_{4}$. Throughout, we frequently use the relationship

[TABLE]

and similarly for $\frac{\partial\xi}{\partial y}$. Note also that we are assuming the existence of these partial derivatives, that is, we do not consider the case of problems with discontinuities.

Specific Volume $\xi=1/\rho$. We use the constitutive relationship for the density $\rho$ in Eq. (1) in the derivation:

[TABLE]

Inserting Eq. (18) into the above, we obtain

[TABLE]

which is quadratic in the learning variables $\xi,v_{x},v_{y}$.

Velocities $v_{x},v_{y}$: We have

[TABLE]

and from Eq. (1) we have that $\dot{g}_{1}=\dot{\rho}=-\nabla\cdot\left(\rho v_{x}\vec{i}+\rho v_{y}\vec{j}\right)$ as well as $\dot{g}_{2}=\dot{(\rho v_{x})}=\nabla\cdot\left(-(\rho v_{x}^{2}+p)\vec{i}-(\rho v_{x}v_{y})\vec{j}+\tau_{xx}\vec{i}+\tau_{xy}\vec{j}\right)$. Thus, we obtain

[TABLE]

and we get a similar expression for $\frac{\partial v_{y}}{\partial t}$. Both dynamics are quadratic in the learning variables $p,v_{x},v_{y},\xi$.

As noted above, Appendix A and Appendix B present the derivations for the chemical source terms, pressure, and chemical species.

4 Numerical Results

We now apply the variable transformations and operator inference framework to learn a predictive ROM from GEMS high-fidelity combustion simulation data.666Code for the operator inference framework is available at https://test.pypi.org/project/operator-inference/ in Python and https://github.com/elizqian/operator-inference in Matlab. Section 4.1 describes the problem setup and GEMS dataset. Section 4.2 discusses implementation details and Section 4.3 presents our numerical results. Additional numerical results can be found in [48].

4.1 GEMS Dataset

The computational domain shown in Figure 1(b) is discretized with $n_{x}=38523$ spatial discretization points. Each CFD state solution thus has dimension $dn_{x}=308184$. The problem considered here has fuel and oxidizer input streams with constant mass flow rates of $5.0\frac{\text{kg}}{\text{s}}$ and $0.37\frac{\text{kg}}{\text{s}}$, respectively. The fuel is composed of gaseous methane and the oxidizer is 42% gaseous O2 and 58% gaseous H2O, as described in [11]. The forcing input Eq. (9) is applied at the right side of the domain. For this simulation, the resulting Reynolds number is about 10,000, defined as $Re=\frac{\rho{v_{x}L}}{\mu}$ where the density $\rho$, horizontal velocity $v_{x}$ and viscosity $\mu$ are evaluated at the inlet of the oxidizer post ($x=-0.04$m in Fig. 1(b)), and the characteristic length $L$ is defined as the height of the oxidizer inlet. The highest Mach number is $\approx 0.25$ and is evaluated inside the oxidizer post (from $x=-0.04$m to [math]m in Fig. 1(b)).

To generate training data, GEMS is simulated for a time duration of 1ms with a time step size of $\Delta t=1\times 10^{-7}$s resulting in $K=10000$ snapshots. The GEMS output is transformed to the variables given in Eq. (17). The recorded snapshot matrix is thus

[TABLE]

Our numerical experiments were parallelized on a cluster with two computing nodes. Each node has two 10-core Intel Xeon-E5 processors (20 cores per node) and 128 GB RAM. The training data generation took approximately 200h in CPU time for the 1ms, 10000 snapshots of high-fidelity CFD data.

The range of variable values for the training data is shown in Table 2. Note that the data covers a wide range of scales. Pressure is of the order $10^{6}$ while species concentrations can be as low as 10*-12*. This large scaling difference presents a challenge when learning models from data. To deal with the numerical issues related to large differences in scaling and small species concentrations and velocities, we scale each variable to the interval $[-1,1]$. Variables are scaled before computing the POD basis and projecting the data.

To obtain snapshots of the projected state time derivative, we approximate the derivative with a five-point approximation $\dot{\widehat{\mathbf{q}}}_{i}=(-{\widehat{\mathbf{q}}}_{i+2}+8{\widehat{\mathbf{q}}}_{i+1}-8{\widehat{\mathbf{q}}}_{i-1}+{\widehat{\mathbf{q}}}_{i-2})/(12\Delta t)$. This approximation is fourth order accurate. The first two and last two time derivatives are computed using first-order forward and backward Euler approximations, respectively.

4.2 Learning a quadratic reduced-order model

To learn the operators of the quadratic ROM, we solve the regularized least squares problem shown in Eq. (16). We are using numpy’s least squares solve numpy.linalg.lstsq. The algorithm is based on the LAPACK routine xGELSD. That routine is based on the SVD, which typically provides a stable implementation. In what follows next, we describe several important implementation details.

4.2.1 Singular value decomposition implementation and POD basis selection

Due to the large size of this dataset, we implement the randomized singular value decomposition algorithm, introduced in [49], to compute the leading 500 singular values and vectors of the snapshot matrix. The randomized singular value decomposition algorithm can be implemented in a scalable way for large datasets as the data does not have to be read into single memory all at once. The POD basis is chosen as the $r$ leading left singular vectors. The dimension $r$ is typically chosen so that the cumulative energy contained in the subspace is greater than a user specified tolerance $\epsilon$, i.e.,

[TABLE]

where $\sigma_{k}^{2}$ are the squared singular values of the data matrix $\mathbf{Q}$. To guide the choice of $r$ we also use the relative projection error

[TABLE]

4.2.2 Removing redundant terms in least squares problem

There are redundant terms that arise in the Kronecker product $\widehat{\mathbf{Q}}\otimes\widehat{\mathbf{Q}}$ in Eq. (14), which can cause the least squares problem to become ill-posed. To see this, consider $q\otimes q$ with $q=[q_{1}\ q_{2}]^{T}$, i.e., we have $q\otimes q=[q_{1}^{2}\ \ q_{2}q_{1}\ \ q_{1}q_{2}\ \ q_{2}^{2}]$. Since we know where the repeated terms occur in the product we merely need to remove the redundant (repeated) terms before we solve the least-squares problem.

Thus, the Kronecker product is replaced with the term

[TABLE]

where $s=\frac{r(r+1)}{2}$. Each vector $\widehat{\mathbf{q}}^{2}_{j}$ is defined, according to [34], as

[TABLE]

and $\widehat{q}_{i,j}$ is the $i$th element of the vector $\widehat{\mathbf{q}}_{j}$. Now, instead of learning the operator $\widehat{\mathbf{H}}\in\mathbb{R}^{r\times r^{2}}$, the operator $\widehat{\mathbf{F}}\in\mathbb{R}^{r\times s}$ is learned, which satisfies the equivalent least squares problem

[TABLE]

The least squares problem is again of the form (15) but we use the data matrix $\mathbf{D}=\begin{bmatrix}\widehat{\mathbf{Q}}^{\top}\quad(\widehat{\mathbf{Q}}^{2})^{\top}\quad\mathbf{U}^{\top}\quad\mathbf{1}_{K}\end{bmatrix}\in\mathbb{R}^{K\times(r+s+m+1)}$ and solve for the operators $\mathbf{O}=[\widehat{\mathbf{A}}\quad\widehat{\mathbf{F}}\quad\widehat{\mathbf{B}}\quad\hat{\mathbf{c}}]\in\mathbb{R}^{r\times(r+s+m+1)}$. Once we have solved for the operator $\widehat{\mathbf{F}}$ we can easily transform it to obtain $\widehat{\mathbf{H}}$.

4.2.3 Regularization

We use an L2-regularization (also known as Tikhonov regularization or ridge regression) to solve the operator inference problem, as shown in Eq. (20). The regularization term introduces a trade-off between operators that fit the data well and operators with small values. This regularization is used to avoid overfitting to the data, which in this setting causes our learned ROMs to be unstable ( solution blow up in finite time). The regularization parameter $\lambda$ affects the performance of this algorithm—we require enough regularization to avoid overfitting, but if $\lambda$ is too large, the data will be poorly fit.

To help determine appropriate values of $\lambda$, we consider the “L-curve" criterion discussed in [50]. The L-curve is a way of visualizing the effects of different values of $\lambda$ on the norm of the residual (data fit) against the norm of the solution. The L-curve criterion recommends choosing a value for $\lambda$ that lies in the corner of the curve, nearest the origin. In our numerical experiments, we compute the L-curve to help determine appropriate values for $\lambda$.

Regularization also helps to reduce the condition number of the regularized least squares data matrix, $[\mathbf{D}\ \ \lambda\mathbf{P}]^{\top}$. In Figure 2 we show the condition number of the data matrix $\mathbf{D}$ from Eq. (20) (with the data already scaled to [-1,1]) when we include $2500,5000,7500$ and $10000$ snapshots in the training set. The condition number is quite large for this application, yet it decreases as we add more training data. This effect is due to the fact that as we add more training data, the first $r$ POD basis vectors become richer. This in turn means that the projected data in $\mathbf{D}$ are richer for a given dimension of the POD basis. Our numerical experiments reinforced this finding by confirming that $10000$ training snapshots were required to achieve a sufficiently rich POD basis to obtain accurate ROMs, an indication of the complexity of the combustor dynamics. While the condition number of the least-squares problem is high, for this example, it remains at manageable levels. However, if the condition number gets much larger, a careful treatment and consideration of proper data sampling and regularization will be required.

4.2.4 Error measures

To evaluate the ROM performance, we define appropriate error metrics for each variable. The error is computed after the ROM solutions have been reconstructed back into the full CFD model dimension and scaled back to the original variable ranges. Recall Table 2, which showed the range of values for each variable. Below we provide details on how error is computed for each variable:

• •

For pressure and temperature, the values are always positive and well above zero, so we use a standard relative error, defined as

[TABLE]

where $\boldsymbol{\zeta}\in\{\mathbf{p},\mathbf{T}\}$ and $\boldsymbol{\zeta}_{\cdot,i}$ denotes the $i$th entry of a vector, and $\boldsymbol{\zeta}_{\textrm{CFD}}$ is the CFD solution variable and $\boldsymbol{\zeta}_{\textrm{ROM}}$ is the solution variable obtained from the ROM simulation.

• •

Due to the small values of species concentrations (on the order of $10^{-12}$), dividing by the true value can skew a small error. Similarly, velocities range from positive to negative, including zero. Thus, for species concentrations and velocities, we use a normalized absolute error, defined as

[TABLE]

where $\boldsymbol{\zeta}\in\{\mathbf{v}_{x},\mathbf{v}_{y},[\mathbf{CH}_{4}],[\mathbf{O}_{2}],[\mathbf{CO}_{2}],[\mathbf{H}_{2}\mathbf{O}]\}$ and $\max_{l}\left(|\boldsymbol{\zeta}_{\textrm{CFD,l}}|\right)$ denotes the maximum entry of $|\boldsymbol{\zeta}_{\textrm{CFD}}|$, i.e., the maximum absolute value over the discretized spatial domain.

4.3 Learned reduced model performance

The given $K=10000$ snapshots representing 1ms of GEMS simulation data (scaled to [-1,1]) are used to learn the ROM. We are also given another 2ms of testing data at the monitor locations shown in Figure 1(b), which we use to assess the predictive capabilities of our learned ROMs beyond the range of training data.

The cumulative energy of the singular values of the snapshot data matrix $\mathbf{Q}$ is shown in Figure 3. The singular values that correspond to a cumulative energy of 98.5% and 99% are indicated with the red triangle and green square, respectively. We use basis sizes of $r=24$, capturing 98.5% of the total energy, and $r=29$, capturing 99% of the total energy.

In Figures 4(a), 4(b) we show the L-curve for each basis size. The L-curve for a basis of $r=24$ is somewhat uninformative in this case. The regularization parameters chosen were $1.0$E+05 and $\lambda=3.0$E+05, which lie in the upper part of the L-curve and are therefore a recommended choice by the L-curve criterion discussed above. For a basis of size of $r=29$, the L-curve indicates a regularization parameter around $\lambda=3.0$E+04. Stable systems are produced for $\lambda=3.0$E+04 and $5.0$E+04.

We simulate the learned ROM for the two model sizes $r=24$ and $r=29$ with the same initial value and time step size ($\Delta t=1\times 10^{-7}$s) as those used to generate the training set. Since the ROM was constructed from data scaled to [-1,1], we solve the ROM in the (scaled) subspace, and then reconstruct the dimensional quantities by reversing the scaling. Figures 5 and 6 compare the time trace of pressure computed by GEMS (our “truth” data) with the ROM predictions for 30000 time steps (10000 time steps used for training, 20000 time steps are pure prediction for the ROM) at the cell located at $(0.0,0.0225)$ in the domain (denoted as monitor location 1 in Figure 1(b)). The performance of the ROM on the training data (first 1ms of data) is good in both cases, although the $r=29$ ROM (Figure 6) is more accurate. For test data predictions beyond the training data, both ROMs yield accurate phase predictions and pressure oscillation amplitudes that are good approximations of the truth.

A Galerkin-projection-based POD method was applied to this GEMS model in Ref. [11]. The authors there found that a large number of modes ($r>100$) was necessary to obtain stable ROMs. Moreover, comparable accuracy, e.g., for the pressure probe predictions shown here, was only achieved with $r=200$ POD modes for the ROM (c.f. Fig. 15 in Ref. [11] with Figures 5–6 herein; both are recording pressure at the same probe location). Consequently, while our learned ROM does not resolve the full flow physics, we do obtain good predictability in time at much lower ROM dimension than in the classical Galerkin-POD approach used in Ref. [11]. The key innovation leading to this improvement is our use of variable transformations to build the ROM over a space for which the transformed governing equations have more polynomial structure. The non-intrusive operator inference approach is an enabler that makes these variable transformations practical from an implementation perspective, since the transformations are applied only to the snapshot data and not to the CFD model itself.

We also compute the average error of each variable over the entire domain at the last time step of the training set (the 10000th time step). The normalized absolute error, defined in Eq. (22), is shown for each species and for $x$ and $y$ velocity in Figure 7. This figure also shows the relative error, defined in Eq. (21), for pressure and temperature. These plots show that overall, the error is decreasing with an increasing basis size. At a basis size of 18, 22, and 24 the systems were unstable ( solution blow up in finite time) and so these basis sizes are excluded from the figure. The cause of this , and the non-monotone error decay, may be due to the fact that the same regularization parameter was used for each of these, $\lambda=3.0$E+04. Ideally, one would pick a parameter specific for the basis size , but here we used the same regularization parameter in order to give a fairer assessment of the ROM performance without using manual tuning to optimize our results. We note, however, that even with manual tuning of the regularization parameters we are not guaranteed monotone state-error decay for strongly nonlinear dynamical system ROMs.

In Figure 8, we show the integrated species concentrations over time. To compute these, at each time step in our simulation, we sum all elements of a species vector. This measure monitors whether our ROM conserves species mass, an important feature of a physically meaningful simulation. As the discretization of the high-fidelity model becomes finer, point-wise error may become large and misleading if the mass is shifted slightly into the neighboring cells. The integrated species concentration complements the evaluation of point-wise errors and provides a global view of the error in the domain. While $\text{CH}_{4}$ conservation is tracked well qualitatively by the ROM, it does show the largest deviation out of the four species. This is a result of $\text{CH}_{4}$ having the sharpest gradients ($\text{CH}_{4}$ concentration ranging from 0 to 1) compared to the other three species, see also Figures 13–16 below, and their respective color bars.

We also compare the state variables over the entire domain predicted by the learned ROM with the given GEMS data at the last time step $K=10000$ (which corresponds to $t=0.0159999$s). We provide the true field, the ROM-predicted field, and an error field for each variable in Figures 9–16. Again, for pressure and temperature, we use a relative error from Eq. (21). For $x$ and $y$ velocity and for species molar concentrations we use a normalized absolute error from Eq. (22). The plots show that the ROM predictions are, as expected, not perfect, but indeed they capture well the overall structure and many details of the solution fields.

Table 3 shows timing results for the ROM generation and simulation, as performed using python 3.6.4 on a dual-core Intel i5 processor with 2.3 GHz and 16GB RAM. We report the following CPU times: solving the operator inference least squares problem (20); ROM runtime for the two different basis sizes for 3ms of real time prediction; and reconstruction of the high-dimensional, unscaled combustion variables. The latter is required since the ROM is evolving the dynamics in the $[-1,1]$ scaled variables, and thus a post-processing step is required to obtain the true magnitudes of the variables. In comparison with the approximately 200h of CPU time needed on a 40-core architecture (more details in Section 4.1) to compute the first 10000 snapshots of data, the ROMs provide five to six orders of magnitude in computational speedup.

5 Conclusion

Operator inference is a data-driven method for learning reduced-order models (ROMs) of dynamical systems with polynomial structure. This paper demonstrates how variable transformations can expose quadratic structure in the nonlinear system of partial differential equations describing a complex combustion model. This quadratic structure is preserved under projection, providing the mathematical justification for learning a quadratic ROM using operator inference. An important feature of the approach is that the learning of the ROM is non-intrusive—it requires state solutions generated by running the high-fidelity combustion model, but it does not require access to the discretized operators of the governing equations. This is important because it means that the variable transformations can be applied as a post-processing step to the simulation data set, but no intrusive modifications are needed to the high-fidelity code. While the quadratic model form is an approximation for this particular application problem, the numerical results show that the learned quadratic ROM can predict relevant quantities of interest and can also conserve species accurately. Many nonlinear equations in scientific and engineering applications admit variable transformations that expose polynomial structure. This combined with the non-intrusive nature of the approach make it a viable option for deriving ROMs for complex nonlinear applications where traditional intrusive model reduction is impractical and/or unreliable. While we improved ROM stability through the presented regularization of the least-squares problem, some of the ROMs were unstable. Future work thus includes devising alternative approaches (see e.g., [23, 24]) and developing theory to guarantee stable learned ROMs.

Appendix A: Lifting chemical source terms

The chemical source terms in $\vec{S}$ in Eq. (4) are $\dot{\omega}_{l}=M_{l}\dot{c}_{l}^{\text{reaction}}$, see Eq. (6). The dynamics for the source terms $\dot{c}_{l}^{\text{reaction}}$ are given by

[TABLE]

To lift this system to polynomial form, we introduce the auxiliary variables

[TABLE]

The source term dynamics (23) are then cubic in $w_{1},w_{3},w_{5}$:

[TABLE]

We next derive the dynamics for the auxiliary variables $w_{1},\ldots,w_{6}$. For instance, we have $\dot{w}_{1}=0.2c_{1}^{0.2-1}\dot{c_{1}}=0.2w_{1}w_{2}\dot{c_{1}}$. Similarly, we obtain the system of auxiliary dynamics:

[TABLE]

The dynamics of the lifted variables are quintic in the variables $w_{1},\ldots,w_{6}$. If we further include an additional auxiliary variable $w_{7}=w_{1}w_{3}w_{5}$, then we obtain the system of equations

[TABLE]

The temperature $T=\frac{p\xi}{R(Y_{l})}$ can be obtained from the states $\xi,p,Y_{l}$ via the ideal gas relationship in Eq. (8).

Appendix B: Equations for pressure and chemical species

Here, we give the complete governing equation for the species molar concentrations $c_{l}$ and pressure $p$. Recall from Section 3.3 the notation used to denote the conservative variables $g_{i}$: $g_{1}=\rho,\ g_{2}=\rho v_{x},\ g_{3}=\rho v_{y},\ g_{4}=\rho e,\ g_{5}=\rho Y_{1},\ g_{6}=\rho Y_{2},\ g_{7}=\rho Y_{3},\ g_{8}=\rho Y_{4}$.

Species molar concentrations $c_{l}$: For the species molar concentration dynamics, we use the relationship $c_{l}=\frac{\rho Y_{l}}{M_{l}}$ from Eq. (5), where the constants $M_{1},\ldots,M_{4}$ are molar masses. We obtain for $l=1,2,\ldots,n_{\text{sp}}$:

[TABLE]

Note that the chemical source terms $\dot{c}_{l}^{\text{reaction}}$ were given in Appendix A, and from Eq. (24) we see that the $\dot{c}_{l}$ are cubic in the auxiliary lifted states $w_{1},\ldots,w_{6}$. The divergence term is

[TABLE]

and since $Y_{l}=M_{l}c_{l}\xi$, we have that

[TABLE]

Overall, we have that

[TABLE]

which is quadratic in the learning variables $v_{x},v_{y},c_{l}$ with the exception of the terms $\rho c_{l}\frac{\partial\xi}{\partial x}$ and $\rho c_{l}\frac{\partial\xi}{\partial y}$. We note that if $\rho$ were included as a lifted variable (in addition to $\xi$), these terms would become quadratic in the lifted state.

Pressure $p$: We start with the energy equation (2). By multiplying with density $\rho$, we have $\rho e=\rho h^{0}-p$, so from the conservation equation (1) for $\rho e$ we obtain

[TABLE]

This directly gives an equation for the time evolution of pressure:

[TABLE]

Moreover, per definition of $h^{0}$ in Eq. (2) and with $c_{l}=\frac{\rho Y_{l}}{M_{l}}$ we have that

[TABLE]

Overall, we have that

[TABLE]

This equation remains nonlinear in our chosen learning variables $\vec{q}_{L}$. In particular, the enthalpies $h_{l}=h_{l}(T)$ and their time derivatives are nonlinear functions of temperature. The other terms show some polynomial structure; for example, in Section 3.3 we showed that $\frac{\partial(v_{x}+v_{y})}{\partial t}$ is quadratic in the learning state variables $p,v_{x},v_{y},\xi$. However, to write this equation exactly in a polynomial form would require introducing a large number of auxiliary variables along with their corresponding dynamics. We have instead chosen to introduce an approximation by learning a ROM in the variables $\vec{q}_{L}$ with quadratic form.

Funding Sources

This work has been supported in part by the Air Force Center of Excellence on Multi-Fidelity Modeling of Rocket Combustor Dynamics award FA9550-17-1-0195, and the Air Force Office of Scientific Research (AFOSR) MURI on managing multiple information sources of multi-physics systems award numbers FA9550-15-1-0038 and FA9550-18-1-0023.