Conservative model reduction for finite-volume models

Introduction

The finite-volume method is commonly employed for discretizing systems of partial differential equations (PDEs) that associate with conservation laws, especially those in fluid dynamics. Rather than operating on the strong form of the PDE, the finite-volume method operates on the integral form of the PDE to numerically enforce conservation over each control volume comprising the computational mesh. Thus, conservation is the primary problem structure imposed by finite-volume discretizations; this contrasts with other discretization techniques that aim to preserve other properties, e.g., variational principles in the case of the finite-element discretizations.

Unfortunately, the computational burden imposed by high-fidelity finite-volume models is often prohibitive, as (1) the fine spatiotemporal resolution typically needed to ensure a verified, validated computational model can lead to extremely large-scale models whose simulations consume months on supercomputers, and (2) many engineering problems are real time or many query in nature. Such problems require the (parameterized) computational model to be simulated rapidly either due to a strict time-to-solution constraint in the case of real-time problems (e.g., model predictive control) or due to the need for hundreds or thousands of simulations in the case of many-query problems (e.g., statistical inversion).

Reduced-order models (ROMs) have been developed to mitigate this burden. These techniques first perform an offline stage during which they execute computationally costly training tasks (e.g., simulating the high-fidelity model for several parameter instances) to compute a low-dimensional ‘trial’ basis for the state. Next, these methods execute a computationally inexpensive online stage during which they rapidly compute approximate solutions for different points in the parameter space by projection: they compute solutions in the span of the trial basis while enforcing the high-fidelity model residual to be orthogonal to the subspace spanned by a low-dimensional ‘test’ basis. In the presence of nonlinearities, these techniques also introduce ‘hyper-reduction’ approximations to ensure the cost of simulating the ROM is independent of the high-fidelity-model dimension.

The most popular model-reduction approach for nonlinear dynamical systems such as those arising from finite-volume discretizations is Galerkin projection , wherein the test basis is set to be equal to the trial basis. The trial basis is often computed via proper orthogonal decomposition (POD) , but it can also be computed via the reduced-basis method; see Refs. , which apply the classical reduced-basis method to finite-volume problems. More recently, the least-squares Petrov–Galerkin (LSPG) projection method was proposed, which has been computationally demonstrated to generate accurate and stable responses for turbulent, compressible flow problems on which Galerkin projection yielded unstable responses. Unfortunately, neither Galerkin nor LSPG projection directly preserves important problem structure related to conservation laws or finite-volume models.

To address this, alternative projection techniques have been developed for improving the performance of reduced-order models when applied to conservation laws, particularly those appearing in fluid dynamics. These include stabilizing inner products applied to finite-difference and finite-element discretizations ; introducing dissipation via closure models or numerical dissipation ; performing nonlinear Galerkin projection based on approximate inertial manifolds ; including a pressure-term representation ; modifying the POD basis by including many modes (such that dissipative modes are captured), changing the norm , enabling adaptivity , or including basis functions that resolve a range of scales or respect the attractor’s power balance ; modifying the projection by adopting a constrained Galerkin , constrained Petrov–Galerkin , or $L^{1}$-norm minimizing projection ; developing approaches tailored to the incompressible Navier–Stokes equations by introducing stabilizations based on supremizer-enriched velocity spaces and a pressure Poisson equation or by modifying the Galerkin projection ; and improving the ROM’s ability to capture shocks . Among these contributions, only a subset is applicable to finite-volume discretizations. Further, no model-reduction method to date has been developed to preserve the structure intrinsic to finite-volume models: conservation. In particular, none of the above methods ensures that conservation holds over any subset of the computational domain, which can lead to spurious growth or dissipation of quantities that should be conserved in principle.

To this end, this work proposes a novel projection scheme for finite-volume models that ensures the reduced-order model is conservative over subdomains of the problem. The approach leverages the minimum-residual formulation of both Galerkin and least-squares Petrov–Galerkin projection by equipping their associated optimization problems with (generally nonlinear) equality constraints that explicitly enforce conservation over subdomains. The resulting conservative reduced-order models can be expressed as the solution to time-dependent saddle-point problems. The approach does not rely on a particular choice of reduced basis, although the reduced basis can affect feasibility of the associated optimization problems. New contributions in this work include:

Conservative Galerkin (Section 4.2) and conservative LSPG (Section 4.3) projection techniques, which ensure that the reduced-order models are conservative over subdomains of the original computational mesh. These methods are equipped with

techniques for handling infeasible constraints (Section 4.4), and

hyper-reduction techniques that respect the underlying finite-volume discretization to handle nonlinearities in the flux and source terms (Section 4.5).

demonstration that conservative Galerkin projection and time discretization are commutative (Theorem 4.3),

sufficient conditions for feasibility of conservative Galerkin (Proposition 5.1) and conservative LSPG (Proposition 5.2) projection,

conditions under which conservative Galerkin and conservative LSPG projection are equivalent (Theorem 5.1), and

a posteriori bounds (Section 5.3) for the error in the quantities conserved over subdomains (Theorem 5.3), in the null space (Lemma 5.1) and row space (Lemma 5.2) of the constraints, in the full state (Theorem 5.2), and in the conserved quantities (Lemma 5.3 and Theorem 5.3).

Numerical experiments on a parameterized quasi-1D Euler equation associated with modeling inviscid compressible flow in a converging–diverging nozzle (Section 6). These experiments demonstrate the merits of the proposed method and illustrate the importance of ensuring reduced-order models are globally conservative.

We remark that this work was first presented publically at the “Recent Developments in Numerical Methods for Model Reduction” workshop at the Institut Henri Poincaré on November 10, 2016.

Other works have also explored formulating reduced-order models that associate with constrained optimization problems. Zimmermann et al. equip equality ‘aerodynamic constraints’ to ROMs applied to steady-state external flows, where the constraints associate with matching experimental data or target performance metrics in a design setting. Recently, Reddy et al. propose equipping the time-discrete Galerkin ROM with inequality constraints that enforce solution positivity or a bound on the gas-void fraction. Relatedly, Fick et al. proposed a modified Galerkin optimization problem applicable to the incompressible Navier–Stokes equations, where the inequality constraints associate with bounds on the generalized coordinates; these bounds correspond to the extreme values of the generalized coordinates arising during the training simulations.

The remainder of this paper is organized as follows. Section 2 describes finite-volume discretizations of conservation laws (Section 2.1) discretized in time with a linear multistep scheme (Section 2.2). Section 3 describes the (standard) nonlinear model-reduction methods of Galerkin (Section 3.1) and LSPG (Section 3.2) projection, as well as their hyper-reduced variants (Section 3.3) and interpretations when applied to finite-volume models (Section 3.4). Section 4 describes the proposed methodology, which is based on enforcing conservation over decompositions (Section 4.1) of the computational mesh. Here, Section 4.2 describes the proposed conservative Galerkin projection technique, Section 4.3 describes the proposed conservative LSPG projection method, Section 4.4 describes approaches for handling constraint infeasibility, Section 4.5 describes the application of hyper-reduction to the constraints that respects the underlying finite-volume discretization, and Section 4.6 describes briefly how the quantities required for the proposed ROMs can be constructed from training data. Next, Section 5 performs analysis, including proving sufficient conditions for feasibility (Section 5.1), providing conditions under with the conservative Galerkin and conservative LSPG models are equivalent (Section 5.2), and deriving local a posteriori error analysis (Section 5.3). Section 6 demonstrates the benefits off the proposed method on a parameterization of the one-dimensional (compressible) Euler equations applied to a converging–diverging nozzle. Finally, Section 7 concludes the paper.

Finite-volume discretization

This work considers parameterized systems of conservation laws. In integral form, the associated governing equations correspond to

where $\delta_{ij}$ denotes the Kronecker delta. In matrix form, Eq. (2.7) becomes

This formulation will be exploited in Section 4, where we introduce the proposed method.

The full-order model ODE (2.5) is typically the starting point for developing reduced-order models for nonlinear dynamical systems. In this work, we exploit the particular structure underlying the dynamical system arising from the definitions of the state (2.3) and velocity (2.4).

2 Time discretization

A time discretization is required to solve (2.5) numerically. For simplicity, we restrict the focus in this work to linear multistep schemes, although other time integrators could be considered; see, e.g., Ref. , which develops LSPG reduced-order models for explicit, fully implicit, and diagonally implicit Runge–Kutta schemes. Applying a linear $k$-step method to numerically solve Eq. (2.5) at a given parameter instance $\boldsymbol{\mu}\in\mathcal{D}$ can be written as

Adams methods. Adams methods consider the integrated form of Eq. (2.5)

and apply a polynomial approximation to the integrand. In particular, the $p$th-order Adams scheme employs coefficients $\alpha_{0}=1$, $\alpha_{1}=-1$, and $\alpha_{j}=0$, $j>1$ and coefficients $\beta_{j}$ that associate with a polynomial interpolation of the integrand. In the explicit ($\beta_{0}=0$) case, these are Adams–Bashforth methods with

where $\boldsymbol{f}(\boldsymbol{\mu}):=(\boldsymbol{f}(\boldsymbol{x}^{0},{t^{0}};\boldsymbol{\mu}),\ldots\boldsymbol{f}(\boldsymbol{x}^{{N_{T}}},{t^{{N_{T}}}};\boldsymbol{\mu}){)}$ and the polynomial approximation (in time) of any time-grid-dependent quantity $\boldsymbol{\xi}:=(\boldsymbol{\xi}^{1},\ldots,\boldsymbol{\xi}^{k})$ using data at $(t^{n},\ldots,t^{n+1-k})$ (with $k\geq 1$) is

In the implicit case (with $\beta_{0}\neq 0$), these are Adams–Moulton methods with coefficients $\beta_{j}$ satisfying

Thus, the time-discrete residual (2.13) becomes

where $I=I^{n-1}_{k}$ in the explicit case and $I=I^{n}_{k+1}$ in the implicit case. Substituting the definitions of the time-discrete state (2.11) and velocity (2.4) in (2.19) yields

Reduced-order models

Applying a linear multistep scheme to integrate Eq. (3.2) in time yields the Galerkin O$\Delta$E

2 LSPG projection

In contrast, LSPG projection associates with a minimum-residual formulation applied to the (time-discrete) O$\Delta$E (2.12), i.e.,

The necessary optimality conditions for problem (3.8) associate with stationarity of the objective function, i.e., the solution $\hat{\boldsymbol{x}}^{n}$ satisfies

Eq. (3.10) reveals that LSPG projection adds the term $\frac{\beta_{0}\Delta t}{\alpha_{0}}\frac{\partial\boldsymbol{f}}{\partial{\boldsymbol{\xi}}}(\boldsymbol{x}^{0}(\boldsymbol{\mu})+\boldsymbol{\Phi}\hat{\boldsymbol{w}},{t^{n}};\boldsymbol{\mu})\boldsymbol{\Phi}$ to the test basis employed by Galerkin projection.

3 Hyper-reduction

In the case of nonlinear dynamical systems, projection is insufficient to yield computational savings, as high-dimensional nonlinear quantities $\boldsymbol{r}$ and $\boldsymbol{r}^{n}$ must be repeatedly computed, projected as $\boldsymbol{\Phi}^{T}\boldsymbol{r}$ and $({\boldsymbol{\Psi}}^{n})^{T}\boldsymbol{r}^{n}$, and differentiated (in the case of implicit time integrators) for Galerkin and LSPG ROMs, respectively. To reduce this computational bottleneck, several ‘hyper-reduction’ techniques have been developed that require computing only a sample of the elements of these nonlinear vector-valued functions. These techniques include collocation , gappy POD , the empirical interpolation method (EIM) , reduced-order quadrature , finite-element subassembly methods , and reduced-basis-sparsification techniques .

respectively. These residual approximations are typically constructed in one of two ways. Later, Section 4.5 proposes a third technique tailored to finite-volume discretizations.

Residual hyper-reduction. This approach amounts to

in the case of gappy POD hyper-reduction, or simply

where $\boldsymbol{A}=(\boldsymbol{P}_{r}\boldsymbol{\Phi}_{r})^{+}\boldsymbol{P}_{r}$ and $\boldsymbol{A}=\boldsymbol{P}_{r}$ in the case of gappy POD and collocation, respectively.

Velocity hyper-reduction. This approach employs an approximated residual constructed from hyper-reduction performed on the velocity vector only, i.e.,

We note that the gappy POD approximations are equivalent to empirical interpolation when the number of samples is equal to the number of reduced-basis elements (i.e., $n_{p,r}=p_{r}$, $n_{p,f}=p_{f}$), as the pseudo-inverse is equal to the inverse and the approximation interpolates the nonlinear function at the sampled elements in this case. Further, the POD–(D)EIM method corresponds to Galerkin projection with gappy POD velocity hyper-reduction and $n_{p,f}=p_{f}$, in which case the hyper-reduced Galerkin ODE becomes

In addition, the GNAT method corresponds to LSPG projection with gappy POD residual hyper-reduction. In principle, the two projection techniques and two hyper-reduction approaches above yield four possible (hyper-reduced) reduced-order models that could be constructed.

4 Lack of conservation

Remarks 4 and 5 demonstrated that Galerkin and LSPG ROMs minimize the violation of conservation in the case of finite-volume models in particular senses; Galerkin performs this minimization at the time-continuous level, while LSPG does so at the time-discrete level. While this is an attractive property, it does not guarantee that the model is conservative in any sense: because the minimum value of the objective functions in Eqs. (3.4) and (3.7) may be non-zero, conservation is generally violated by each of these approaches. We interpret this as violating the structure intrinsic to finite-volume models. This provides the motivation for this work: we aim to develop reduced-order models that ensure the resulting model is conservative globally and—more generally—over subdomains.

Proposed method

This section describes the proposed method, which equips the optimization problems characterizing the online ROM solution with equality constraints that explicitly enforce conservation over subdomains. The approach requires no modification to the offline stage except when hyper-reduction is applied to the nonlinear terms appearing in the constraints. Section 4.1 introduces the concept of conservation over subdomains, Section 4.2 introduces conservative Galerkin projection, Section 4.3 describes conservative LSPG projection, Section 4.4 described approaches for handling infeasibility, and Section 4.5 describes hyper-reduction techniques applicable to objective function and constraint, and Section 4.6 describes snapshot-based (offline) training procedures that may be used for generating the reduced-basis matrices required by the method.

Enforcing conservation (2.1) on each subdomain in the decomposed mesh yields

Critically, noting that $\bar{\boldsymbol{B}}={\bar{\boldsymbol{C}}}\boldsymbol{B}$ due to the fact that neighboring control volumes have outward unit normals of opposite sign along a shared face, we have

Thus, conservation on the decomposed mesh $\bar{\mathcal{M}}$ given an underlying finite-volume discretization on mesh $\mathcal{M}$ can be expressed as

Applying a linear multistep scheme to discretize (4.12) in time yields

Note that the decomposed ODE (4.12) and decomposed O$\Delta$E (4.14) are underdetermined, as they comprise $\bar{N}$ equations in $N(\geq\bar{N})$ unknowns.

We now demonstrate that conservation that is enforced over a decomposed mesh automatically leads to conservation over a coarser mesh that embeds the decomposed mesh.

and satisfaction of time-discrete conservation on $\bar{\mathcal{M}}$ (i.e., Eq. (4.14)) implies satisfaction of time-discrete conservation on $\bar{\bar{\mathcal{M}}}$, i.e.,

Thus, Eqs. (4.15) and (4.16) can be rewritten as

which are clearly satisfied if Eqs. (4.12) and (4.14) are satisfied, respectively. $\square$

The full-order model satisfies time-continuous and time-discrete conservation over any decomposed mesh.

Proof

Proof

We now derive the proposed conservative Galerkin and conservative LSPG projection techniques, which equip their associated optimization problems with equality constraints that enforce conservation over the decomposed mesh $\bar{\mathcal{M}}$.

2 Conservative Galerkin projection

Equivalently, the conservative Galerkin generalized coordinates $\frac{d\hat{\boldsymbol{x}}}{dt}\left(\boldsymbol{x}^{0}(\boldsymbol{\mu})+\boldsymbol{\Phi}\hat{\boldsymbol{x}},t;\boldsymbol{\mu}\right)$ satisfy

We now provide a finite-volume interpretation of the conservative Galerkin model, define the feasible set, and provide an algebraic description of the solution.

If Problem (4.21) is feasible, then the solution is unique and satisfies the time-dependent saddle-point problem

The Lagrangian associated with problem (4.21) can be written as

We note that problem (4.21) corresponds to a convex linear least-squares problem with linear equality constraints; thus, $(\frac{d\hat{\boldsymbol{x}}}{dt},\frac{d\boldsymbol{\lambda}_{\text{G}}}{dt})$ is a unique solution if and only if it satisfies the stationarity conditions

Now, substituting Eqs. (4.30) with $\hat{\boldsymbol{v}}_{1}$ defined in (4.32) into Problem (4.21) yields an unconstrained optimization problem in $\hat{\boldsymbol{v}}_{2}$ only, i.e., $\left[\frac{d\hat{\boldsymbol{x}}}{dt}\right]_{2}$ is the solution to

Applying Eqs.(4.34), and (4.36) to Eq. (4.33) yields

We now show that the conservative Galerkin velocity can be expressed as the orthogonal projection of the standard Galerkin velocity onto the feasible set.

If Problem (4.20) is feasible, then the solution corresponds to the orthogonal projection of the standard Galerkin velocity (3.2) onto the feasible set, i.e.,

Proof

We first identify the feasible set from Eqs. (4.33) and (4.34) as

Of course, numerically solving the conservative Galerkin ROM ODE, requires introducing a time integrator. Applying a linear multistep scheme to solve Eq. (4.23) characterizing the conservative Galerkin ROM ODE yields at time instance $n$ yields the conservative Galerkin ROM O$\Delta$E

We now demonstrate that conservative Galerkin projection and time discretization are commutative.

Further, performing conservative Galerkin projection on Eq. (4.41) and subsequently applying time discretization yields the same model as first applying time discretization on Eq. (4.41) and subsequently performing conservative Galerkin projection.

Proof

The first part of the theorem can be derived by noticing that substituting Eq. (3.1) in (4.41) and premultiplying by $\begin{bmatrix}\boldsymbol{\Phi}&\boldsymbol{0}\\ \boldsymbol{0}&\boldsymbol{I}\end{bmatrix}$ yields the conservative Galerkin saddle-point system (4.23). Then, applying a linear multistep scheme to solve Eq. (4.23) yields the conservative Galerkin ROM O$\Delta$E (4.40) above. Now, applying a linear multistep scheme to integrate (4.41) in time yields

Because applying conservative Galerkin projection to Eq. (4.42) yields Eq. (4.40), we conclude that conservative Galerkin projection and time discretization are commutative. $\square$

3 Conservative LSPG projection

Equivalently, the conservative LSPG generalized coordinates $\hat{\boldsymbol{x}}^{n}$ satisfy

We now provide a finite-volume interpretation of the proposed model, define the feasible set, and provide an algebraic description of the solution.

If Problem (4.44) is feasible, then a solution exists and satisfies the nonlinear saddle-point problem

Defining the Lagrangian associated with problem (4.43) as

the solution $(\hat{\boldsymbol{x}}^{n},\boldsymbol{\lambda}_{\text{P}}^{n})$ satisfies the first-order necessary optimality conditions associated with problem (4.43), i.e., $\partial\mathcal{L}_{L}^{n}/\partial\hat{\boldsymbol{z}}(\hat{\boldsymbol{x}}^{n},\boldsymbol{\lambda}_{\text{P}}^{n};\boldsymbol{\mu})=\boldsymbol{0}$ and ${\partial\mathcal{L}_{L}^{n}}/{\partial\boldsymbol{\gamma}}(\hat{\boldsymbol{x}}^{n},\boldsymbol{\lambda}_{\text{P}}^{n};\boldsymbol{\mu})=\boldsymbol{0}$, which—using the definition of the test basis in Eq. (3.10)—are equivalent to Eqs. (4.46). $\square$

Any appropriate optimization algorithm could be applied to solve minimization problem (4.43) characterizing the conservative LSPG ROM at each time instance. In this work, we propose solving problem (4.43) using the sequential quadratic programming (SQP) method with the Gauss–Newton Hessian approximation. This amounts to applying Newton’s method (with globalization) to the first-order necessary optimality conditions (4.46) and neglecting the term involving differentiation of the test basis $\boldsymbol{\Psi}^{n}(\hat{\boldsymbol{x}}^{n};\boldsymbol{\mu})$. After choosing an initial guess $\hat{\boldsymbol{x}}^{n(0)}$, this approach leads to the following iterations for $k=0,\ldots,K$

4 Handling infeasibility

Of course, the optimization problems characterizing conservative Galerkin projection (i.e., problems (4.20)–(4.21)) and conservative LSPG projection (i.e., problems (4.43)–(4.44)) may not be feasible for arbitrary decomposed meshes $\bar{\mathcal{M}}$ and reduced basis matrices $\boldsymbol{\Phi}$. For example, if the decomposed mesh corresponds to the original mesh (i.e., $\bar{\mathcal{M}}=\mathcal{M}$) and the reduced basis is low-dimensional (i.e., $p\ll N$), then the constraints in these problems correspond to exactly satisfying the full-order-model equations over a low-dimensional subspace; it is likely impossible to do so.

Coarsen the decomposed mesh. First, the number of constraints can be reduced by coarsening the decomposed mesh, i.e., replace $\bar{\mathcal{M}}$ by another decomposed mesh characterized by fewer subdomains ${N_{\bar{\Omega}}}$. As this reduces the number of constraints, the likelihood of feasibility increases, although feasibility remains not guaranteed. This procedure can be repeated until the decomposed mesh leads to a nonempty feasible set or a decomposed mesh characterized by one subdomain (${N_{\bar{\Omega}}}=1$) is infeasible.

If a decomposed mesh leading to feasibility is constructed via coarsening, the conservative reduced-order model can be redefined using the new decomposed mesh and the reduced-order-model simulation can be either (1) reinitialized and restarted from $t=0$, or (2) resumed from the time instance $t^{n}$ where infeasibility was detected. The first approach facilitates analysis, as the reduced-order-model trajectory association with a fixed decomposed mesh, while the latter precludes the need to re-simulate any part of the time interval. Further, if the new decomposed mesh is a decomposition of the previous decomposed mesh, and the previous decomposed mesh is non-overlapping, then conservation over the new decomposed mesh holds over the first part of the time interval (see Theorem 4.1). We note that this approach is not guaranteed to ensure feasibility, as it is possible for infeasibility to exist even in the case of ${N_{\bar{\Omega}}}=1$.

Penalty formulation. Alternatively, infeasibility can be addressed by including the constraints in the objective function via penalization. In the case of conservative Galerkin projection, problem (4.21) is reformulated as

while the conservative LSPG projection problem (4.44) is reformulated as

5 Hyper-reduction

To enable hyper-reduction for the proposed conservative reduced-order models, in addition to approximating the nonlinear objective functions that appear in optimization problems (4.21) and (4.44) as previously described in Section 3.3, we must also approximate the nonlinear constraints ${\bar{\boldsymbol{C}}}\boldsymbol{r}(\boldsymbol{\Phi}\hat{\boldsymbol{v}},\boldsymbol{x}^{0}+\boldsymbol{\Phi}\hat{\boldsymbol{x}},t;\boldsymbol{\mu})=\boldsymbol{0}$ and ${\bar{\boldsymbol{C}}}{\boldsymbol{r}}^{n}(\boldsymbol{x}^{0}(\boldsymbol{\mu})+\boldsymbol{\Phi}\hat{\boldsymbol{z}};\boldsymbol{\mu})=\boldsymbol{0}$. To accomplish this, we propose applying hyper-reduction to the nonlinear residuals that appears in the constraints, i.e., the constraints become

In addition to the two forms of hyper-reduction introduced in Section 3.3, we also propose a third type that leverages the underlying finite-volume discretization of the governing equations:

Flux and source hyper-reduction. This approach respects the underlying decomposition of the velocity vector. It adopts the same residual approximation (3.16)–(3.17) as velocity hyper-reduction (approach 2 in Section 3.3), but employs separate approximations for each term comprising the velocity, i.e.,

One can consider a hierarchy of models that employ objective functions and constraints, each of which may or may not employ one of the three proposed hyper-reduction techniques. For this purpose, we define the Tier-1 and Tier-2 Galerkin and LSPG objective functions as

and the Tier-0 (unconstrained), Tier-1, and Tier-2 Galerkin and LSPG constraints as

Then, we say the Tier A-B Galerkin ROM solution $\frac{d\hat{\boldsymbol{x}}}{dt}$ is the solution to

and the Tier A-B LSPG ROM solution $\hat{\boldsymbol{x}}^{n}$ is the solution to

Note that Tier $i$-0 models correspond to the (original) unconstrained models, Tier $i$-1 models enforce conservation over subdomains, and Tier $i$-2 models enforce approximate conservation over subdomains. The penalty-method variants of the Tier A-B Galerkin and LSPG ROMs are, respectively,

We note that the computational cost incurred by evaluating the constraints is often significantly lower than the cost of evaluating the objective function. For example, for a linear or zero source term, the only nonlinear contribution to the constraints arises from the face flux along the boundary of the subdomains comprising the decomposed mesh. For a small number of subdomains (e.g., global conservation with ${N_{\bar{\Omega}}}=1$), this requires computing only a small number of the elements of the face-flux vector $\boldsymbol{h}$, even without hyper-reduction. Thus, applying hyper-reduction to the objective function is generally more important for computational-cost reduction than applying hyper-reduction to the constraints, i.e., Tier 2–1 ROMs may be preferable to Tier 2–2 ROMs, as their cost is often similar and the former strictly enforces conservation.

6 Snapshot-based training

Here, we propose to construct the reduced-basis matrices $\boldsymbol{\Phi}$, $\boldsymbol{\Phi}_{r}$, $\boldsymbol{\Phi}_{f}$, $\boldsymbol{\Phi}_{h}$, and $\boldsymbol{\Phi}_{s}$ during the offline stage using proper orthogonal decomposition (POD). In particular, given a set of training parameter instances $\mathcal{D}_{\text{train}}:=\{\boldsymbol{\mu}^{1}_{\text{train}},\ldots,\boldsymbol{\mu}^{n_{\text{train}}}_{\text{train}}\}\subset\mathcal{D}$, we execute training simulations from which we compute ‘data tensors’

The reduced-basis matrix associated with each data tensor can be computed as the dominant left singular vectors of its mode-1 unfolding; for example, the state basis $\boldsymbol{\Phi}\equiv[\boldsymbol{\phi}_{1}\ \cdots\ \boldsymbol{\phi}_{p}]$ is computed as

where $\boldsymbol{X}(\boldsymbol{\mu}):=[\boldsymbol{x}^{1}(\boldsymbol{\mu})\ \cdots\ \boldsymbol{x}^{{N_{T}}}(\boldsymbol{\mu})]$ is often referred to as the ‘snapshot matrix’.

Further, we propose to construct the sampling matrices $\boldsymbol{P}_{r}$, $\boldsymbol{P}_{f}$, $\boldsymbol{P}_{h}$, and $\boldsymbol{P}_{s}$ using the sample-mesh greedy method presented in Ref. , which allows for oversampling to enable least-squares regression via gappy POD and also constructs a ‘sample mesh’ wherein all residual elements associated with a given control volume are sampled. However, rather than constructing each of these sampling matrices independently, we propose to construct $\boldsymbol{P}_{r}$ according to the greedy method executed with basis $\boldsymbol{\Phi}_{r}$ and subsequently set $\boldsymbol{P}_{f}=\boldsymbol{P}_{s}=\boldsymbol{P}_{r}$. Further, we construct $\boldsymbol{P}_{h}$ to select the faces associated with the control volumes sampled by $\boldsymbol{P}_{r}$; this corresponds to selecting the sampling matrix $\boldsymbol{P}_{h}$ with the maximum number of rows such that $\boldsymbol{P}_{r}\boldsymbol{B}=\boldsymbol{P}_{r}\boldsymbol{B}\boldsymbol{P}_{h}^{T}\boldsymbol{P}_{h}$.

Analysis

This section performs analysis of the proposed conservative Galerkin and conservative LSPG techniques. For simplicity, we focus on the models without hyper-reduction; the hyper-reduced variants of the results can be derived in a similar manner by making the obvious substitutions.

We first derive sufficient conditions under which the optimization problems characterizing conservative Galerkin and conservative LSPG projection are feasible.

Proof

Case 1. If an explicit scheme is employed, then $\beta_{0}=0$ and the feasible set becomes

2 Equivalence conditions

We now derive conditions under which conservative Galerkin and conservative LSPG projection are equivalent.

The discrete-time conservative Galerkin ROM solution is equivalent to the conservative LSPG solution if either (1) an explicit scheme is employed or (2) the limit $\Delta t\rightarrow 0$ is taken. Further, under these conditions, the Lagrange multipliers are related as

If either $\beta_{0}=0$ or the limit $\Delta t\rightarrow 0$ is taken, the second term vanishes such that we have

This expression is equivalent to $a\sum_{j=0}^{k}\alpha_{j}\boldsymbol{\Phi}^{T}{\bar{\boldsymbol{C}}}^{T}\boldsymbol{\lambda}_{\text{G}}^{n-j}$ with $a=\alpha_{0}$ if Eq. (5.4) holds. $\square$

3 Error analysis

We now derive several a posteriori error bounds for (components of) the solution computed by the proposed conservative model-reduction methods. We employ some of the same techniques used for error analysis in Ref. . For notational simplicity, we drop dependence of the operators on the parameters $\boldsymbol{\mu}$.

We begin by writing the discrete equations characterizing the full-order, Galerkin, and LSPG models as

We also assume Lipschitz continuity of $\boldsymbol{f}$ in its first argument:

To simplify notation, we define the Galerkin and LSPG operators as

respectively, and the Galerkin and LSPG state-space errors at time instance $n$ as

respectively. Because the time instance of the first and second arguments of $\boldsymbol{f}$ always match for linear multistep schemes, we omit the second argument (time) from $\boldsymbol{f}$ in the remainder of this section. All norms in this section correspond to the Euclidean norm, i.e., $\|\cdot\|=\|\cdot\|_{2}$.

Note that ${\bar{\boldsymbol{C}}}^{+}{\bar{\boldsymbol{C}}}={\boldsymbol{V}}_{\bar{\boldsymbol{C}}}{\boldsymbol{V}}_{\bar{\boldsymbol{C}}}^{T}$.

If ${\bf A_{1}}$ holds and $\Delta t<|\alpha_{0}^{n}|/(|\beta_{0}^{n}|\kappa)$, then

Using $\Delta t<|\alpha^{n}_{0}|/(|\beta_{0}^{n}|\kappa)$, we have

Combining inequalities (5.25) and (5.27) yields the final result (5.20). The Galerkin counterpart (5.19) can be derived by following the same steps with the Galerkin operators. $\square$

If ${\bf A_{1}}$ holds and $\Delta t<|\alpha_{0}^{n}|/(|\beta_{0}^{n}|\kappa)$, then

Proof

If ${\bf A_{1}}$ holds and $\Delta t<|\alpha_{0}^{n}|/(|\beta_{0}^{n}|\kappa)$, then

Proof

Now, using $\Delta t<|\alpha_{0}^{n}|/(|\beta_{0}^{n}|\kappa)$ yields

The error in the conserved quantities can be bounded as

Proof

The result can be obtained trivially by subtracting Eq. (5.11) from the premultiplication of Eq. (5.8) by ${\bar{\boldsymbol{C}}}$ and applying the triangle inequality. $\square$

Proof

Subtracting Eq. (5.11) from the premultiplication of Eq. (5.8) by ${\bar{\boldsymbol{C}}}$ yields

Applying the triangle inequality and following the above steps yields

Combining inequalities (5.53) and (5.55) produces the final result. $\square$

Numerical experiments

This section compares the performance of several reduced-order models on a parameterization of the quasi-1D Euler equations applied to supersonic flow in a converging–diverging nozzle.

We consider a parameterized quasi-1D Euler equation associated with modeling inviscid compressible flow in a one-dimensional converging–diverging nozzle with a continuously varying cross-sectional area [40, Chapter 13]; Figure 2 depicts the problem geometry.

In integral form, the governing equations are:

$\forall\omega\subseteq\Omega=[0,L]$. Thus, the governing system of nonlinear partial differential equations is consistent with the conservation-law formulation in Eq. (2.1) with $d=1$ spatial dimension, $n_{u}=3$ conserved variables corresponding to density $u_{1}={A}\rho$, momentum $u_{2}={A}\rho u$, and energy density $u_{3}={A}e$. The flux corresponds to ${g}_{1}={A}\rho u$, ${g}_{2}={A}{(}\rho u^{2}+p{)}$, and ${g}_{3}={A}(e+p)u$, and the source corresponds to $s_{1}=s_{3}=0$ and $s_{2}=p\frac{\partial A}{\partial x}$. In addition, we have $p=(\gamma-1)\rho\epsilon$, $\epsilon=\frac{e}{\rho}-\frac{u^{2}}{2}$, and assume a perfect gas (i.e., $p=\rho RT$). Here, $\rho$ denotes density, $u$ denotes velocity, $p$ denotes pressure, $\epsilon$ denotes potential energy per unit mass, $e$ denotes total energy density, $\gamma$ denotes the specific heat ratio, and $A$ denotes the converging–diverging nozzle cross-sectional area. We employ a specific heat ratio of $\gamma=1.3$ and a specific gas constant of $R=355.4$ $\text{m}^{2}/\text{s}^{2}/\text{K}$. The spatial domain is $\Omega=[0,L]$ with $L=0.25$ m. The cross-sectional area $A(x)$ is determined by a cubic spline interpolation over the points

The initial flow field is created in several steps. First, the following isentropic relations are used to generate a zero pressure-gradient flow field at the inlet ($x=0$ m) and the outlet ($x=0.25$ m):

where $c$ denotes the speed of sound, the total temperature is $T_{t}=2800$ K, and the total pressure is $p_{t}=2.068\times 10^{6}$ $\text{N}/\text{m}^{2}$.

2 Compared methods

These experiments compare the following methods, which employ the Tier A-B notation established in Section 4.5:

FOM. This model corresponds to the full-order model, i.e., the solution satisfying Eq. (2.5).

Galerkin. This model corresponds to the Tier 1-0 Galerkin ROM.

LSPG. This model corresponds to the Tier 1-0 LSPG ROM.

LSPG-FV. This model corresponds to the Tier 1-1 LSPG ROM, which is conservative.

GNAT-FV. This model corresponds to the Tier 2-1 LSPG ROM, which is conservative. While objective function is approximated in the same way as in the GNAT method above, no hyper-reduction is applied to the constraints.

and the mean-squared and time-instantaneous error in the globally conserved variables

All timings are obtained by performing calculations on an Intel(R) Xeon(R) CPU E5-2670 @ 2.60 GHz, 31.4 GB RAM using the MORTestbed in Matlab.

3 GNAT-FV(X) snapshot study

This section assesses the effect of snapshot-collection method on the performance of the GNAT-FV(X) method; all subsequent experiments employ the snapshot-collection method yielding the best performance.

We set the number of control volumes to ${N_{\Omega}}=100$ such that $N={N_{\Omega}}n_{u}=300$, the reduced-basis dimensions to $p=5$ (which corresponds to a relative statistical energy of 99.78%) and $p_{r}=p_{h}$, and employ a sample mesh of 20 control volumes, which corresponds to $n_{p,r}=n_{p,f}=n_{p,h}=60$. We employ a penalty parameter of $\rho=10^{3}$, which is used by infeasibility-handling approach 2 when $\bar{N}>p$. In this setting we vary the number of constraints $\bar{N}$ and flux-basis dimension $p_{h}$ and report the relative mean-squared violation in global conservation over the time interval, i.e., the value of $\mathcal{E}_{\boldsymbol{r},\text{global}}$ for the given reduced-order model divided by the value of $\mathcal{E}_{\boldsymbol{r},\text{global}}$ for the (unconstrained) GNAT model; note that this value is zero if the constraint-approximation error is zero and a feasible solution is computed at each time instance.

Figure 3 reports the results for this experiment and elucidates several trends. First, Figure 3(a) shows that the GNAT-FV model—for which the constraints are enforced exactly—yields near-exact satisfaction of the conservation laws for $\bar{N}<p$; this implies that a feasible solution was computed at every time instance of that simulation.

Second, we note that for $1<\bar{N}/p<2$, the GNAT-FV model yields approximate but accurate satisfaction of the conservation laws, as the relative value of $\mathcal{E}_{\boldsymbol{r},\text{global}}$ is less than $10^{-2}$ in these cases.

Third, Figures (3(b))–(3(f)) show that the best results for the GNAT-FV(X) method are obtained for X=LSPG-FV (Figure 3(e)) and X=GNAT-FV (Figure 3(f)); these techniques yield relative values of $\mathcal{E}_{\boldsymbol{r},\text{global}}$ for the GNAT-FV(X) model less than $10^{-2}$ for $\bar{N}/p<2$ in almost all cases. This result is sensible, as the training simulations corresponding to (constrained) LSPG-FV and GNAT-FV are ‘closer’ to the (constrained) GNAT-FV(X) simulation relative to the (unconstrained) FOM, LSPG, and GNAT simulations. However, in these cases, the relative value of $\mathcal{E}_{\boldsymbol{r},\text{global}}$ for $\bar{N}<p$ is small, but not close to machine zero as in the GNAT-FV case because the constraints are approximated. Thus, these methods—while having a cost independent of $N$ due to the introduction of hyper-reduction—are only approximately conservative.

Fourth, we note that the GNAT-FV(LSPG-FV) and GNAT-FV(LSPG-FV) results are insensitive to the flux-basis dimension $p_{h}$ for $p_{h}$ sufficiently large ($p_{h}>12$).

Finally, we note that while the GNAT-FV(LSPG-FV) and GNAT-FV(GNAT-FV) models yield similar accuracy, the latter method incurs a lower training cost, as the former incurs training simulations with the (Tier 1-1) LSPG-FV model, while the latter incurs training simulations with the (Tier 2-1) GNAT-FV model. Thus, the only GNAT-FV(X) method we consider in subsequent experiments is the GNAT-FV(GNAT-FV) approach.

4 Penalty-parameter study

This section assesses the effect of the penalty parameter $\rho$ employed by infeasibility-handling approach 2 when $\bar{N}>p$ on the performance of the (constrained) ROMs LSPG-FV, GNAT-FV, and GNAT-FV(GNAT-FV). All subsequent experiments employ the penalty parameter yielding the best performance.

We again set the number of control volumes to ${N_{\Omega}}=100$ such that $N={N_{\Omega}}n_{u}=300$, the reduced-basis dimensions to $p=5$ and $p_{r}=p_{h}=p_{s}=20$ and again employ a sample mesh of 20 control volumes, which corresponds to $n_{p,r}=n_{p,f}=n_{p,h}=60$. We vary the number of constraints $\bar{N}$ and penalty parameter $\rho$ and report the mean-squared state-space error ${\mathcal{E}_{\boldsymbol{x}}}$ and the (absolute) mean-squared violation in global conservation over the time interval $\mathcal{E}_{\boldsymbol{r},\text{global}}$. We note that a penalty value of $\rho=\infty$ corresponds to minimizing the norm of the constraints only (i.e., the objective function is ignored).

Figure 4 reports the results for this experiment. First, we note that values $\rho\in\{10,10^{2},10^{3}\}$ yield similar performance, which outperforms the other tested values. In particular, values of $\rho\in\{1,\infty\}$ often yield unstable responses, while $\rho=10$ almost always yields larger errors than employing $\rho\in\{10,10^{2},10^{3}\}$.

Second, we note that nearly all cases outperform the unconstrained model, characterized by $\rho=0$; this implies that employing the proposed constraints can improve accuracy, even if the constraints are employed in a penalty formulation rather than as strictly enforced constraints.

Third, we observe that the two reported metrics are often correlated: larger values of mean-squared violation in global conservation $\mathcal{E}_{\boldsymbol{r},\text{global}}$ typically implies larger values of the relative mean-squared state-space error ${\mathcal{E}_{\boldsymbol{x}}}$. This lends credibility to the proposed technique, which aims to reduce the violation in global conservation, as it suggests that enforcing this constraint (or employing it as a penalty in the objective function) can lead to more accurate ROMs.

In subsequent experiments, we employ a penalty-parameter value of $\rho=10^{3}$.

5 State-basis-dimension study

This section assesses the effect of basis dimension $p$ on the proposed methods. We again employ ${N_{\Omega}}=100$ control volumes in the finite-volume discretization, set reduced-basis dimensions to $p_{r}=p_{h}=p_{s}=20$, employ a sample mesh with 20 control volumes, and set the penalty parameter to $\rho=10^{3}$. We vary both the state-basis dimension $p$ and the number of constraints $\bar{N}$ the relative mean-squared state-space error ${\mathcal{E}_{\boldsymbol{x}}}$ and the (absolute) mean-squared violation in global conservation over the time interval $\mathcal{E}_{\boldsymbol{r},\text{global}}$.

Figure 5 reports the results. First, and most importantly, we note that Figures 5(a), 5(c), and 5(e) show that the introduction of constraints yields the most significant improvements for the smallest basis dimension $p=5$. In these cases, the relative mean-squared state-space error ${\mathcal{E}_{\boldsymbol{x}}}$ is reduced by over an order of magnitude for all ROMs, as the unconstrained ROMs yield errors exceeding 30%, while their constrained counterparts employing $\bar{N}=3$ (i.e., global conservation with $\bar{\mathcal{M}}={\bar{\mathcal{M}}}_{\text{global}}$) all yield errors less than 2%. In contrast, for $p\geq 7$, the unconstrained ROMs are already quite accurate, with errors already less then 2%; incorporating constraints in these cases does yield accuracy improvements in most cases, although these improvements are less dramatic. Because the most significant improvements were obtained by enforcing global conservation with $\bar{N}=3$, subsequent experiments employ ROMs that enforce global conservation by using a decomposed mesh of $\bar{\mathcal{M}}={\bar{\mathcal{M}}}_{\text{global}}$.

Second, Figures 5(b) and 5(d) show that the LSPG-FV and GNAT-FV models produce near-exact satisfaction of the conservation laws for $\bar{N}<p$; this implies that a feasible solution was computed at every time instance of the corresponding simulations. In contrast, Figure 5(f) shows that the GNAT-FV(GNAT-FV) ROM is only approximately conservative. Nonetheless, this approximate conservation does not adversely impact the actual errors produced by the ROM, as the errors reported in Figures 5(c) and 5(e) are nearly identical in all cases. So, while applying hyper-reduction to the constraints results in a loss of numerically exact satisfaction of global conservation, the results are extremely similar to the case where the constraints are applied exactly.

6 Comparison across all methods

This section assesses the relative performance of the methods over time; all ROMs that employ constraints enforce global conservation, i.e., $\bar{N}=3$ and $\bar{\mathcal{M}}={\bar{\mathcal{M}}}_{\text{global}}$.

Figure 6 reports the results. First, we note that the errors $\varepsilon_{\boldsymbol{x}}^{n}$ and $\varepsilon_{\boldsymbol{x},\text{global}}^{n}$ exhibit the same trends in all cases; this suggests that enforcing global conservation—which leads to lower errors in the globally conserved quantities by construction—is an effective approach for also reducing the error in the state itself. This also supports previous observations that enforcing global conservation rather than employing a penalty approach leads to smaller errors in most cases.

Second, we observe that the FOM, LSPG-FV, and GNAT-FV models all lead to global-conservation violations $\varepsilon_{\boldsymbol{x},\text{global}}^{n}$ near zero as expected. In contrast, the GNAT-FV(GNAT-FV) approach only approximately satisfies global conservation due the introduction of hyper-reduction to the constraints; however, this has no noticeable effect on its response, as the errors reported for GNAT-FV and GNAT-FV(GNAT-FV) are nearly identical in Figures 6(a)–6(d).

Third, we notice that the conservative methods LSPG-FV and GNAT-FV, as well as the approximately conservative method GNAT-FV(GNAT-FV), all yield significantly lower errors than the unconstrained methods Galerkin, LSPG, and GNAT. Further, these unconstrained methods yield significant violation in global conservation.

Table 1 reports the timings for these methods. We first note that the LSPG ROM does not have a valid timing for either problem, as the associated simulations yield negative pressures and thus do not successfully run for the entire time interval (see premature termination in Figure 6). Second, while all other ROMs produce a speedup relative to the FOM, methods that employ hyper-reduction for the objective function (GNAT, GNAT-FV) produce more significant speedups; further applying hyper-reduction to the constraint (GNAT-FV(GNAT-FV)) improves the speedup further.

To enable an objective comparison of the ROM methods, we compare their performance across a wide variation of all method parameters. We subject each model to a parameter study wherein each model parameter is varied between the limits specified in Table 2. From these results, we then construct a Pareto front for each method, which is characterized by the method parameters that minimize the competing objectives of error and wall time.

Figure 7 reports these Pareto fronts, where both the mean-squared state-space error ${\mathcal{E}_{\boldsymbol{x}}}$ and mean-squared violation in global conservation $\mathcal{E}_{\boldsymbol{r},\text{global}}$ are considered as error measures, as well as an ‘overall’ Pareto front that selects the Pareto-optimal methods across all parameter variations. Note that this figure reports the relative wall time with respect to that of the FOM simulation; relative wall times less than one imply the ROM yields a speedup. Here, Figure 7(a) shows that the GNAT-FV(GNAT-FV) method is always Pareto dominant for error measure ${\mathcal{E}_{\boldsymbol{x}}}$, as no other method is both less expensive and more accurate for any tested parameter combination. The method that performs second best is the proposed GNAT-FV method, which exactly enforces constraints; note that it is only slightly more expensive than the GNAT-FV(GNAT-FV) method, as the benefit of performing hyper-reduction on the residual appearing in the constraints with $\bar{N}$ small is much less significant than the benefit of performing hyper-reduction on the residual appearing in the objective function when $\bar{N}$ is small. In particular, note that the Pareto-optimal parameter combinations for the LSPG-FV method yield similar accuracy to the Pareto-optimal GNAT-FV(GNAT-FV) points, but incur significantly larger wall times. Figure 7(b) shows that GNAT-FV(GNAT-FV) is Pareto optimal for error measure $\mathcal{E}_{\boldsymbol{r},\text{global}}$ for relative wall times less than 0.28, but GNAT-FV, which enforces constraints exactly, is Pareto optimal for larger relative wall times, yielding near-zero violations in global conservation. We emphasize that both conservative variants of the GNAT method (i.e., GNAT-FV and GNAT-FV(GNAT-FV)) outperform the original GNAT approach, and the conservative variant of the LSPG method (i.e., LSPG-FV) outperforms the original LSPG method; this demonstrates the benefit of the proposed method and the performance improvement gained by enforcing conservation. In particular, note that the introduction of constraints does not adversely affect ROM wall-time performance; in fact, LSPG-FV has better wall-time performance relative to the the LSPG method. This occurs because global conservation corresponds to only $\bar{N}=3$ constraints in this case, and because these constraints lead to improved accuracy and thus promote convergence, the associated simulations require fewer iterations to solve the optimization problem at each time instance. We also note that hyper-reduction is needed to realize significant speedups: Pareto-optimal parameter combinations for ROMs employing hyper-reduction lead to relative wall times less than 0.36, while Pareto-optimal parameter combinations for ROMs without hyper-reduction yield relative wall times exceeding 0.6.

Conclusions

This work proposed two model-reduction methods for finite-volume models that enforce conservation over subdomains: conservative Galerkin and conservative LSPG projection. These methods associate with optimization problems characterized by a minimum-residual objective function and nonlinear equality constraints formulated at the time-continuous and time-discrete levels, respectively. We equipped these methods with techniques for handling infeasible constraints, and we also developed hyper-reduction methods to ensure low-cost ROM simulations in the presence of nonlinear flux or source terms.

We performed analysis that demonstrated commutativity of conservative Galerkin projection and time discretization, developed sufficient conditions for feasibility, demonstrated conditions under which conservative Galerkin and conservative LSPG models are equivalent, and derived a posteriori error bounds. Numerical experiments on a model problem highlighted the benefit of conservative projection, and also demonstrated that enforcing global conservation led to the most accurate results.

Future work involves implementing the proposed techniques in a production-level computational fluid-dynamics code, demonstrating the methods on truly large-scale finite-volume models, and investigating combining the methodology with space–time projection approaches , as these techniques have demonstrated error bounds that grow slowly in time.

Acknowledgments

We thank Matthew Barone and Irina Tezaur for insightful conversations related to structure-preserving model reduction in fluid dynamics. We also thank J. Nathan Kutz for his help in forging the collaboration. This work was funded by Sandia’s Laboratory Directed Research and Development (LDRD) program under Project #190968. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525.

References