Data-driven reduced-order models via regularised Operator Inference for a single-injector combustion process

ABSTRACT

This paper derives predictive reduced-order models for rocket engine combustion dynamics via Operator Inference, a scientific machine learning approach that blends data-driven learning with physics-based modelling. The non-intrusive nature of the approach enables variable transformations that expose system structure. The specific contribution of this paper is to advance the formulation robustness and algorithmic scalability of the Operator Inference approach. Regularisation is introduced to the formulation to avoid over-fitting. The task of determining an optimal regularisation is posed as an optimisation problem that balances training error and stability of long-time integration dynamics. A scalable algorithm and open-source implementation are presented, then demonstrated for a single-injector rocket combustion example. This example exhibits rich dynamics that are difficult to capture with state-of-the-art reduced models. With appropriate regularisation and an informed selection of learning variables, the reduced-order models exhibit high accuracy in re-predicting the training regime and acceptable accuracy in predicting future dynamics, while achieving close to a million times speedup in computational cost. When compared to a state-of-the-art model reduction method, the Operator Inference models provide the same or better accuracy at approximately one thousandth of the computational cost.

KEYWORDS:

• Model reduction
• Operator Inference
• scientific machine learning
• combustion
• data-driven reduced model

Introduction

The emerging field of scientific machine learning brings together the perspectives of physics-based modelling and data-driven learning. In the field of fluid dynamics, physics-based modelling and simulation have played a critical role in advancing scientific discovery and driving engineering innovation in domains as diverse as biomedical engineering (Yin et al. 2010; Nordsletten et al. 2011), geothermal modelling (O'Sullivan et al. 2001; Cui et al. 2011), and aerospace (Spalart and Venkatakrishnan 2016). These advances are based on decades of mathematical and algorithmic developments in computational fluid dynamics (CFD). Scientific machine learning builds upon these rigorous physics-based foundations while seeking to exploit the flexibility and expressive modelling capabilities of machine learning (Baker et al. 2019). This paper presents a scientific machine learning approach that blends data-driven learning with the theoretical foundations of physics-based model reduction. This creates the capability to learn predictive reduced-order models (ROMs) that provide approximate predictions of complex physical phenomena while exhibiting several orders of magnitude computational speedup over CFD.

Projection-based model reduction considers the class of problems for which the governing equations are known and for which we have a high-fidelity (e.g. CFD) model (Antoulas 2005; Benner et al. 2015). The goal is to derive a ROM that has lower complexity and yields accurate solutions with reduced computation time. Projection-based approaches define a low-dimensional manifold on which the dominant dynamics evolve. This manifold may be defined as a function of the operators of the high-fidelity model, as in interpolatory methods that employ a Krylov subspace (Bai 2002; Freund 2003), or it may be determined empirically from representative high-fidelity simulation data, as in the proper orthogonal decomposition (POD) (Lumley 1967; Sirovich 1987; Berkooz et al. 1993). The POD has been particularly successful in fluid dynamics, dating back to the early applications in unsteady flows and turbulence modelling (Sirovich 1987; Deane et al. 1991; Gatski and Glauser 1992), and in unsteady fluid-structure interaction (Dowell and Hall 2001).

Model reduction methods have advanced to include error estimation (Veroy et al. 2003; Grepl and Patera 2005; Veroy and Patera 2005; Rozza et al. 2008) and to address parametric and nonlinear problems (Barrault et al. 2004; Astrid et al. 2008; Chaturantabut and Sorensen 2010; Carlberg et al. 2013), yet the intrusive nature of the methods has limited their impact in practical applications. When legacy or commercial codes are used, as is often the case for CFD applications, it can be difficult or impossible to implement classical projection-based model reduction. Black-box surrogate modelling instead derives the ROM by fitting to simulation data; such methods include response surfaces and Gaussian process models, long used in engineering, as well as machine learning surrogate models. These methods are powerful and often yield good results, but since the approximations are based on generic data-fit representations, they are not equipped with the guarantees (e.g. stability guarantees, error estimators) that accompany projection-based ROMs. Nonlinear system identification techniques seek to illuminate the black box by discovering the underlying physics of a system from data (Brunton et al. 2016). However, when the governing dynamics are known and simulation data are available, reduced models may be directly tailored to the specific dynamics without access to the details of the large-scale CFD code.

This paper presents a non-intrusive alternative to black-box surrogate modelling. We use the Operator Inference method of Peherstorfer and Willcox (2016) to learn a ROM from simulation data; the structure of the ROM is defined by known governing equations combined with the theory of projection-based model reduction. Our approach may be termed ‘glass-box modelling’, which we define as the situation where the form of the targeted dynamics is known (here via the partial differential equations that define the problem of interest), but we do not have internal access to the CFD code that produces the simulation data. That is, we know what dynamics to expect, but we may only calibrate our models using outputs of the full-order CFD model. This glass-box setting is in contrast to black-box modelling approaches which do not exploit knowledge of the governing equations. We build on our prior work in Swischuk et al. (2020) by formally introducing regularisation to the Operator Inference approach. Regularisation is critical to avoid overfitting for problems with complex dynamics, as is the case for the combustion example considered here. A second contribution of this paper is a scalable implementation of the approach, which is available via an open-source implementation.

Methodology

This section begins with an overview of the Operator Inference approach. We then augment Operator Inference with a new regularisation formulation, posed as an optimisation problem, and present a complete algorithm for regularisation selection and model learning. The section concludes with a scalable implementation of the algorithm, which can then be applied to CFD problems of high dimension.

Operator Inference

We target problems governed by systems of nonlinear partial differential equations. Consider the governing equations of the system of interest written, after spatial discretisation, in semi-discrete form (1) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{q}\left(t\right)=\mathbf{c}+\mathbf{A}\mathbf{q}\left(t\right)+\mathbf{H}\left(\mathbf{q}\right(t)\otimes \mathbf{q}(t\left)\right)+\mathbf{B}\mathbf{u}\left(t\right),\mathbf{q}\left(t_0\right)={\mathbf{q}}_0,t\in [t_0,t_f],$(1) where $\mathbf{q}\left(t\right)\in {\mathbb{R}}^n$ is the state vector discretised over n spatial points at time t, $\mathbf{u}\left(t\right)\in {\mathbb{R}}^m$ denotes the m inputs at time t, typically related to boundary conditions or forcing terms, $t_0$ and $t_f$ are respectively the initial and final time, and ${\mathbf{q}}_0$ is the given initial condition. We refer to Equation (1(1) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{q}\left(t\right)=\mathbf{c}+\mathbf{A}\mathbf{q}\left(t\right)+\mathbf{H}\left(\mathbf{q}\right(t)\otimes \mathbf{q}(t\left)\right)+\mathbf{B}\mathbf{u}\left(t\right),\mathbf{q}\left(t_0\right)={\mathbf{q}}_0,t\in [t_0,t_f],$(1) ) as the full-order model (FOM) and note that it has been written to have a polynomial structure: $\mathbf{c}\in {\mathbb{R}}^n$ are constant terms; $\mathbf{A}\mathbf{q}\left(t\right)$ are the terms that are linear in the state $\mathbf{q}\left(t\right)$, with the discretised operator $\mathbf{A}\in {\mathbb{R}}^{n\times n}$; $\mathbf{H}\left(\mathbf{q}\right(t)\otimes \mathbf{q}(t\left)\right)$ are the terms that are quadratic in $\mathbf{q}\left(t\right)$, with $\mathbf{H}\in {\mathbb{R}}^{n\times n^2}$; and $\mathbf{B}\mathbf{u}\left(t\right)$ are the terms that are linear in the input $\mathbf{u}\left(t\right)$, with $\mathbf{B}\in {\mathbb{R}}^{n\times m}$. This polynomial structure arises in three ways: (1) it may be an attribute of the governing equations; (2) it may be exposed via variable transformations; or (3) it may be derived by introducing auxiliary variables through the process of lifting. As examples of each: (1) the incompressible Navier-Stokes equations have a quadratic form; (2) the Euler equations can be transformed to quadratic form by using pressure, velocity, and specific volume as state variables; (3) a nonlinear tubular reactor model with Arrhenius reaction terms can be written in quadratic form via the lifting transformation shown in Kramer and Willcox (2019) that introduces six auxiliary variables. Higher-order terms may be included in the formulation as well, but in this work we focus on a quadratic model due to the nature of the driving application.

A projection-based reduced-order model (ROM) of Equation (1(1) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{q}\left(t\right)=\mathbf{c}+\mathbf{A}\mathbf{q}\left(t\right)+\mathbf{H}\left(\mathbf{q}\right(t)\otimes \mathbf{q}(t\left)\right)+\mathbf{B}\mathbf{u}\left(t\right),\mathbf{q}\left(t_0\right)={\mathbf{q}}_0,t\in [t_0,t_f],$(1) ) preserves the polynomial structure (Benner et al. 2015). Approximating the high-dimensional state $\mathbf{q}\left(t\right)$ in a low-dimensional basis $\mathbf{V}\in {\mathbb{R}}^{n\times r}$, with $r\ll n$, we write $\mathbf{q}\left(t\right)\approx \mathbf{V}{\mathbf{q}}_r\left(t\right)$, where ${\mathbf{q}}_r\left(t\right)\in {\mathbb{R}}^r$ is the reduced state. Using a Galerkin projection, this yields the intrusive ROM of Equation (1(1) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{q}\left(t\right)=\mathbf{c}+\mathbf{A}\mathbf{q}\left(t\right)+\mathbf{H}\left(\mathbf{q}\right(t)\otimes \mathbf{q}(t\left)\right)+\mathbf{B}\mathbf{u}\left(t\right),\mathbf{q}\left(t_0\right)={\mathbf{q}}_0,t\in [t_0,t_f],$(1) ): $\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{q}}_r\left(t\right)={\mathbf{c}}_r+{\mathbf{A}}_r{\mathbf{q}}_r\left(t\right)+{\mathbf{H}}_r\left({\mathbf{q}}_r\right(t)\otimes {\mathbf{q}}_r(t\left)\right)+{\mathbf{B}}_r\mathbf{u}\left(t\right),{\mathbf{q}}_r\left(t_0\right)={\mathbf{V}}^{\mathrm{\top }}{\mathbf{q}}_0,t\in [t_0,t_f],$where ${\mathbf{c}}_r={\mathbf{V}}^{\mathrm{\top }}\mathbf{c}\in {\mathbb{R}}^r$, ${\mathbf{A}}_r={\mathbf{V}}^{\mathrm{\top }}\mathbf{A}\mathbf{V}\in {\mathbb{R}}^{r\times r}$, ${\mathbf{H}}_r={\mathbf{V}}^{\mathrm{\top }}\mathbf{H}(\mathbf{V}\otimes \mathbf{V})\in {\mathbb{R}}^{r\times r^2}$, and ${\mathbf{B}}_r={\mathbf{V}}^{\mathrm{\top }}\mathbf{B}\in {\mathbb{R}}^{r\times m}$ are the ROM operators corresponding to the FOM operators $\mathbf{c}$, $\mathbf{A}$, $\mathbf{H}$, and $\mathbf{B}$, respectively. The ROM is intrusive because computing these ROM operators requires access to the discretised FOM operators, which typically entails intrusive queries and/or access to source code.

The non-intrusive Operator Inference (OpInf) approach proposed by Peherstorfer and Willcox (2016) parallels the intrusive projection-based ROM setting, but learns ROMs from simulation data without direct access to the FOM operators. Recognising that the intrusive ROM has the same polynomial form as Equation (1(1) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{q}\left(t\right)=\mathbf{c}+\mathbf{A}\mathbf{q}\left(t\right)+\mathbf{H}\left(\mathbf{q}\right(t)\otimes \mathbf{q}(t\left)\right)+\mathbf{B}\mathbf{u}\left(t\right),\mathbf{q}\left(t_0\right)={\mathbf{q}}_0,t\in [t_0,t_f],$(1) ), OpInf uses a data-driven regression approach to derive a ROM of Equation (1(1) $\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{q}\left(t\right)=\mathbf{c}+\mathbf{A}\mathbf{q}\left(t\right)+\mathbf{H}\left(\mathbf{q}\right(t)\otimes \mathbf{q}(t\left)\right)+\mathbf{B}\mathbf{u}\left(t\right),\mathbf{q}\left(t_0\right)={\mathbf{q}}_0,t\in [t_0,t_f],$(1) ) as (2) $\frac{\mathrm{d}}{\mathrm{d}t}\hat{\mathbf{q}}\left(t\right)=\hat{\mathbf{c}}+\hat{\mathbf{A}}\hat{\mathbf{q}}\left(t\right)+\hat{\mathbf{H}}\left(\hat{\mathbf{q}}\right(t)\otimes \hat{\mathbf{q}}(t\left)\right)+\hat{\mathbf{B}}\mathbf{u}\left(t\right),\hat{\mathbf{q}}\left(t_0\right)={\mathbf{V}}^{\mathrm{\top }}{\mathbf{q}}_0,t\in [t_0,t_f],$(2) where $\hat{\mathbf{c}}\in {\mathbb{R}}^r$, $\hat{\mathbf{A}}\in {\mathbb{R}}^{r\times r}$, $\hat{\mathbf{H}}\in {\mathbb{R}}^{r\times r\left(r+1\right)/2}$, and $\hat{\mathbf{B}}\in {\mathbb{R}}^{r\times m}$ are determined by solving a data-driven regression problem, and $\hat{\mathbf{q}}\left(t\right)\in {\mathbb{R}}^r$ is the state of the OpInf ROM.¹

OpInf solves a regression problem to find reduced operators that yield the ROM that best matches projected snapshot data in a minimum-residual sense. Mathematically, OpInf solves the least-squares problem (3) $\underset{\hat{\mathbf{c}},\hat{\mathbf{A}},\hat{\mathbf{H}},\hat{\mathbf{B}}}{min}\sum _{j=0}^{k-1}{∥\hat{\mathbf{c}}+\hat{\mathbf{A}}{\hat{\mathbf{q}}}_j+\hat{\mathbf{H}}({\hat{\mathbf{q}}}_j\otimes {\hat{\mathbf{q}}}_j)+\hat{\mathbf{B}}{\mathbf{u}}_j-{\dot{\hat{\mathbf{q}}}}_j∥}_2^2,$(3) where $\{{\hat{\mathbf{q}}}_j{\}}_{j=0}^{k-1}$ is the dataset used to drive the learning with ${\hat{\mathbf{q}}}_j$ denoting a reduced-state snapshot at timestep j, $\{{\dot{\hat{\mathbf{q}}}}_j{\}}_{j=0}^{k-1}$ are the associated time derivatives, and $\{{\mathbf{u}}_j{\}}_{j=0}^{k-1}$ is the collection of inputs corresponding to the data with ${\mathbf{u}}_j\equiv \mathbf{u}\left(t_j\right)$. To generate the dataset, we employ the following steps: (1) Collect a set of k high-fidelity state snapshots $\{{\mathbf{z}}_j{\}}_{j=0}^{k-1}\subset {\mathbb{R}}^{\tilde{n}}$ by solving the original high-fidelity model at times $\{t_j{\}}_{j=0}^{k-1}$ with inputs $\{{\mathbf{u}}_j{\}}_{j=0}^{k-1}$. (2) Apply a variable transformation ${\mathbf{q}}_j=\mathcal{T}\left({\mathbf{z}}_j\right)$ to obtain snapshots of the transformed variables. Here $\mathcal{T}:{\mathbb{R}}^{\tilde{n}}\to {\mathbb{R}}^n$ is the map representing a reversible transformation (e.g. from density to specific volume) or a lifting transformation (Qian et al. 2020). (3) Compute the proper orthogonal decomposition (POD) basis $\mathbf{V}$ of the transformed snapshots.² (4) Project the transformed snapshots onto the POD subspace as ${\hat{\mathbf{q}}}_j={\mathbf{V}}^{\mathrm{\top }}{\mathbf{q}}_j$. (5) Estimate projected time derivative information $\{\dot{\hat{\mathbf{q}}}{\}}_{j=0}^{k-1}$. The training period $[t_0,t_{k-1}]$ for which we have data is a subset of the full time domain of interest $[t_0,t_f]$; results from the ROM over $[t_k,t_f]$ will be entirely predictive.

Equation (3(3) $\underset{\hat{\mathbf{c}},\hat{\mathbf{A}},\hat{\mathbf{H}},\hat{\mathbf{B}}}{min}\sum _{j=0}^{k-1}{∥\hat{\mathbf{c}}+\hat{\mathbf{A}}{\hat{\mathbf{q}}}_j+\hat{\mathbf{H}}({\hat{\mathbf{q}}}_j\otimes {\hat{\mathbf{q}}}_j)+\hat{\mathbf{B}}{\mathbf{u}}_j-{\dot{\hat{\mathbf{q}}}}_j∥}_2^2,$(3) ) can also be written in matrix form as (4) $\underset{\mathbf{O}}{min}{∥\mathbf{D}{\mathbf{O}}^{\mathrm{\top }}-{\mathbf{R}}^{\mathrm{\top }}∥}_F^2,$(4) where $\begin{array}{rl}\mathbf{O}& =\left[\hat{\mathbf{c}}\hat{\mathbf{A}}\hat{\mathbf{H}}\hat{\mathbf{B}}\right]\in {\mathbb{R}}^{r\times d(r,m)},\left(\mathrm{u}\mathrm{n}\mathrm{k}\mathrm{n}\mathrm{o}\mathrm{w}\mathrm{n}\mathrm{o}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{r}\mathrm{s}\right)\\ \mathbf{D}& =\left[{\mathbf{1}}_k{\hat{\mathbf{Q}}}^{\mathrm{\top }}(\hat{\mathbf{Q}}\otimes \hat{\mathbf{Q}}{)}^{\mathrm{\top }}{\mathbf{U}}^{\mathrm{\top }}\right]\in {\mathbb{R}}^{k\times d(r,m)},\left(\mathrm{k}\mathrm{n}\mathrm{o}\mathrm{w}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{t}\mathrm{a}\right)\\ \hat{\mathbf{Q}}& =\left[{\hat{\mathbf{q}}}_0{\hat{\mathbf{q}}}_1\cdots {\hat{\mathbf{q}}}_{k-1}\right]\in {\mathbb{R}}^{r\times k},\left(\mathrm{s}\mathrm{n}\mathrm{a}\mathrm{p}\mathrm{s}\mathrm{h}\mathrm{o}\mathrm{t}\mathrm{s}\right)\\ \mathbf{R}& =\left[{\dot{\hat{\mathbf{q}}}}_0{\dot{\hat{\mathbf{q}}}}_1\cdots {\dot{\hat{\mathbf{q}}}}_{k-1}\right]\in {\mathbb{R}}^{r\times k},\left(\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{s}\right)\\ \mathbf{U}& =\left[{\mathbf{u}}_0{\mathbf{u}}_1\cdots {\mathbf{u}}_{k-1}\right]\in {\mathbb{R}}^{m\times k},\left(\mathrm{i}\mathrm{n}\mathrm{p}\mathrm{u}\mathrm{t}\mathrm{s}\right)\end{array}$and where $d(r,m)=1+r+r\left(r+1\right)/2+m$ and ${\mathbf{1}}_k\in {\mathbb{R}}^k$ is the length-k column vector with all entries set to unity. The OpInf least-squares problem Equation (3(3) $\underset{\hat{\mathbf{c}},\hat{\mathbf{A}},\hat{\mathbf{H}},\hat{\mathbf{B}}}{min}\sum _{j=0}^{k-1}{∥\hat{\mathbf{c}}+\hat{\mathbf{A}}{\hat{\mathbf{q}}}_j+\hat{\mathbf{H}}({\hat{\mathbf{q}}}_j\otimes {\hat{\mathbf{q}}}_j)+\hat{\mathbf{B}}{\mathbf{u}}_j-{\dot{\hat{\mathbf{q}}}}_j∥}_2^2,$(3) ) is therefore linear in the coefficients of the unknown ROM operators $\hat{\mathbf{c}}$, $\hat{\mathbf{A}}$, $\hat{\mathbf{H}}$, and $\hat{\mathbf{B}}$.

The OpInf approach permits us to compute the ROM operators $\hat{\mathbf{c}}$, $\hat{\mathbf{A}}$, $\hat{\mathbf{H}}$, and $\hat{\mathbf{B}}$ without explicit access to the original high-dimensional operators $\mathbf{c}$, $\mathbf{A}$, $\mathbf{H}$, and $\mathbf{B}$. This point is key since we apply variable transformations only to the snapshot data, not to the operators or the underlying model. Thus, even in a setting where deriving a classical intrusive ROM might be possible, the OpInf approach enables us to work with variables other than those used for the original high-fidelity discretisation. We will see the importance of this later for a reacting flow application. We also note that under some conditions, OpInf recovers the intrusive POD ROM (Peherstorfer and Willcox 2016; Peherstorfer 2020).

Regularisation

The problem in Equation (4(4) $\underset{\mathbf{O}}{min}{∥\mathbf{D}{\mathbf{O}}^{\mathrm{\top }}-{\mathbf{R}}^{\mathrm{\top }}∥}_F^2,$(4) ) is generally overdetermined (i.e. $k>d(r,m)$), but is also noisy due to errors in the numerically estimated time derivatives $\mathbf{R}$, model mis-specification (e.g. if the system is not truly quadratic), and truncated POD modes that leave some system dynamics unresolved. The ROMs resulting from Equation (4(4) $\underset{\mathbf{O}}{min}{∥\mathbf{D}{\mathbf{O}}^{\mathrm{\top }}-{\mathbf{R}}^{\mathrm{\top }}∥}_F^2,$(4) ) can thus suffer from overfitting the operators to the data and therefore exhibit poor predictive performance over the time domain of interest $[t_0,t_f]$.

To avoid overfitting, we introduce a Tikhonov regularisation (Tikhonov and Arsenin 1977) to Equation (4(4) $\underset{\mathbf{O}}{min}{∥\mathbf{D}{\mathbf{O}}^{\mathrm{\top }}-{\mathbf{R}}^{\mathrm{\top }}∥}_F^2,$(4) ), which then becomes (5) $\underset{\mathbf{O}}{min}{∥\mathbf{D}{\mathbf{O}}^{\mathrm{\top }}-{\mathbf{R}}^{\mathrm{\top }}∥}_F^2+{∥\mathbf{\Gamma }{\mathbf{O}}^{\mathrm{\top }}∥}_F^2,$(5) where $\mathbf{\Gamma }\in {\mathbb{R}}^{d(r,m)\times d(r,m)}$ is a full-rank regulariser. The minimiser of Equation (5(5) $\underset{\mathbf{O}}{min}{∥\mathbf{D}{\mathbf{O}}^{\mathrm{\top }}-{\mathbf{R}}^{\mathrm{\top }}∥}_F^2+{∥\mathbf{\Gamma }{\mathbf{O}}^{\mathrm{\top }}∥}_F^2,$(5) ) satisfies the modified normal equations (6) $\left({\mathbf{D}}^{\mathrm{\top }}\mathbf{D}+{\mathbf{\Gamma }}^{\mathrm{\top }}\mathbf{\Gamma }\right){\mathbf{O}}^{\mathrm{\top }}={\mathbf{D}}^{\mathrm{\top }}{\mathbf{R}}^{\mathrm{\top }},$(6) which admit a unique solution since ${\mathbf{D}}^{\mathrm{\top }}\mathbf{D}+{\mathbf{\Gamma }}^{\mathrm{\top }}\mathbf{\Gamma }$ is symmetric positive definite.

An $L_2$ regulariser $\mathbf{\Gamma }=\lambda \mathbf{I}$, $\lambda >0$ and $\mathbf{I}$ the identity matrix, penalises each entry of the inferred ROM operators $\hat{\mathbf{c}}$, $\hat{\mathbf{A}}$, $\hat{\mathbf{H}}$, and $\hat{\mathbf{B}}$, thereby driving the ROM toward the globally stable system $\frac{\mathrm{d}}{\mathrm{d}t}\hat{\mathbf{q}}\left(t\right)=\mathbf{0}$. Since the entries of the quadratic operator $\hat{\mathbf{H}}$ have a different scaling than entries of the other operators, we construct a diagonal regulariser $\mathbf{\Lambda }({\lambda }_1,{\lambda }_2)$, with ${\lambda }_1,{\lambda }_2>0$, such that the operator entries are penalised by ${\lambda }_1$, except for the entries of $\hat{\mathbf{H}}$, which are penalised by ${\lambda }_2$. That is, with $\mathbf{\Gamma }=\mathbf{\Lambda }({\lambda }_1,{\lambda }_2)$, Equation (5(5) $\underset{\mathbf{O}}{min}{∥\mathbf{D}{\mathbf{O}}^{\mathrm{\top }}-{\mathbf{R}}^{\mathrm{\top }}∥}_F^2+{∥\mathbf{\Gamma }{\mathbf{O}}^{\mathrm{\top }}∥}_F^2,$(5) ) can be expressed as $\underset{\hat{\mathbf{c}},\hat{\mathbf{A}},\hat{\mathbf{H}},\hat{\mathbf{B}}}{min}\left[{∥\mathbf{D}{\left[\hat{\mathbf{c}}\hat{\mathbf{A}}\hat{\mathbf{H}}\hat{\mathbf{B}}\right]}^{\mathrm{\top }}-{\mathbf{R}}^{\mathrm{\top }}∥}_F^2+{\lambda }_1\left(\parallel \hat{\mathbf{c}}{\parallel }_2^2+\parallel \hat{\mathbf{A}}{\parallel }_F^2+\parallel \hat{\mathbf{B}}{\parallel }_F^2\right)+{\lambda }_2\parallel \hat{\mathbf{H}}{\parallel }_F^2\right].$The scalar hyperparameters ${\lambda }_1$ and ${\lambda }_2$, which balance the minimisation between the data fit and the regularisation, must be chosen with care. The ideal regulariser produces a ROM that minimises some error metric over the full time domain $[t_0,t_f]$; however, since data are only available for the smaller training domain $[t_0,t_{k-1}]$, we choose ${\lambda }_1$ and ${\lambda }_2$ so that the resulting ROM minimises error over $[t_0,t_{k-1}]$ while maintaining a bound on the integrated POD coefficients over $[t_0,t_f]$. That is, we require the state $\hat{\mathbf{q}}\left(t\right)=\left[{\hat{q}}_1\right(t\left){\hat{q}}_2\right(t)\cdots {\hat{q}}_r(t){]}^{\mathrm{\top }}$ produced by integrating Equation (2(2) $\frac{\mathrm{d}}{\mathrm{d}t}\hat{\mathbf{q}}\left(t\right)=\hat{\mathbf{c}}+\hat{\mathbf{A}}\hat{\mathbf{q}}\left(t\right)+\hat{\mathbf{H}}\left(\hat{\mathbf{q}}\right(t)\otimes \hat{\mathbf{q}}(t\left)\right)+\hat{\mathbf{B}}\mathbf{u}\left(t\right),\hat{\mathbf{q}}\left(t_0\right)={\mathbf{V}}^{\mathrm{\top }}{\mathbf{q}}_0,t\in [t_0,t_f],$(2) ) to satisfy $\left|{\hat{q}}_i\right(t\left)\right|\le B$, $i=1,\dots ,r$ and $t\in [t_0,t_f]$, for some B>0. This in turn ensures a bound on the magnitude of the entries of the high-dimensional state $\mathbf{q}\left(t\right)=\mathbf{V}\hat{\mathbf{q}}\left(t\right)$: (7) $\left|q_i\right(t\left)\right|=\left|\sum _{j=1}^r{V}_{ij}{\hat{q}}_j\left(t\right)\right|\le \sum _{j=1}^r\left|V_{ij}\right|\left|{\hat{q}}_j\left(t\right)\right|\le B\sum _{j=1}^r\left|V_{ij}\right|,i=1,2,\dots ,n,$(7) where $V_{ij}$ is the ith element of the jth POD basis vector. Note that the bound B may be chosen with the intent of imposing a particular bound on the $\left|q_i\right|$ since the sums $\sum _{j=1}^r\left|V_{ij}\right|$ can be precomputed. For example, our regularisation strategy provides a computationally efficient way to impose the temperature limiters proposed by Huang et al. (2019).

Algorithm 1 details our regularised OpInf procedure, in which we choose B as a multiple of the maximum absolute entry of the projected training data $\hat{\mathbf{Q}}$. This particular strategy for selecting ${\lambda }_1$ and ${\lambda }_2$ could be replaced with a cross-validation or resampling grid search technique, but our experiments in this vein did not produce robust results. The training error $\parallel \hat{\mathbf{Q}}-{\tilde{\mathbf{Q}}}_{:,:k}\parallel$ in step 13 may compare $\hat{\mathbf{Q}}$ and ${\tilde{\mathbf{Q}}}_{:,:k}$ directly in the reduced space (e.g. with a matrix norm or an $L^p\left(\right[t_0,t_{k-1}\left]\right)$ norm), or it may be replaced with a more targeted comparison of some quantity of interest.

Table

Algorithm 1: Operator Inference with regularisation selection
1:	procedure RegOpInf(snapshots $\mathbf{Z}\in {\mathbb{R}}^{\tilde{n}\times k}$, inputs $\mathbf{U}\in {\mathbb{R}}^{m\times k}$, final time $t_f>t_0$, reversible map $\mathcal{T}:{\mathbb{R}}^{\tilde{n}}\to {\mathbb{R}}^n$, reduced dimension $r\in \mathbb{N}$, bound margin $\tau \ge 1$ )
2:	$\mathbf{Q}\leftarrow \mathcal{T}\left(\mathbf{Z}\right)$	$▹$ Map native variables to learning variables (columnwise).
3:	$\mathbf{V}\leftarrow$pod($\mathbf{Q},r$)	$▹$ Compute a rank-$r$ POD basis from transformed snapshots.
4:	$\hat{\mathbf{Q}}\leftarrow V^{\mathrm{\top }}\mathbf{Q}$	$▹$ Project snapshots to the $r$-dimensional POD subspace.
5:	$\mathbf{R}\leftarrow {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}\hat{\mathbf{Q}}}$	$▹$ Estimate time derivatives of the projected snapshots.
6:	$B\leftarrow \tau \underset{i,j}{max}\left|{\hat{\mathbf{Q}}}_{ij}\right|$	$▹$ Select a bound to require for integrated POD modes.
7:	procedure TrainError(${\lambda }_1,{\lambda }_2$)
8:	$\hat{\mathbf{c}},\hat{\mathbf{A}},\hat{\mathbf{H}},\hat{\mathbf{B}}\leftarrow$ solve Eq. (5) with regulariser $\mathbf{\Gamma }=\mathbf{\Lambda }\left({\lambda }_1,{\lambda }_2\right)$
9:	$\tilde{\mathbf{Q}}\leftarrow$ integrate Eq. (2) with $\hat{\mathbf{c}},\hat{\mathbf{A}},\hat{\mathbf{H}},\hat{\mathbf{B}}$ from ${\hat{q}}_0={\hat{\mathbf{Q}}}_{:,0}$ over $\left[t_0,t_f\right]$
10:	if $\underset{i,j}{max}\left|{\tilde{\mathbf{Q}}}_{ij}\right|>B$ then
11:	return $\mathrm{\infty }$	$▹$ Disqualify ROMs that violate the bound.
12:	else
13:	return $\parallel \hat{\mathbf{Q}}-{\tilde{\mathbf{Q}}}_{:,:k}\parallel$	$▹$ Compute ROM error over $\left[t_0,t_{k-1}\right]$.
14:	${\lambda }_1^{\ast },{\lambda }_2^{\ast }\leftarrow \mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}$ TrainError$\left({\lambda }_1,{\lambda }_2\right)$
15:	$\hat{\mathbf{c}},\hat{\mathbf{A}},\hat{\mathbf{H}},\hat{\mathbf{B}}\leftarrow$ solve Eq. (5) with optimal regulariser $\mathbf{\Gamma }=\mathbf{\Lambda }\left({\lambda }_1^{\ast },{\lambda }_2^{\ast }\right)$
16:	return $\hat{\mathbf{c}},\hat{\mathbf{A}},\hat{\mathbf{H}},\hat{\mathbf{B}}$

The regularisation approach of Algorithm 1 may be viewed as a stabilisation method since it selects a ROM with reasonable behaviour over a given time domain. It should be noted, however, that this method does not modify an existing ROM to achieve stability, different from other stabilisation methods such as eigenvalue reassignment (Kalashnikova et al. 2014; Rezaian and Wei 2020). The optimisation is driven by a penalisation but has no built-in constraints, which is a major advantage in terms of the computational cost, but the resulting ROMs are not guaranteed to preserve properties such as energy conservation. Adding constraints to Operator Inference to target conservation, similar to the work in Carlberg et al. (2018), is a subject for possible future work.

Scalable implementation

The steps of Algorithm 1 are highly modular and amenable to large-scale problems. The variable transformation $\mathbf{Q}=\mathcal{T}\left(\mathbf{Z}\right)$ of step 2 consists of $O\left(nk\right)$ elementary computations on the original snapshot matrix $\mathbf{Z}\in {\mathbb{R}}^{\tilde{n}\times k}$. To compute the rank-r POD basis $\mathbf{V}\in {\mathbb{R}}^{n\times r}$ in step 3, we use a randomised singular-value decomposition algorithm requiring $O(rnk+r^2(n+k\left)\right)$ operations (Halko et al. 2011); since $r\ll k,n$, the leading order behaviour is $O\left(rnk\right)$. The projection $\hat{\mathbf{Q}}={\mathbf{V}}^{\mathrm{\top }}\mathbf{Q}$ in step 4 costs about 2nk operations. Note that $\hat{\mathbf{Q}}\in {\mathbb{R}}^{r\times k}$ is small compared to $\mathbf{Q}\in {\mathbb{R}}^{n\times k}$, as typically $r\sim {10}^0$–${10}^3$ and $n\sim {10}^4$–${10}^9$. The time derivatives $\mathbf{R}\in {\mathbb{R}}^{r\times k}$ in step 5 may be provided by the full-order solver or estimated, e.g. with finite differences of the columns of $\hat{\mathbf{Q}}$. In the latter case, the cost is $O\left(rk\right)$. Step 6 is a simple $O\left(rk\right)$ selection. The computational cost of steps 2–6 is therefore $O\left(rnk\right)$.

The subroutine TrainError$:{\mathbb{R}}^2\to \mathbb{R}$ comprising steps 7–13 must be evaluated for several choices of ${\lambda }_1$ and ${\lambda }_2$. In step 8, we solve the minimisation problem Equation (5(5) $\underset{\mathbf{O}}{min}{∥\mathbf{D}{\mathbf{O}}^{\mathrm{\top }}-{\mathbf{R}}^{\mathrm{\top }}∥}_F^2+{∥\mathbf{\Gamma }{\mathbf{O}}^{\mathrm{\top }}∥}_F^2,$(5) ) via the linear system Equation (6(6) $\left({\mathbf{D}}^{\mathrm{\top }}\mathbf{D}+{\mathbf{\Gamma }}^{\mathrm{\top }}\mathbf{\Gamma }\right){\mathbf{O}}^{\mathrm{\top }}={\mathbf{D}}^{\mathrm{\top }}{\mathbf{R}}^{\mathrm{\top }},$(6) ). Forming ${\mathbf{D}}^{\mathrm{\top }}\mathbf{D}\in {\mathbb{R}}^{d(r,m)\times d(r,m)}$ and ${\mathbf{D}}^{\mathrm{\top }}{\mathbf{R}}^{\mathrm{\top }}\in {\mathbb{R}}^{d(r,m)\times r}$ costs $O\left(d\right(r,m{)}^{2}k+rd(r,m)k)=O(r^{4}k)$ operations, but these can be precomputed once and reused in each subroutine evaluation; for specific ${\lambda }_1$ and ${\lambda }_2$, we form $\mathbf{\Gamma }\in {\mathbb{R}}^{d(r,m)\times d(r,m)}$ and compute ${\mathbf{D}}^{\mathrm{\top }}\mathbf{D}+{\mathbf{\Gamma }}^{\mathrm{\top }}\mathbf{\Gamma }\in {\mathbb{R}}^{d(r,m)\times d(r,m)}$ with only $2d(r,m)$ additional operations since $\mathbf{\Gamma }$ is diagonal. Solving Equation (6(6) $\left({\mathbf{D}}^{\mathrm{\top }}\mathbf{D}+{\mathbf{\Gamma }}^{\mathrm{\top }}\mathbf{\Gamma }\right){\mathbf{O}}^{\mathrm{\top }}={\mathbf{D}}^{\mathrm{\top }}{\mathbf{R}}^{\mathrm{\top }},$(6) ) via the Cholesky decomposition has a leading order cost of $\frac{1}{3}\left(d\right(r,m{)}^3)=\frac{1}{24}{r}^6$ operations (Demmel 1997), but this estimate can be improved upon with iterative methods for larger r if desired. The ROM integration in step 9 can be carried out with any time-stepping scheme; for explicit methods, evaluating the ROM at a single point, i.e. computing the right-hand side of Equation (2(2) $\frac{\mathrm{d}}{\mathrm{d}t}\hat{\mathbf{q}}\left(t\right)=\hat{\mathbf{c}}+\hat{\mathbf{A}}\hat{\mathbf{q}}\left(t\right)+\hat{\mathbf{H}}\left(\hat{\mathbf{q}}\right(t)\otimes \hat{\mathbf{q}}(t\left)\right)+\hat{\mathbf{B}}\mathbf{u}\left(t\right),\hat{\mathbf{q}}\left(t_0\right)={\mathbf{V}}^{\mathrm{\top }}{\mathbf{q}}_0,t\in [t_0,t_f],$(2) ), costs $O\left(r^3\right)$ operations. The total cost each evaluation of TrainError$:{\mathbb{R}}^2\to \mathbb{R}$ is therefore dominated by the cost of solving Equation (6(6) $\left({\mathbf{D}}^{\mathrm{\top }}\mathbf{D}+{\mathbf{\Gamma }}^{\mathrm{\top }}\mathbf{\Gamma }\right){\mathbf{O}}^{\mathrm{\top }}={\mathbf{D}}^{\mathrm{\top }}{\mathbf{R}}^{\mathrm{\top }},$(6) ), which is independent of n and k.

Finally, the minimisation in step 14 is carried out with a derivative-free search method, which enables fewer total evaluations of the subroutine than a fine grid search. However, a coarse grid search is useful for identifying appropriate initial guesses for ${\lambda }_1$ and ${\lambda }_2$.

Results

This section applies regularised OpInf to a single-injector combustion problem, studied previously by Swischuk et al. (2020), on the two-dimensional computational domain shown in Figure 1. We first describe the governing dynamics, a set of high-fidelity data obtained from a CFD code, and the variable transformations used to produce training data for learning reduced models with Algorithm 1. The resulting OpInf ROM performance is then analysed and compared to a state-of-the-art intrusive model reduction method.

Figure 1. The computational domain for the single-injector combustion problem. On the left, monitor locations for numerical results. On the right, a typical temperature field demonstrating the complexity and nonlinear nature of the problem.

Problem setup

The combustion dynamics for this problem are governed by conservation laws (8) $\frac{\mathrm{\partial }{\overrightarrow{q}}_{\mathrm{c}}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\overrightarrow{K}-{\overrightarrow{K}}_{\mathrm{v}})=\overrightarrow{S},$(8) where ${\overrightarrow{q}}_{\mathrm{c}}=[\rho \rho v_x\rho v_y\rho e\rho Y_1\rho Y_2\rho Y_3\rho Y_4{]}^{\mathrm{\top }}$ are the conservative variables, $\overrightarrow{K}$ are the inviscid flux terms, ${\overrightarrow{K}}_{\mathrm{v}}$ are the viscous flux terms, and $\overrightarrow{S}$ are the source terms. Here ρ is the density $\left[\frac{\mathrm{k}\mathrm{g}}{{\mathrm{m}}^3}\right]$, $v_x$ and $v_y$ are the x and y velocity $\left[\frac{\mathrm{m}}{\mathrm{s}}\right]$, e is the total energy $\left[\frac{\mathrm{J}}{{\mathrm{m}}^3}\right]$, and $Y_{\ell }$ is the ℓth species mass fraction with $\ell =1,2,3,4$. The chemical species are CH${}_4$, O${}_2$, H${}_2$O, and CO${}_2$, which follow a global one-step chemical reaction ${\mathrm{C}\mathrm{H}}_4+2{\mathrm{O}}_2\to {\mathrm{C}\mathrm{O}}_2+2{\mathrm{H}}_2\mathrm{O}$ (Westbrook and Dryer 1981). See Harvazinski et al. (2015) for more details on the governing equations.

At the downstream end of the combustor, we impose a non-reflecting boundary condition while maintaining the chamber pressure via (9) $p_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}}\left(t\right)=p_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k},\mathrm{r}\mathrm{e}\mathrm{f}}\left(1+0.1\sin \left(2\pi ft\right)\right)$(9) where $p_{\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k},\mathrm{r}\mathrm{e}\mathrm{f}}={10}^6\mathrm{P}\mathrm{a}$ and $f=\mathrm{5,000}\mathrm{H}\mathrm{z}$. The top and bottom wall boundary conditions are no-slip conditions, and for the upstream boundary we impose a constant mass flow at the inlets. Swischuk et al. (2020) show that if the governing equations (8(8) $\frac{\mathrm{\partial }{\overrightarrow{q}}_{\mathrm{c}}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\overrightarrow{K}-{\overrightarrow{K}}_{\mathrm{v}})=\overrightarrow{S},$(8) ) are transformed to be written in the specific volume variables, many of the terms take a quadratic form. Following that idea, we choose as learning variables the transformed and augmented state $\overrightarrow{q}=[p{v}_x{v}_{y}T\xi c_1{c}_2{c}_3{c}_4{]}^{\mathrm{\top }}$ where p is the pressure $\left[\mathrm{P}\mathrm{a}\right]$, T is the temperature $\left[\mathrm{K}\right]$, $\xi =1/\rho$ is the specific volume $[\frac{{\mathrm{m}}^3}{\mathrm{k}\mathrm{g}}]$, and $c_1,\dots ,c_4$ are the species molar concentrations $[\frac{\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}}{{\mathrm{m}}^3}]$ given by $c_{\ell }=\rho Y_{\ell }/M_{\ell }$ with $M_{\ell }$ the molar mass of the ℓth species $[\frac{\mathrm{g}}{\mathrm{m}\mathrm{o}\mathrm{l}}]$. As shown in Swischuk et al. (2020), the equations for $v_x$, $v_y$, and ξ are exactly quadratic in $\overrightarrow{q}$, while the remaining equations are quadratic with some non-polynomial terms in $\overrightarrow{q}$. Note that, differently from Swischuk et al. (2020), $\overrightarrow{q}$ here is chosen to contain specific volume, pressure, and temperature, even though only two of the three quantities are needed to fully define the high-fidelity model (and the equation of state then defines the third). We augment the learning variables in this way because doing so exposes the quadratic form while also directly targeting the variables that are of primary interest for assessing ROM performance. In particular, the resulting ROMs provide more accurate predictions of temperature when temperature is included explicitly as a learning variable. We can do this since the transformations are applied only to the snapshot data, not to the CFD model itself. Indeed, constructing the full-order spatial operators in these transformed coordinates would be both impractical and inexact. This flexibility to learn from transformed snapshots instead of transformed full-order operators is a major advantage of the non-intrusive OpInf approach in comparison to traditional intrusive projection-based model reduction methods.

To generate high-fidelity training data, we use the finite-volume based General Equation and Mesh Solver (GEMS) (Harvazinski et al. 2015) to solve for the variables $\left[p{v}_x{v}_{y}T{Y}_1{Y}_2{Y}_3{Y}_4\right]$ over $n_x=\mathrm{38,523}$ cells, resulting in snapshots with $8n_x=\mathrm{308,184}$ entries each. The snapshots are computed for 60,000 time steps beyond the initial condition with a temporal discretisation of $\delta t={10}^{-7}\mathrm{s}$, from $t_0=0.015\mathrm{s}$ to $t_f=0.021\mathrm{s}$. The computational cost of computing this dataset is approximately 1200 CPU hours on two computing nodes, each of which contains two Haswell CPUs at 2.60 GHz and 20 cores per node.

Scaling is an essential aspect of successful model reduction and is particularly critical for this problem due to the wide range of scales across variables. After transforming the GEMS snapshot data to the learning variables $\overrightarrow{q}$, the species molar concentrations are scaled to $[0,1]$, and all other variables are scaled to $[-1,1]$. This scaling ensures that null velocities and null molar concentrations are preserved. For example, some upstream regions of the injector have zero methane concentration at all times. By construction, the POD basis vectors and thus the ROM predictions will preserve those zero concentration values.

We implement Algorithm 1 via the rom_operator_inference Python package,³ which is built on NumPy, SciPy, and scikit-learn (Walt et al. 2011; Pedregosa et al. 2011; Virtanen et al. 2020). The time derivatives in step 5 of Algorithm 1 are estimated with fourth-order finite differences, and the least-squares problem in step 8 are solved by applying the LAPACK routine POSV to Equation (6(6) $\left({\mathbf{D}}^{\mathrm{\top }}\mathbf{D}+{\mathbf{\Gamma }}^{\mathrm{\top }}\mathbf{\Gamma }\right){\mathbf{O}}^{\mathrm{\top }}={\mathbf{D}}^{\mathrm{\top }}{\mathbf{R}}^{\mathrm{\top }},$(6) ). The learned ROMs are integrated in step 9 with the explicit, adaptive, fourth-order Runge-Kutta scheme RK45, and the error evaluation of step 13 uses the $L^2\left(\right[t_0,t_{k-1}\left]\right)$ norm in the reduced space. To minimise the function TrainError$:{\mathbb{R}}^2\to \mathbb{R}$ in step 14, we find an initial estimate of the minimiser via a coarse grid search, then refine the result with a Nelder-Mead search method (Nelder and Mead 1965). The code and details are publicly available at https://github.com/Willcox-Research-Group/ROM-OpInf-Combustion-2D.

Sensitivity to training data

We study the sensitivity of our approach to the training data by varying the number of snapshots used to compute the POD basis and learn the OpInf ROM. Specifically, we consider the three cases where we use the first $k=\mathrm{10,000}$, $k=\mathrm{20,000}$, and $k=\mathrm{30,000}$ snapshots from GEMS as training data sets. In each case we compute the POD basis and, to select an appropriate reduced dimension r, the cumulative energy based on the POD singular values: ${\mathcal{E}}_r=\left(\sum _{j=1}^r{\sigma }_j^2\right)/\left(\sum _{j=1}^k{\sigma }_j^2\right)$, where $\{{\sigma }_j{\}}_{j=1}^k$ are the singular values of the learning variable snapshot matrix $\mathbf{Q}$ (see Figure 2). Specifically, we choose the minimal integer r such that ${\mathcal{E}}_r$ is greater than a fixed energy threshold. Table 1 shows that r is increasing linearly with the number of snapshots, indicating that the basis is not being saturated as additional information is incorporated. This is an indication of the challenging nature of this application, due to the rich and complex dynamics. Table 1 also shows the column dimension $d(r,m)$ of the data matrix $\mathbf{D}\in {\mathbb{R}}^{k\times d(r,m)}$ in the Operator Inference problem, which grows quadratically with r; choosing r so that $k\gg d(r,m)$ helps ensure that $\mathbf{D}$ has full rank.

Figure 2. The POD singular values for varying size of the snapshot training set.

Table 1. The basis size r required to exceed a given cumulative energy level ${\mathcal{E}}_r$ increases linearly with the number of snapshots k in the training set; the dimension $d=d(r,m)$ increases quadratically with r.

	Cumulative energy ${\mathcal{E}}_r$
	0.985	0.990	0.995
$k=\mathrm{10,000}$	r = 22, d = 277	r = 27, d = 407	r = 36, d = 704
$k=\mathrm{20,000}$	r = 43, d = 991	r = 53, d = 1486	r = 72, d = 2701
$k=\mathrm{30,000}$	r = 66, d = 2279	r = 82, d = 3487	r = 110, d = 6217

Figure 3 plots the GEMS and OpInf ROM results for pressure and x-velocity predictions over time at two of the monitor locations in Figure 1.⁴ While it can be misleading to assess accuracy based on predictions at a single spatial point, these plots reveal several representative insights. First, each OpInf ROM faithfully reconstructs the training data but has some discrepancies in the prediction regime. Second, the pressure and x-velocity frequencies are well captured throughout the time domain, but the amplitudes are sometimes less accurate in the prediction regime. The effects of the 5000 Hz downstream pressure forcing are clearly visible in the pressure. Third, we see the importance of the training data–as the amount of training data increases, the ROM predictions change significantly and generally (but not always) improve. This is yet another indication of the complexity of the dynamics we are aiming to approximate.

Figure 3. Traces of pressure and x-velocity through time at monitor locations 1 and 3 of Figure 1, respectively, for $k=\mathrm{10,000}$ (top row), $k=\mathrm{20,000}$ (middle row), and $k=\mathrm{30,000}$ (bottom row) training snapshots. The vertical black lines separate the training and prediction periods. For each choice of k, the error over the training domain is low, but increasing k has a significant impact on the behaviour in the testing domain. Here r is chosen in each case so that ${\mathcal{E}}_r>0.985$.

The effectiveness of the low-dimensional basis differs among the state variables. Let $\mathbf{z}\left(t\right)\in {\mathbb{R}}^{8n_x}$ denote the GEMS simulation data at time t with components $\left[s_1\right(t)\cdots s_{n_x}(t){]}^{\mathrm{\top }}\in {\mathbb{R}}^{n_x}$ corresponding to a single state variable (e.g. pressure or temperature), and let $\mathcal{T}:{\mathbb{R}}^{8n_x}\to {\mathbb{R}}^{9n_x}$ be the preprocessing variable transformation and scaling with reverse operator ${\mathcal{T}}^{\ast }:{\mathbb{R}}^{9n_x}\to {\mathbb{R}}^{8n_x}$ (i.e. ${\mathcal{T}}^{\ast }\left(\mathcal{T}\right(\mathbf{z}\left(t\right)\left)\right)=\mathbf{z}\left(t\right)$ for all t). Given the basis $\mathbf{V}\in {\mathbb{R}}^{9n_x\times r}$ for the scaled learning variables, the projection of the GEMS data onto the basis is ${\mathcal{T}}^{\ast }\left(\mathbf{V}{\mathbf{V}}^{\mathrm{\top }}\mathcal{T}\right(\mathbf{z}\left(t\right)\left)\right).$ Let $\left[s_1^{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}\right(t)\cdots s_{n_x}^{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}(t){]}^{\mathrm{\top }}\in {\mathbb{R}}^{n_x}$ denote the components of the projection corresponding to the state variable of interest. Next, letting $\tilde{\mathbf{q}}\left(t\right)$ denote the ROM state (the result of integrating Equation (2(2) $\frac{\mathrm{d}}{\mathrm{d}t}\hat{\mathbf{q}}\left(t\right)=\hat{\mathbf{c}}+\hat{\mathbf{A}}\hat{\mathbf{q}}\left(t\right)+\hat{\mathbf{H}}\left(\hat{\mathbf{q}}\right(t)\otimes \hat{\mathbf{q}}(t\left)\right)+\hat{\mathbf{B}}\mathbf{u}\left(t\right),\hat{\mathbf{q}}\left(t_0\right)={\mathbf{V}}^{\mathrm{\top }}{\mathbf{q}}_0,t\in [t_0,t_f],$(2) )), the reconstructed ROM prediction in the original state variables is ${\mathcal{T}}^{\ast }\left(\mathbf{V}\tilde{\mathbf{q}}\right(t\left)\right).$ Let $\left[s_1^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}\right(t)\cdots s_{n_x}^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}(t){]}^{\mathrm{\top }}\in {\mathbb{R}}^{n_x}$ be the components of this reconstruction corresponding to the state variable of interest. We define the spatially averaged relative errors for the projection and the ROM predictions as $s_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}\mathrm{e}\mathrm{r}\mathrm{r}}\left(t\right)=\frac{1}{n_x}\sum _{i=1}^{n_x}\frac{\left|s_i^{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{j}}\right(t)-s_i(t\left)\right|}{\left|s_i\right(t\left)\right|},s_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{r}}\left(t\right)=\frac{1}{n_x}\sum _{i=1}^{n_x}\frac{\left|s_i^{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}\right(t)-s_i(t\left)\right|}{\left|s_i\right(t\left)\right|},$respectively. Figure 4 plots these errors against time for pressure and temperature using a POD basis with r = 43 vectors computed from the first $k=\mathrm{20,000}$ training snapshots. While both error measures for pressure remain low throughout the full time interval, both the projection errors and the ROM prediction errors for temperature increase significantly at the end of the training regime. The temperature profile is influenced by both the advective flow dynamics and the local chemical reactions, which in combination lead to a highly nonlinear and multiscale behaviour that is difficult to represent with the POD basis after the training period. However, while the ROM struggles to accurately predict the detailed temperature variations pointwise, it does adequately predict the general trends of temperature evolution beyond the training horizon. Figure 4 shows the time-averaged temperature profiles for the GEMS data and for an OpInf ROM, suggesting that the ROM captures the time-averaged behaviour of the temperature dynamics.

Figure 4. Spatially averaged relative errors in time for pressure and temperature (left) and temperature profiles for the GEMS dataset and an OpInf ROM averaged over $\mathrm{60,000}$ time steps (right). The basis comprises the r = 43 dominant singular vectors of $k=\mathrm{20,000}$ training snapshots.

Figure 5 plots CH${}_4$ and CO${}_2$ concentrations integrated over the spatial domain. These measures give a more global sense of the ROM predictive accuracy and the predicted chemical reaction rate. In each case, the ROMs are able to accurately re-predict the training data and capture much of the overall system behaviour in the prediction phase, with slightly more training error as the number of snapshots increases.

Figure 5. Integrated molar concentrations for CH${}_4$ and CO${}_2$, computed over the spatial domain for each point in time, for $k=\mathrm{10,000}$ (top row), $k=\mathrm{20,000}$ (middle row), and $k=\mathrm{30,000}$ (bottom row) training snapshots. These statistical features further highlight the effect of increasing k. As in Figure 3, r is chosen so that ${\mathcal{E}}_r>0.985$.

Comparison to POD-DEIM

We now compare regularised OpInf to a state-of-the-art nonlinear model reduction method that uses a least-squares Petrov-Galerkin POD projection coupled with the discrete empirical interpolation method (DEIM) (Chaturantabut and Sorensen 2010), as implemented for the same combustion problem in Huang et al. (2019) and Huang et al. (2018). This POD-DEIM method is intrusive–it requires nonlinear residual evaluations of the GEMS code at sparse discrete interpolation points. This also increases the computational cost of solving the POD-DEIM ROM in comparison to the OpInf ROM: integrating a POD-DEIM ROM with r = 70 for 6000 time steps of size $\delta t={10}^{-6}\mathrm{s}$ takes approximately 30 min on two nodes, each with two Haswell CPUs processors at $2.60\mathrm{G}\mathrm{H}\mathrm{z}$ and 20 cores per node; for OpInf, using Python 3.6.9 and a single CPU on an AMD EPYC 7702 64-core processor at $3.3\mathrm{G}\mathrm{H}\mathrm{z}$ with $2.1\mathrm{T}\mathrm{B}$ RAM, we solve Equation (5(5) $\underset{\mathbf{O}}{min}{∥\mathbf{D}{\mathbf{O}}^{\mathrm{\top }}-{\mathbf{R}}^{\mathrm{\top }}∥}_F^2+{∥\mathbf{\Gamma }{\mathbf{O}}^{\mathrm{\top }}∥}_F^2,$(5) ) with $k=\mathrm{20,000}$ training snapshots and r = 43 POD modes in approximately $0.6\mathrm{s}$ and integrate the resulting OpInf ROM for 60,000 time steps of size $\delta t={10}^{-7}\mathrm{s}$ in approximately $0.4\mathrm{s}$. While these measurements are made on different hardware, and though the execution time for POD-DEIM can be improved with optimal load balancing, the difference in execution times (30 min versus 1 s) is representative and illustrates one of the advantages over POD-DEIM of the polynomial form employed in the OpInf approach.

Figure 6 compares select GEMS outputs to POD-DEIM and OpInf ROM outputs, with each ROM trained on $k=\mathrm{20,000}$ training snapshots. As before, the OpInf ROM dimension $r=43$ is chosen such that ${\mathcal{E}}_r>0.985$; the POD-DEIM ROM, which uses an entirely different basis than the OpInf approach, requires r = 70 vectors to achieve the same level of cumulative energy. Both approaches maintain appropriate pressure oscillation frequencies, and while neither model accurately predicts the global species concentration dynamics after the training period, the OpInf ROM reconstructs the training data more faithfully than the POD-DEIM ROM. Huang et al. (2019) show similar results for the same POD-DEIM model with 1 ms of training and 1 ms of prediction; here we are using 2 ms of training and 4 ms of prediction. Note from Figures 3 and 5 that the OpInf ROMs achieve excellent prediction results for the 1 ms period following the training.

Figure 6. Pressure trace at monitor location 3 of Figure 1 (top) and spatially integrated O${}_2$ and CO${}_2$ molar concentrations (bottom), computed by GEMS, a POD-DEIM ROM, and an OpInf ROM. Both ROMs use $k=\mathrm{20,000}$ training snapshots with r chosen so that ${\mathcal{E}}_r>0.985$.

Figures 7 and 8 show, respectively, full-domain results for the temperature and molar concentration of CH${}_4$. The figures show the solution at time instants within the training regime, at the end of the training regime, and into the prediction regime. As with the point traces shown earlier, we see that the ROMs have impressive accuracy over the training region, but lose accuracy as they attempt to predict dynamics beyond the training horizon. However, many of the coherent features are reasonably predicted, especially the recirculation zone dynamics near the dump plane ($x=0$ in Figure 1) shown in the temperature fields. Significantly, both ROMs maintain appropriate temperature ranges throughout the prediction phase. The POD-DEIM ROM explicitly enforces such limits by reconstructing the full solution at each time step, constraining the temperature to a desired range, and projecting the result back to the reduced space (see Section IV.G of Huang et al. 2019); in contrast, the OpInf ROM selects a regularisation that results in bounded behaviour due to the criteria $\left|{\hat{q}}_i\right(t\left)\right|\le B$ (see Equation (7(7) $\left|q_i\right(t\left)\right|=\left|\sum _{j=1}^r{V}_{ij}{\hat{q}}_j\left(t\right)\right|\le \sum _{j=1}^r\left|V_{ij}\right|\left|{\hat{q}}_j\left(t\right)\right|\le B\sum _{j=1}^r\left|V_{ij}\right|,i=1,2,\dots ,n,$(7) )). In other words, POD-DEIM limits the temperature in the online phase, while OpInf builds a similar constraint into the offline phase.

Figure 7. Temperature fields produced by GEMS (left column), POD-DEIM (middle column), and OpInf (right column), where each ROMs uses $k=\mathrm{20,000}$ training snapshots, with r chosen so that ${\mathcal{E}}_r>0.985$. Each row shows results for a given time, with an increment of $0.0005\mathrm{s}$ between rows. The training period ends at $t=0.0170\mathrm{s}$ (fourth row); $t=0.0175\mathrm{s}$ (last row) is well into the prediction regime.

Figure 8. Molar concentrations of CH${}_4$ produced by GEMS (left column), POD-DEIM (middle column), and OpInf (right column), where each ROMs uses $k=\mathrm{20,000}$ training snapshots, with r chosen so that ${\mathcal{E}}_r>0.985$. Each row shows results for a given time, with an increment of $0.0005\mathrm{s}$ between rows. The training period ends at $t=0.0170\mathrm{s}$ (fourth row); $t=0.0175$ (last row) is well into the prediction regime.

Conclusions

The presented scientific machine learning approach is broadly applicable to problems where the governing equations are known but access to the high-fidelity simulation code is limited. The approach is computationally as accessible as black-box surrogate modelling while achieving the accuracy of intrusive projection-based model reduction. While the conclusions drawn from the numerical studies apply to the single-injector combustion example, they are relevant and likely apply to other problems. First, the quality and quantity of the training data are critical to the success of the method. Second, regularisation is essential to avoid overfitting. Third, achieving a low error over the training regime is not necessarily indicative of a reduced model with good predictive capability. This emphasises the importance of the training data. Fourth, physical quantities that exhibit large-scale coherent structures (e.g. pressure) are more accurately predicted by a reduced-order model than quantities that exhibit multiscale behaviour (e.g. temperature, species concentrations). Fifth, a significant advantage of the data-driven learning aspects of the approach is that the reduced model may be derived in any variables. This includes the possibility to include redundancy in the learning variables (e.g. to include both pressure and temperature). Overall, this paper illustrates the power and effectiveness of learning from data through the lens of physics-based models as a physics-grounded alternative to black-box machine learning.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This work has been supported in part by the Air Force Center of Excellence on Multi-Fidelity Modeling of Rocket Combustor Dynamics (Air Force Office of Scientific Research) under award FA9550-17-1-0195, and the U.S. Department of Energy AEOLUS MMICC center under award DE-SC0019303.

Notes

1 From here on we use $\hat{\mathbf{q}}\otimes \hat{\mathbf{q}}$ to indicate a compact Kronecker product with only the $\left(\begin{array}{c}r+1\\ 2\end{array}\right)=\frac{1}{2}r(r+1)$ unique quadratic terms (${\hat{q}}_1^2$, ${\hat{q}}_1{\hat{q}}_2$, ${\hat{q}}_1{\hat{q}}_3$, …); for matrices, the product is applied column-wise.

2 Or any other low-dimensional basis as desired.

3 See https://github.com/Willcox-Research-Group/rom-operator-inference-Python3.

4 See https://github.com/Willcox-Research-Group/ROM-OpInf-Combustion-2D for additional results.