Data-driven reduced-order models via regularized operator inference for a single-injector combustion process
aa r X i v : . [ c s . C E ] A ug Data-driven reduced-order models via regularized operator inferencefor a single-injector combustion process
Shane A. McQuarrie a , Cheng Huang b , Karen Willcox a a Oden Institute for Computational Engineering and Sciences, University of Texas at Austin,USA; b Department of Aerospace Engineering, University of Michigan, USA
ARTICLE HISTORY
Compiled August 10, 2020
ABSTRACT
This paper derives predictive reduced-order models for rocket engine combustiondynamics via a scientific machine learning approach, based on Operator Inference,that blends data-driven learning with physics-based modeling. The non-intrusive na-ture of the approach enables variable transformations that expose system structure.The specific contribution of this paper is to advance the formulation robustness andalgorithmic scalability of the approach. Regularization is introduced to the formu-lation to avoid over-fitting. The task of determining an optimal regularization isposed as an optimization problem that balances training error and stability of long-time integration dynamics. A scalable algorithm and open-source implementationare 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 regularization and an informed selection oflearning variables, the reduced-order models exhibit high accuracy in re-predictingthe training regime and acceptable accuracy in predicting future dynamics, whileachieving close to a million times speedup in computational cost. When compared toa state-of-the-art model reduction method, the Operator Inference models providethe same or better accuracy at approximately one thousandth of the computationalcost.
KEYWORDS
Model reduction; operator inference; scientific machine learning; combustion;data-driven reduced model.
1. Introduction
The emerging field of scientific machine learning brings together the perspectivesof physics-based modeling and data-driven learning. In the field of fluid dynamics,physics-based modeling and simulation have played a critical role in advancing scien-tific discovery and driving engineering innovation in domains as diverse as biomedicalengineering (Yin et al. 2010; Nordsletten et al. 2011), geothermal modeling (O’Sullivanet al. 2001; Cui et al. 2011), and aerospace (Spalart and Venkatakrishnan 2016). Theseadvances are based on decades of mathematical and algorithmic developments in com-putational fluid dynamics (CFD). Scientific machine learning builds upon these rig-orous physics-based foundations while seeking to exploit the flexibility and expressivemodeling capabilities of machine learning (Baker et al. 2019). This paper presents a
CONTACT Karen E. Willcox. Email: [email protected] cientific machine learning approach that blends data-driven learning with the theo-retical foundations of physics-based model reduction. This creates the capability tolearn predictive reduced-order models (ROMs) that provide approximate predictionsof complex physical phenomena while exhibiting several orders of magnitude compu-tational speedup over CFD.Projection-based model reduction considers the class of problems for which thegoverning 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 lowercomplexity and yields accurate solutions with reduced computation time. Projection-based approaches define a low-dimensional manifold on which the dominant dynamicsevolve. 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-fidelitysimulation data, as in the proper orthogonal decomposition (POD) (Lumley 1967;Sirovich 1987; Berkooz et al. 1993). The POD has been particularly successful influid dynamics, dating back to the early applications in unsteady flows and turbulencemodeling (Sirovich 1987; Deane et al. 1991; Gatski and Glauser 1992), and in unsteadyfluid-structure interaction (Dowell and Hall 2001).Model reduction methods have advanced to include error estimation (Veroy et al.2003; Veroy and Patera 2005; Grepl and Patera 2005; Rozza et al. 2008) and to ad-dress parametric and nonlinear problems (Barrault et al. 2004; Astrid et al. 2008;Chaturantabut and Sorensen 2010; Carlberg et al. 2013), yet the intrusive nature ofthe methods has limited their impact in practical applications. When legacy or com-mercial codes are used, as is often the case for CFD applications, it can be difficult orimpossible to implement classical projection-based model reduction. Black-box surro-gate modeling instead derives the ROM by fitting to simulation data; such methodsinclude response surfaces and Gaussian process models, long used in engineering, aswell as machine learning surrogate models. These methods are powerful and oftenyield good results, but since the approximations are based on generic data-fit repre-sentations, they are not equipped with the guarantees (e.g., stability guarantees, errorestimators) that accompany projection-based ROMs. Nonlinear system identificationtechniques seek to illuminate the black box by discovering the underlying physics of asystem from data (Brunton et al. 2016). However, when the governing dynamics areknown and simulation data are available, reduced models may be directly tailored tothe specific dynamics without access to the details of the large-scale CFD code.Building on the work in Swischuk et al. (2020), this paper presents a non-intrusivealternative to black-box surrogate modeling. We use the Operator Inference method ofPeherstorfer and Willcox (2016) to learn a ROM from simulation data; the structureof the ROM is defined by the known governing equations combined with the theory ofprojection-based model reduction. Our approach may be termed ‘glass-box modeling’,as the targeted dynamics are known via the governing partial differential equationsthat define the problem of interest but the inner workings of the CFD code are notaccessed. This paper extends our prior work by formally introducing regularization tothe approach, which is critical to avoid overfitting for problems with complex dynamics,as is the case for the combustion example considered here. A second contribution ofthis paper is a scalable implementation of the approach, which is available via anopen-source implementation. Section 2 presents the methodology and regularizationapproach and describes the scalable implementation. Section 3 presents numericalresults for a single-injector combustion problem and Section 4 concludes the paper.2 . Methodology
This section begins with an overview the Operator Inference approach in Section 2.1.Section 2.2 augments Operator Inference with a new regularization formulation, posedas an optimization problem, and presents a complete algorithm for regularizationselection and model learning. In Section 2.3, we discuss a scalable implementationof the algorithm, which can then be applied to CFD problems of high dimension.
We target problems governed by systems of nonlinear partial differential equations.Consider the governing equations of the system of interest written, after spatial dis-cretization, in semi-discrete formdd t q ( t ) = c + Aq ( t ) + H ( q ( t ) ⊗ q ( t )) + Bu ( t ) , q ( t ) = q , t ∈ [ t , t f ] , (1)where q ( t ) ∈ R n is the state vector discretized over n points in space at time t , u ( t ) ∈ R m denotes the m inputs at time t , typically related to boundary conditions or forcingterms, t and t f are respectively the initial and final time, and q is the given initialcondition. We refer to Eq. (1) as the full-order model (FOM) and note that it has beenwritten to have a polynomial structure: c ∈ R n are constant terms; Aq ( t ) are the termsthat are linear in the state q ( t ), with the discretized operator A ∈ R n × n ; H ( q ( t ) ⊗ q ( t ))are the terms that are quadratic in q ( t ), with H ∈ R n × n ; and Bu ( t ) are the terms thatare linear in the input u ( t ), with B ∈ R n × m . This polynomial structure arises in threeways: (1) it may be an attribute of the governing equations; (2) it may be exposed viavariable transformations; or (3) it may be derived by introducing auxiliary variablesthrough the process of lifting. As examples of each: (1) the incompressible Navier-Stokes equations have a quadratic form; (2) the Euler equations can be transformedto quadratic form by using pressure, velocity, and specific volume as state variables;(3) a nonlinear tubular reactor model with Arrhenius reaction terms can be writtenin quadratic form via the lifting transformation shown in Kramer and Willcox (2019)that introduces six auxiliary variables.A projection-based reduced-order model (ROM) of Eq. (1) preserves the polynomialstructure (Benner et al. 2015). Approximating the high-dimensional state q in a low-dimensional basis V ∈ R n × r , with r ≪ n , we write q ( t ) ≈ Vq r ( t ), where q r ( t ) ∈ R r is the reduced state. Using a Galerkin projection, this yields the intrusive ROM ofEq. (1):dd t q r ( t ) = c r + A r q r ( t ) + H r ( q r ( t ) ⊗ q r ( t )) + B r u ( t ) , q r ( t ) = V ⊤ q , t ∈ [ t , t f ] , where c r = V ⊤ c ∈ R r , A r = V ⊤ AV ∈ R r × r , H r = V ⊤ H ( V ⊗ V ) ∈ R r × r , and B r = V ⊤ B ∈ R r × m are the ROM operators corresponding to the FOM operators c , A , H , and B , respectively. The ROM is intrusive because computing these ROMoperators requires access to the discretized FOM operators, which typically entailsintrusive queries and/or access to source code.The non-intrusive Operator Inference (OpInf) approach proposed by Peherstorferand Willcox (2016) parallels the intrusive projection-based ROM setting, but learnsROMs from simulation data without direct access to the FOM operators. Recognizing3hat the intrusive ROM has the same polynomial form as Eq. (1), OpInf uses a data-driven regression approach to derive a ROM of Eq. (1) asdd t b q ( t ) = b c + b A b q ( t ) + b H ( b q ( t ) ⊗ b q ( t )) + b Bu ( t ) , b q ( t ) = V ⊤ q , t ∈ [ t , t f ] , (2)where b c ∈ R r , b A ∈ R r × r , b H ∈ R r × ( r +12 ), and b B ∈ R r × m are determined by solving adata-driven regression problem, and b q ( t ) ∈ R r is the state of the OpInf ROM. OpInf solves a regression problem to find reduced operators that yield the ROM thatbest matches projected snapshot data in a minimum-residual sense. Mathematically,OpInf solves the least-squares problemmin b c , b A , b H , b B k − X j =0 (cid:13)(cid:13)(cid:13)b c + b A b q j + b H ( b q j ⊗ b q j ) + b Bu j − ˙ b q j (cid:13)(cid:13)(cid:13) , (3)where { b q j } k − j =0 is the dataset used to drive the learning with b q j denoting a reduced-state snapshot at timestep j , { ˙ b q j } k − j =0 are the associated time derivatives, and { u j } k − j =0 is the collection of inputs corresponding to the data with u j ≡ u ( t j ). To generate thedataset, we employ the following steps: (1) Collect a set of k high-fidelity state snap-shots { x j } k − j =0 ⊂ R ˜ n by solving the original high-fidelity model at times { t j } k − j =0 withinputs { u j } k − j =0 . (2) Apply a variable transformation q j = T ( x j ) to obtain snapshotsof the transformed variables. Here T : R ˜ n → R n is the map representing an invert-ible transformation (e.g., from density to specific volume) or a lifting transformation(Qian et al. 2020). (3) Compute the proper orthogonal decomposition (POD) basis V of the transformed snapshots. (4) Project the transformed snapshots onto the PODsubspace as b q j = V ⊤ q j . (5) Estimate projected time derivative information { ˙ b q } k − j =0 .The training period [ t , t k − ] for which we have data is a subset of the full time domainof interest [ t , t f ]; results from the ROM over [ t k , t f ] will be entirely predictive.Eq. (3) can also be written in matrix form asmin O (cid:13)(cid:13)(cid:13) DO ⊤ − R ⊤ (cid:13)(cid:13)(cid:13) F , (4)where O = (cid:2) b c b A b H b B (cid:3) ∈ R r × d ( r,m ) , (unknown operators) D = (cid:2) k b Q ⊤ ( b Q ⊗ b Q ) ⊤ U ⊤ (cid:3) ∈ R k × d ( r,m ) , (known data) b Q = (cid:2) b q b q · · · b q k − (cid:3) ∈ R r × k , (snapshots) R = (cid:2) ˙ b q ˙ b q · · · ˙ b q k − (cid:3) ∈ R r × k , (time derivatives) U = (cid:2) u u · · · u k − (cid:3) ∈ R m × k , (inputs)and where d ( r, m ) = 1 + r + (cid:0) r +12 (cid:1) + m and k ∈ R k is the length- k column vector withall entries set to unity. The OpInf least-squares problem Eq. (3) is therefore linear in From here on we use b q ⊗ b q to indicate a compact Kronecker product with only the (cid:0) r +12 (cid:1) = r ( r + 1) uniquequadratic terms (ˆ q , ˆ q ˆ q , ˆ q ˆ q , . . . ); for matrices, the product is applied column-wise. Or any other low-dimensional basis as desired. b c , b A , b H , and b B .Note that the OpInf approach permits us to compute the ROM operators b c , b A , b H ,and b B without explicit access to the original high-dimensional operators c , A , H , and 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 derivinga classical intrusive ROM might be possible, the OpInf approach enables us to workwith variables other than those used for the original high-fidelity discretization. InSection 3 we will see the importance of this for a reacting flow application. We alsonote that under some conditions, OpInf recovers the intrusive POD ROM (Peherstorferand Willcox 2016; Peherstorfer 2019). Eq. (4) decouples into r independent linear least-squares problems, one for each ofthe rows of O (Peherstorfer and Willcox 2016). Each problem is generally overdeter-mined, but is also noisy due to error in the numerically estimated time derivatives R , model mis-specification (e.g., if the system is not truly quadratic), and truncatedPOD modes that leave some system dynamics unresolved. The ROMs resulting fromEq. (4) can thus suffer from overfitting the operators to the data and therefore exhibitpoor predictive performance over the time domain of interest [ t , t f ].To avoid overfitting, we introduce a Tikhonov regularization to the sub-problemsin Eq. (4), which then becomemin o i k Do i − r i k + k Γo i k = (cid:13)(cid:13)(cid:13)(cid:13)(cid:20) DΓ (cid:21) o i − (cid:20) r i (cid:21)(cid:13)(cid:13)(cid:13)(cid:13) , i = 1 , . . . , r, (5)where o i ∈ R d ( r,m ) is the i th row of O (the i th column of O ⊤ ), r i ∈ R k is the i th rowof R , and Γ ∈ R d ( r,m ) × d ( r,m ) is a full-rank regularizer. Each regularized sub-problemin Eq. (5) admits a unique solution since, by construction, the augmented data matrix (cid:2) D ⊤ Γ ⊤ (cid:3) ⊤ ∈ R ( k + d ( r,m )) × d ( r,m ) is taller than it is wide and has full column rank.An L regularizer Γ = λ I , λ > I the identity matrix, penalizes each entry ofthe inferred ROM operators b c , b A , b H , and b B , thereby driving the ROM toward theglobally stable system dd t b q ( t ) = . An appropriate value for the scalar hyperparameter λ , which balances the minimization between the data fit and the regularization, mustbe chosen with care. The ideal λ produces a ROM that minimizes some error metricover the full time domain [ t , t f ]; however, since data are only available for the smallertraining domain [ t , t k − ], we choose λ so that the resulting ROM minimizes error over[ t , t k − ] while maintaining a bound on the integrated POD coefficients over [ t , t f ].That is, we require that any b q = (cid:2) ˆ q ˆ q · · · ˆ q r (cid:3) ⊤ produced by the ROM trainedwith parameter λ satisfy max i | ˆ q i | ≤ B for some B >
0. This in turn ensures a boundon the magnitude of the entries of the high-dimensional state q = V b q : | q i | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r X j =1 V ij ˆ q j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ r X j =1 | V ij | | ˆ q j | ≤ B r X j =1 | V ij | , i = 1 , , . . . , n, where V ij is the i th element of the j th POD basis vector. Note that B could be chosenwith the intent of imposing a particular bound on the | q i | since the sums P rj =1 | V ij | can be precomputed. 5lgorithm 1 details our regularized OpInf procedure, in which we choose B as a mul-tiple of the maximum absolute entry of the projected training data b Q . This particularstrategy for selecting λ could be replaced with a cross-validation or resampling gridsearch technique, but our experiments in this vein did not produce robust results. Thetraining error k b Q − e Q : , : k k in step 14 may compare b Q and e Q : , : k directly in the reducedspace (e.g., with a matrix norm or an L p ([ t , t k − ]) norm), or it may be replaced witha more targeted comparison of some quantity of interest. Algorithm 1
Regularized Operator Inference with regularization selection procedure RegOpInf (snapshots X ∈ R ˜ n × k , inputs U ∈ R m × k , final time t f > t ,invertible map T : R ˜ n → R n , reduced dimension r ∈ N , regularization parameterinterval [ λ low , λ high ] ⊂ R + , bound margin τ ≥ Q ← T ( X ) ⊲ Map native variables to learning variables (columnwise). V ← pod( Q , r ) ⊲ Compute a rank- r POD basis from transformed snapshots. b Q ← V ⊤ Q ⊲ Project snapshots to the r -dimensional POD subspace. R ← dd t b Q ⊲ Estimate time derivatives of the projected snapshots. b q ← b Q : , ⊲ The initial condition for Eq. (2) is the first projected snapshot. B ← τ max i,j | b Q ij | ⊲ Select a bound to require for integrated POD modes. procedure TrainError ( λ ) b c , b A , b H , b B ← solve Eq. (5) with data b Q , R , U and regularizer λ I e Q ← integrate Eq. (2) with b c , b A , b H , b B from b q over [ t , t f ] if max i,j | e Q ij | > B then return ∞ ⊲ Disqualify ROMs that violate the bound. else return k b Q − e Q : , : k k ⊲ Compute ROM error over [ t , t k − ]. λ ∗ ← argmin λ ∈ [ λ low ,λ high ] TrainError ( λ ) ⊲ Find the λ with least training error. b c , b A , b H , b B ← solve Eq. (5) with data b Q , R , U and optimal regularizer λ ∗ I return b c , b A , b H , b B The steps of Algorithm 1 are highly modular and amenable to large problems. Thevariable transformation Q = T ( X ) of step 2 consists of O ( nk ) elementary computa-tions on the original snapshot matrix X ∈ R ˜ n × k . To compute the rank- r POD basis V ∈ R n × r in step 3, we use a randomized singular-value decomposition (SVD) algo-rithm requiring O ( rnk + r ( n + k )) operations (Halko et al. 2011); since r ≪ k, n , theleading order behavior is O ( rnk ). The projection b Q = V ⊤ Q in step 4 costs 2 nk oper-ations. Note that b Q ∈ R r × k is small compared to Q ∈ R n × k , as typically r ∼ –10 and n ∼ –10 . The time derivatives R ∈ R r × k in step 5 may be provided by thefull-order solver or estimated, e.g., with finite differences of the columns of b Q . In thelatter case, the cost is O ( rk ). The computational cost of steps 2–6 is therefore O ( rnk ).Solving the regularized problems Eq. (5) in steps 9 and 16 requires computing theSVD of the augmented data matrix (cid:2) D ⊤ Γ ⊤ (cid:3) ⊤ — equivalently, the eigendecom-position of D ⊤ D + Γ ⊤ Γ ∈ R d ( r,m ) × d ( r,m ) — which costs O ( d ( r, m ) ) = O ( r ) oper-ations (Demmel 1997). This SVD is then applied to each of the augmented vectors6 igure 1. The computational domain for the single-injector combustion problem. On the left, monitor lo-cations for numerical results. On the right, a typical temperature field demonstrating the complexity andnonlinear nature of the problem. (cid:2) r ⊤ i ⊤ (cid:3) ⊤ to produce the unique solution. The ROM integration in step 10 can becarried out with any time-stepping scheme; for explicit methods, evaluating the ROMat a single point, i.e., computing the right-hand side of Eq. (2), costs O ( r ) opera-tions. The total cost of evaluating the subroutine TrainError : R → R is therefore O ( r ) which, importantly, is independent of n . The minimization in step 15 is carriedout with a bisection-type search method, which enables fewer total evaluations of thesubroutine than a blind grid search. However, an initial coarse grid search is useful foridentifying an appropriate search interval [ λ low , λ high ].
3. Results
This section applies regularized OpInf to a single-injector combustion problem, studiedpreviously by Swischuk et al. (2020), on the two-dimensional computational domainshown in Figure 1. Section 3.1 describes the governing dynamics, a set of high-fidelitydata obtained from a CFD code, and the variable transformations used to producetraining data for learning reduced models with Algorithm 1. The resulting OpInf ROMperformance is analyzed in Section 3.2 and compared to a state-of-the-art intrusivemodel reduction method in Section 3.3.
The combustion dynamics for this problem are governed by conservation laws ∂~q c ∂t + ∇ · ( ~K − ~K v ) = ~S, (6)where ~q c = (cid:2) ρ ρv x ρv y ρe ρY ρY ρY ρY (cid:3) ⊤ are the conservative vari-ables, ~K are the inviscid flux terms, ~K v are the viscous flux terms, and ~S are thesource terms. Here ρ is the density [ kgm ], v x and v y are the x and y velocity [ ms ], e isthe total energy [ Jm ], and Y ℓ is the ℓ th species mass fraction with ℓ = 1 , , ,
4. Thechemical species are CH , O , H O, and CO , which follow a global one-step chemicalreaction CH + 2O → CO + 2H O (Westbrook and Dryer 1981). See Harvazinskiet al. (2015) for more details on the governing equations.At the downstream end of the combustor, we impose a non-reflecting boundary7ondition while maintaining the chamber pressure via p back ( t ) = p back,ref (1 + 0 . πf t )) (7)where p back,ref = 10 Pa and f = 5 ,
000 Hz. The top and bottom wall boundaryconditions are no-slip conditions, and for the upstream boundary we impose a constantmass flow at the inlets.Swischuk et al. (2020) show that if the governing equations (6) are transformed tobe 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 augmentedstate ~q = (cid:2) p v x v y T ξ c c c c (cid:3) ⊤ where p is the pressure [Pa], T is thetemperature [K], ξ = 1 /ρ is the specific volume [ m kg ], and c , . . . , c are the speciesmolar concentrations [ kmolm ] given by c ℓ = ρY ℓ /M ℓ with M ℓ the molar mass of the ℓ th species [ gmol ]. As shown in Swischuk et al. (2020), the equations for v x , v y , and ξ are exactly quadratic in ~q , while the remaining equations are quadratic with somenon-polynomial terms in ~q . Note that, differently from Swischuk et al. (2020), ~q hereis chosen to contain specific volume, pressure, and temperature, even though onlytwo of the three quantities are needed to fully define the high-fidelity model (andthe equation of state then defines the third). We augment the learning variables inthis way because doing so exposes the quadratic form while also directly targeting thevariables that are of primary interest for assessing ROM performance. In particular, theresulting ROMs provide more accurate predictions of temperature when temperatureis included explicitly as a learning variable. We can do this since the transformationsare applied only to the snapshot data, not to the CFD model itself. This flexibility isa major advantage of the non-intrusive OpInf approach in comparison to traditionalintrusive projection-based model reduction methods.To generate high-fidelity training data, we use the finite-volume based GeneralEquation and Mesh Solver (GEMS) (Harvazinski et al. 2015) to solve for the variables (cid:2) p v x v y T Y Y Y Y (cid:3) over n x = 38 ,
523 cells, resulting in snapshots with8 n x = 308 ,
184 entries each. The snapshots are computed for 60,000 time steps beyondthe initial condition with a temporal discretization of δt = 10 − s, from t = 0 .
015 sto t f = 0 .
021 s. The computational cost of computing this dataset is approximately1,200 CPU hours on two computing nodes, each of which contains two Haswell CPUsat 2 .
60 GHz and 20 cores per node.Scaling is an essential aspect of successful model reduction and is particularly criticalfor this problem due to the wide range of scales across variables. After transformingthe GEMS snapshot data to the learning variables ~q , the species molar concentrationsare scaled to [0 , − , rom operator inference Python package, which is built on NumPy, SciPy, and scikit-learn (Walt et al. 2011; Virtanen et al. 2020;Buitinck et al. 2013). The time derivatives in step 5 of Algorithm 1 are estimated withfourth-order finite differences, and the least-squares problems in step 9 are solved bythe divide-and-conquer LAPACK routine DGELSD. The learned ROMs are integratedin step 10 with the explicit, adaptive, fourth-order Runge-Kutta scheme RK45, and See https://github.com/Willcox-Research-Group/rom-operator-inference-Python3 . Singular Value Index − S i n g u l a r V a l u e s k = 10000 k = 20000 k = 30000 Figure 2.
The POD singular values for varying sizeof the snapshot training set.
Table 1.
The basis size r required to reacha given cumulative energy level E r increaseslinearly with the number of snapshots in thetraining set, k . Cumulative energy E r .
985 0 .
990 0 . k = 10 , r = 22 r = 28 r = 37 k = 20 , r = 44 r = 54 r = 73 k = 30 , r = 68 r = 83 r = 112 the error evaluation of step 14 uses the L ([ t , t k − ]) norm in the reduced space. Tominimize the scalar function TrainError ( λ ) in step 15, we use Brent’s method (Brent2002) with enforced bounds [ λ low , λ high ]. The code and details are publicly availableat https://github.com/Willcox-Research-Group/ROM-OpInf-Combustion-2D . We study the sensitivity of our approach to the training data by varying the numberof snapshots used to compute the POD basis and learn the OpInf ROM. Specifically,we consider the three cases where we use the first k = 10 , k = 20 , k = 30 ,
000 snapshots from GEMS as training data sets. In each case we computethe POD basis, and to select an appropriate reduced dimension r , we compute thecumulative energy based on the POD singular values: E r = (cid:16)P rj =1 σ j (cid:17).(cid:16)P kj =1 σ j (cid:17) ,where { σ j } kj =1 are the singular values of the learning variable snapshot matrix Q (see Figure 2). Specifically, we choose the minimal r such that E r is greater than afixed threshold energy. Table 1 shows that r is increasing linearly with the number ofsnapshots, indicating that the basis is not being saturated as additional informationis incorporated. This is an indication of the challenging nature of this application, dueto the rich and complex dynamics. It is also an indication of the importance of havingsufficient training data.Figure 3 plots the GEMS and OpInf ROM results for pressure, temperature and x -velocity predictions over time at the monitor locations in Figure 1. While it canbe misleading to assess accuracy based on predictions at a single spatial point, theseplots reveal several representative insights. First, we see the importance of the trainingdata—as the amount of training data increases, the ROM predictions change signifi-cantly and generally (but not always) improve. This is yet another indication of thecomplexity of the dynamics we are aiming to approximate. Second, we see that goodROM performance over the training region is not sufficient to produce good predictivebehavior. In particular, the ROM with only k = 10 ,
000 training snapshots can accu-rately re-predict the training dynamics, but shows large errors as it predicts furtherbeyond the training horizon. Third, we see that pressure and x -velocity are better ap-proximated by the ROM than is temperature. The effects of the 5,000 Hz downstreampressure forcing are clearly visible in the pressure and x -velocity traces. The temper- See https://github.com/Willcox-Research-Group/ROM-OpInf-Combustion-2D for additional results. . . × r = 22 Location 1
GEMS ROM × r = 22 Location 2 × r = 22 Location 3 . . P r e ss u r e [ P a ] × r = 44 T e m p e r a t u r e [ K ] × r = 44 x - v e l o c i t y [ m / s ] × r = 44 .
015 0 .
017 0 .
019 0 . Time [s] . . × r = 68 .
015 0 .
017 0 .
019 0 . Time [s] × r = 68 .
015 0 .
017 0 .
019 0 . Time [s] × r = 68 Figure 3.
Traces of p , T , and v x through time at monitor locations 1, 2, and 3 of Figure 1, respectively, for k = 10 ,
000 (top row), k = 20 ,
000 (middle row), and k = 30 ,
000 (bottom row) training snapshots. The verticalblack lines separate the training and prediction periods. For each choice of k the error over the training domainis low, but increasing k has a significant impact on the error in the testing domain. Here r is chosen in eachcase so that E r > . ature dynamics, on the other hand, are more irregular, and exhibit short-term swingsexceeding 1,000 K. This is because the temperature profile is influenced by both theflow dynamics and the local chemical reactions, which in combination lead to a highlynonlinear and multiscale behavior. The ROM cannot accurately predict the detailedvariations, but with k = 30 ,
000 training snapshots, the ROM does an adequate job ofpredicting the general trends of temperature evolution beyond the training horizon.Figure 4 plots integrated quantities as a function of time: the temperature averagedover the spatial domain, and the CH and O concentrations integrated over the spatialdomain. These measures give a more global sense of the ROM predictive accuracy. Ineach case, the ROMs are able to accurately re-predict the training data and capturemuch of the overall system behavior in the prediction phase. We now compare regularized OpInf to a state-of-the-art nonlinear model reductionmethod that uses a least-squares Petrov-Galerkin POD projection coupled with thediscrete empirical interpolation method (DEIM) (Chaturantabut and Sorensen 2010),as implemented for the same combustion problem in Huang et al. (2019, 2018). ThisPOD-DEIM method is intrusive—it requires nonlinear residual evaluations of theGEMS code at sparse discrete interpolation points. This also increases the compu-tational cost of solving the POD-DEIM ROM in comparison to the OpInf ROM:integrating a POD-DEIM ROM with r = 70 for 6,000 time steps of size δt = 10 − stakes approximately 30 minutes on two nodes, each with two Haswell CPUs processorsat 2 .
60 GHz and 20 cores per node; for OpInf, using Python 3 . . . . k = 20 ,
000 training snapshots and r = 44 POD modes in approximately3 . δt = 10 − sin approximately 0 . . . × r = 22 GEMS ROM . . × r = 22 . . × r = 22 . . Sp a t i a ll y A v e r ag e d T e m p e r a t u r e [ K ] × r = 44 . . C H C o n ce n tr a t i o n I n t e g r a l[ k m o l / m ] × r = 44 . . O C o n ce n tr a t i o n I n t e g r a l[ k m o l / m ] × r = 44 .
015 0 .
017 0 .
019 0 . Time [s] . . × r = 68 .
015 0 .
017 0 .
019 0 . Time [s] . . × r = 68 .
015 0 .
017 0 .
019 0 . Time [s] . . × r = 68 Figure 4.
Spatially averaged temperature and molar concentration integrals for CH and O , computed overthe spatial domain for each point in time, for k = 10 ,
000 (top row), k = 20 ,
000 (middle row), and k = 30 , k .As in Figure 3, r is chosen so that E r > . .
015 0 .
017 0 .
019 0 . Time [s] . . P r e ss u r e [ P a ] × Location 4
GEMS POD-DEIM, r = 70 OpInf, r = 44 .
015 0 .
017 0 .
019 0 . Time [s] . . Sp a t i a ll y A v e r ag e d T e m p e r a t u r e [ K ] × Figure 5.
Pressure trace at monitor location 4 of Figure 1 and the spatially averaged temperature, computedby GEMS, a POD-DEIM ROM, and an OpInf ROM. Both ROMs use k = 20 ,
000 training snapshots with r chosen so that E r > . and though the execution time for POD-DEIM can be improved with optimal loadbalancing, the difference in execution times (30 minutes versus 0.5 seconds) is rep-resentative and illustrates one of the advantages over POD-DEIM of the polynomialform employed in the OpInf approach.Figure 5 compares select GEMS outputs to POD-DEIM and OpInf ROM outputs,with each ROM trained on k = 20 ,
000 training snapshots and the ROM dimension r chosen to achieve a cumulative energy E r of 0 . . The figures show the solution at time instantswithin the training regime, at the end of the training regime, and into the predictionregime. As with the point traces shown earlier, we see that the ROMs have impres-11 a) t = t , = 0 . t = t , = 0 . t = t , = 0 . Figure 6.
Pressure p produced by GEMS (top row), POD-DEIM (middle row), and OpInf (bottom row), attimes within the training regime (left column), at the end of the training regime (middle column), and into theprediction regime (right column), with k = 20 ,
000 training snapshots and r chosen so that E r > . sive accuracy over the training region, but lose accuracy as they attempt to predictdynamics beyond the training horizon. Pressure is again well predicted, but the tem-perature and CH concentration fields predicted at t = 0 . x = 0 in Figure 1) shown inthe temperature fields.
4. Conclusions
The presented scientific machine learning approach is broadly applicable to problemswhere the governing equations are known but access to the high-fidelity simulationcode is limited. The approach is computationally as accessible as black-box surrogatemodeling while achieving the accuracy of intrusive projection-based model reduction.While the conclusions drawn from the numerical studies apply to the single-injectorcombustion example, they are relevant and likely apply to other problems. First, thequality and quantity of the training data are critical to the success of the method.Second, regularization is essential to avoid overfitting. Third, achieving a low errorover the training regime is not necessarily indicative of a reduced model with goodpredictive capability. This emphasizes the importance of the training data. Fourth,physical quantities that exhibit large-scale coherent structures (e.g., pressure) are moreaccurately predicted by a reduced-order model than quantities that exhibit multiscalebehavior (e.g., temperature, species concentrations). Fifth, a significant advantageof the data-driven learning aspects of the approach is that the reduced model maybe derived in any variables. This includes the possibility to include redundancy inthe learning variables (e.g., to include both pressure and temperature). Overall, thispaper illustrates the power and effectiveness of learning from data through the lens ofphysics-based models as a physics-grounded alternative to black-box machine learning.12 a) t = t , = 0 . t = t , = 0 . t = t , = 0 . Figure 7.
Temperature T produced by GEMS (top row), POD-DEIM (middle row), and OpInf (bottom row),at times within the training regime (left column), at the end of the training regime (middle column), and intothe prediction regime (right column), with k = 20 ,
000 training snapshots and r chosen so that E r > . t = t , = 0 . t = t , = 0 . t = t , = 0 . Figure 8.
Molar concentrations of CH produced by GEMS (top row), POD-DEIM (middle row), and OpInf(bottom row), at times within the training regime (left column), at the end of the training regime (middlecolumn), and into the prediction regime (right column), with k = 20 ,
000 training snapshots and r chosen sothat E r > . cknowledgements This work has been supported in part by the Air Force Center of Excellence on Multi-Fidelity Modeling of Rocket Combustor Dynamics under award FA9550-17-1-0195, andthe US Department of Energy AEOLUS MMICC center under award DE-SC0019303.
References repl M, Patera A. 2005. A posteriori error bounds for reduced-basis approximations ofparametrized parabolic partial differential equations. ESAIM-Mathematical Modelling andNumerical Analysis (M2AN). 39(1):157–181.Halko N, Martinsson PG, Tropp JA. 2011. Finding structure with randomness: Probabilisticalgorithms for constructing approximate matrix decompositions. SIAM Review. 53(2):217–288.Harvazinski ME, Huang C, Sankaran V, Feldman TW, Anderson WE, Merkle CL, Talley DG.2015. Coupling between hydrodynamics, acoustics, and heat release in a self-excited unstablecombustor. Physics of Fluids. 27(4):045102.Huang C, Duraisamy K, Merkle CL. 2019. Investigations and improvement of robustness ofreduced-order models of reacting flow. AIAA Journal. 57(12):5377–5389.Huang C, Xu J, Duraisamy K, Merkle C. 2018. Exploration of reduced-order models for rocketcombustion applications. In: 2018 AIAA Aerospace Sciences Meeting; Orlando, FL. PaperAIAA-2018-1183.Kramer B, Willcox K. 2019. Nonlinear model order reduction via lifting transformations andproper orthogonal decomposition. AIAA Journal. 57(6):2297–2307.Lumley J. 1967. The structures of inhomogeneous turbulent flow. Atmospheric Turbulence andRadio Wave Propagation:166–178.Nordsletten D, McCormick M, Kilner P, Hunter P, Kay D, Smith N. 2011. Fluid–solid couplingfor the investigation of diastolic and systolic human left ventricular function. InternationalJournal for Numerical Methods in Biomedical Engineering. 27(7):1017–1039.O’Sullivan M, Pruess K, Lippmann M. 2001. State of the art of geothermal reservoir simulation.Geothermics. 30:395–429.Peherstorfer B. 2019. Sampling low-dimensional Markovian dynamics for pre-asymptoticallyrecovering reduced models from data with operator inference. arXiv:190811233.Peherstorfer B, Willcox K. 2016. Data-driven operator inference for nonintrusive projection-based model reduction. Computer Methods in Applied Mechanics and Engineering. 306:196–215.Qian E, Kramer B, Peherstorfer B, Willcox K. 2020. Lift & Learn: Physics-informed machinelearning for large-scale nonlinear dynamical systems. Physica D: Nonlinear Phenomena.406:132401.Rozza G, Huynh D, Patera A. 2008. Reduced basis approximation and a posteriori error esti-mation for affinely parametrized elliptic coercive partial differential equations: applicationto transport and continuum mechanics. Archives of Computational Methods in Engineering.15(3):229–275. Available from: http://dx.doi.org/10.1007/s11831-008-9019-9 .Sirovich L. 1987. Turbulence and the dynamics of coherent structures. I. Coherent structures.Quarterly of Applied Mathematics. 45(3):561–571.Spalart P, Venkatakrishnan V. 2016. On the role and challenges of cfd in the aerospace industry.The Aeronautical Journal. 120(1223):209.Swischuk R, Kramer B, Huang C, Willcox K. 2020. Learning physics-based reduced-ordermodels for a single-injector combustion process. AIAA Journal. 58(6):2658–2672.Veroy K, Patera A. 2005. Certified real-time solution of the parametrized steady incompressibleNavier-Stokes equations: Rigorous reduced-basis a posteriori error bounds. InternationalJournal for Numerical Methods in Fluids. 47:773–788.Veroy K, Prud’homme C, Rovas D, Patera A. 2003. A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differentialequations. In: Proceedings of the 16th AIAA Computational Fluid Dynamics Conference;Orlando, FL. Paper AIAA-2003-3847.Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, Burovski E,Peterson P, Weckesser W, Bright J, et al. 2020. SciPy 1.0: Fundamental Algorithms forScientific Computing in Python. Nature Methods. 17:261–272.Walt Svd, Colbert SC, Varoquaux G. 2011. The NumPy array: A structure for efficient nu-merical computation. Computing in Science & Engineering. 13(2):22–30.Westbrook CK, Dryer FL. 1981. Simplified reaction mechanisms for the oxidation of hydro- arbon fuels in flames. Combustion Science and Technology. 27(1-2):31–43.Yin Y, Choi J, Hoffman EA, Tawhai MH, Lin CL. 2010. Simulation of pulmonary air flow witha subject-specific boundary condition. Journal of Biomechanics. 43(11):2159–2163.arbon fuels in flames. Combustion Science and Technology. 27(1-2):31–43.Yin Y, Choi J, Hoffman EA, Tawhai MH, Lin CL. 2010. Simulation of pulmonary air flow witha subject-specific boundary condition. Journal of Biomechanics. 43(11):2159–2163.