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 sirovich1987tad3 ; deane1991low ; ma2002low , wherein the test basis is set to be equal to the trial basis. The trial basis is often computed via proper orthogonal decomposition (POD) POD , but it can also be computed via the reduced-basis method; see Refs. haasdonk2008FV ; haasdonkExplicit ; haasdonk2008reduced , which apply the classical reduced-basis method to finite-volume problems. More recently, the least-squares Petrov–Galerkin (LSPG) projection method CarlbergGappy ; carlbergJCP ; carlbergGalDiscOpt 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 rowley2004mrc and finite-element discretizations barone2009stable ; kalashnikova2010stability ; introducing dissipation via closure models aubry1988dynamics ; sirisup2004spectral ; Bergmann2009516 ; wang2012proper ; san2013proper or numerical dissipation iollo2000stability ; performing nonlinear Galerkin projection based on approximate inertial manifolds marion1989nonlinear ; shen1990long ; jolly1991preserving ; including a pressure-term representation noack2005need ; galletti2004low ; modifying the POD basis by including many modes (such that dissipative modes are captured), changing the norm iollo2000stability , enabling adaptivity Bergmann2009516 ; carlberg2014adaptive , or including basis functions that resolve a range of scales balajewicz2012stabilization or respect the attractor’s power balance balajewicz2013low ; modifying the projection by adopting a constrained Galerkin reddy2017constrained ; fick2017reduced , constrained Petrov–Galerkin Fang2013540 , or -norm minimizing projection abgrall2015robust ; developing approaches tailored to the incompressible Navier–Stokes equations by introducing stabilizations based on supremizer-enriched velocity spaces and a pressure Poisson equation stabile2017finite ; stabile2017advances or by modifying the Galerkin projection lorenzi2016pod ; and improving the ROM’s ability to capture shocks ohlberger2013nonlinear ; gerbeau2014approximated ; carlberg2014adaptive ; taddei2015reduced . 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
Analysis, which includes:
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. zimmermann2014reduced 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. reddy2017constrained 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. fick2017reduced 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.
In this work, matrices are denoted by capitalized bold letters, vectors by lowercase bold letters, and scalars by unbolded letters. The columns of a matrixare denoted by , with such that . The scalar-valued matrix elements are denoted by such that , . A superscript denotes the value of a variable at that time instance, e.g., is the value of at time , where is the time step.
2 Finite-volume discretization
This work considers parameterized systems of conservation laws. In integral form, the associated governing equations correspond to
which is solved in time domain with final time , and a (parameterized) initial condition denoted by such that . Here, with denotes any subset of the spatial domain of interest with , whose boundary is , denotes integration with respect to the boundary, , denotes the th conserved variable (per unit volume); , , denotes the flux associated with the th conserved variable (per unit area per unit time); denotes the outward unit normal to ; , denotes the source associated with the th conserved variable (per unit volume per unit time); and denotes the parameter domain. We assume the domain is independent of the parameters for notational simplicity.
2.1 Spatial discretization
We consider the particular case where the governing equations (2.1) have been discretized in space by a finite-volume method. This implies that the spatial domain has been partitioned into a mesh of non-overlapping (closed, connected) control volumes , such that , which intersect only on their ()-dimensional interface, i.e., for , where , . We define the mesh as , and we denote the boundary of the th control volume by . The th control-volume boundary is partitioned into a set of faces444We note that this is a set of faces for , faces for , or simply extremities for . denoted by such that . Then the full set of faces within the mesh is . Applying Eq. (2.1) to each control volume in the mesh yields
where and denotes the unit normal to control volume . Finite-volume schemes complete the spatial discretization by introducing a state vector with whose elements comprise
denotes a mapping from conservation-law index and control-volume index to degree of freedom, and a velocity vectorwith whose elements consist of
for . Here, , denotes the approximated (or reconstructed) flux associated with the th conserved variable (per unit area per unit time); and , denotes the approximated source associated with the th conserved variable (per unit volume per unit time), which may arise, e.g., from applying a quadrature rule to evaluate the integral. We emphasize that both the approximated flux and approximated source will in general depend on the entire state vector , e.g., due to high-order flux reconstructions or reactions, respectively.
Substituting , , and in Eq. (2.2) and dividing by yields
. This is a parameterized system of nonlinear ordinary differential equations (ODEs) characterizing an initial value problem, which we consider to be our full-order model (FOM).
Remark 1 (Full-order model ODE: finite-volume interpretation)
From the definitions of the state (2.3) and velocity (2.4), the full-order-model ODE residual element can be interpreted as the (normalized) rate of violation of conservation in variable in control volume at time instance under one approximation: the flux and source terms are approximated using the finite-volume discretization (i.e., , and ).
Remark 2 (Flux velocity from face fluxes)
The elements of the flux velocity can be computed from a vector of face fluxes , whose elements are
where denotes the unit normal assigned to face (using any convention) and denotes a mapping from conservation-law index and face index to degrees of freedom defined on the faces. This mapping is provided by
where the elements of are
where 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.
2.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. (carlbergGalDiscOpt, ), which develops LSPG reduced-order models for explicit, fully implicit, and diagonally implicit Runge–Kutta schemes. Applying a linear -step method to numerically solve Eq. (2.5) at a given parameter instance can be written as
where denotes the time step, denotes the numerical approximation to , i.e.,
where denotes the numerical approximation to . The coefficients and define a particular linear multistep scheme, and is necessary for consistency, and the method is implicit if . For notational simplicity, we employ a uniform time grid , with and . The fully discrete full-order model, which is sometimes denoted as the FOM OE, is characterized by the following system of algebraic equations to be solved at each time instance :
where denotes the linear multistep residual, which is defined as
The unknown vector can be interpreted as
and denotes an approximation to the th conserved variable when evaluating the residual (2.13) at the th time instance.
Adams methods. Adams methods consider the integrated form of Eq. (2.5)
and apply a polynomial approximation to the integrand. In particular, the th-order Adams scheme employs coefficients , , and , and coefficients
that associate with a polynomial interpolation of the integrand. In the explicit () case, these are Adams–Bashforth methods with
where and the polynomial approximation (in time) of any time-grid-dependent quantity using data at (with ) is
In the implicit case (with ), these are Adams–Moulton methods with coefficients satisfying
Thus, the time-discrete residual (2.13) becomes
Remark 3 (Full-order model OE: finite-volume interpretation for Adams methods)
Eq. (2.20) shows that the full-order-model OE residual element in the case of Adams methods can be interpreted as the (normalized) violation of conservation in variable in control volume over time interval under two approximations: (1) the flux and source terms are approximated using the finite-volume discretization (i.e., , and ), and (2) a polynomial interpolation is used to approximate the integrand for time integration.
3 Reduced-order models
During the online stage, projection-based reduced-order models compute an approximate solution that lies in a low-dimensional affine trial subspace , i.e.,
where is the reduced-basis matrix of dimension , which we assume without loss of generality satisfies , denotes the generalized coordinates, and denotes the range of a matrix . This basis can be computed in a variety of ways during the offline stage, e.g., eigenmode analysis, POD POD , or the reduced-basis method prud2002reliable ; rozza2007reduced . Substituting the approximation into governing equations (2.5) yields an overdetermined system of equations in unknowns. To compute a unique solution, reduced-order models must enforce the residual to be orthogonal to a -dimensional test subspace. Galerkin and LSPG projection differ in their choices of this subspace; each choice leads to an approximate solution that exhibits a particular notion of optimality.
3.1 Galerkin projection
Galerkin projection employs a test subspace of and thus enforces the residual to be orthogonal to , i.e., the Galerkin ODE is
Applying a linear multistep scheme to integrate Eq. (3.2) in time yields the Galerkin OE
As demonstrated, e.g., in Ref. carlbergGalDiscOpt , Galerkin projection exhibits continuous optimality if the reduced basis is orthogonal, i.e., , as the Galerkin ROM computes the approximated velocity that minimizes the -norm of the FOM ODE residual (2.5) over , i.e.,
denotes the FOM ODE residual.
Remark 4 (Galerkin ROM ODE: finite-volume interpretation)
From the time-continuous optimality of the Galerkin ROM ODE (3.5) and the finite-volume interpretation of the FOM ODE in Remark 1, the Galerkin ROM ODE (3.2) can be interpreted as minimizing the sum of squared (normalized) rates of violation of conservation across all variables , and control volumes , at time instance under one approximation: the flux and source terms are approximated using the finite-volume discretization (i.e., , and ).
3.2 LSPG projection
In contrast, LSPG projection associates with a minimum-residual formulation applied to the (time-discrete) OE (2.12), i.e.,
The necessary optimality conditions for problem (3.8) associate with stationarity of the objective function, i.e., the solution satisfies
where the LSPG test basis is
Eq. (3.10) reveals that LSPG projection adds the term to the test basis employed by Galerkin projection.
Remark 5 (Lspg Rom OE: finite-volume interpretation for Adams methods)
From the time-discrete optimality of the LSPG ROM OE (3.8) and the finite-volume interpretation of the FOM OE for Adams methods in Remark 3, the LSPG ROM OE (3.8) can be interpreted as minimizing the sum of squared (normalized) violation of conservation across all variables , and control volumes , over time interval under two approximations: (1) the flux and source terms are approximated using the finite-volume discretization (i.e., , and ), and (2) a polynomial interpolation is used to approximate the integrand for time integration.
In the case of nonlinear dynamical systems, projection is insufficient to yield computational savings, as high-dimensional nonlinear quantities and must be repeatedly computed, projected as and , 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 astrid2007mpe ; ryckelynck2005phm ; LeGresleyThesis , gappy POD sirovichOrigGappy ; bos2004als ; astrid2007mpe ; CarlbergGappy ; carlbergJCP , the empirical interpolation method (EIM) barrault2004eim ; chaturantabut2010journal ; galbally2009non ; drohmannEOI ; HAntil_MHeinkenschloss_DCSorensen_2013a , reduced-order quadrature HAntil_SField_RHNochetto_MTiglio_2013 , finite-element subassembly methods an2008optimizing ; farhat2014dimensional , and reduced-basis-sparsification techniques carlberg2012spd .
In the present context, hyper-reduction can be achieved by replacing the residuals appearing in the objective functions of (3.5) and (3.8) by and , respectively, such that the hyper-reduced optimization problems for Galerkin and LSPG projection become
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
in the case of collocation. Here, denotes a sampling matrix comprising selected rows of the identity matrix, while denotes a -dimensional reduced-basis matrix constructed for the residual, a superscript denotes the Moore–Penrose pseudoinverse, and denotes the set of full-column rank matrices (the non-compact Stiefel manifold). This approach has the advantage of associating hyper-reduced optimization problems (3.11) and (3.12) with a weighted-norm variant of the original optimization problems (3.5) and(3.8), i.e.,
where and 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.,
in the case of gappy POD or collocation, respectively. Here, denotes a sampling matrix comprising selected rows of the identity matrix, while denotes a -dimensional reduced-basis matrix constructed for the velocity. This approach has the advantage of limiting the hyper-reduction approximation to the nonlinear component of the residual.
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., , ), 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 chaturantabut2010journal corresponds to Galerkin projection with gappy POD velocity hyper-reduction and , in which case the hyper-reduced Galerkin ODE becomes
In addition, the GNAT method CarlbergGappy ; carlbergJCP 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.
3.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.
4 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.
4.1 Domain decomposition
To begin, we decompose the mesh into subdomains, each of which comprises the union of control volumes. That is, we define a decomposed mesh of subdomains , with . We note that the subdomains need not be non-overlapping, closed, or connected. Denoting the boundary of the th subdomain by , we have , with representing the set of faces belonging to the th subdomain. We denote the full set of faces within the decomposed mesh by . Figure 1 depicts several decompositions that satisfy the above conditions. We emphasize that the subdomains can overlap, their union need not correspond to the global domain, and the global domain can be considered by employing , which is characterized by subdomain that corresponds to the global domain, i.e., and , as depicted in Figure 0(c).
Enforcing conservation (2.1) on each subdomain in the decomposed mesh yields
where denotes the unit normal to subdomain . We propose applying a finite-volume discretization to Eq. (4.1) that operates on the decomposed mesh . That is, we introduce a ‘decomposed’ state vector with and elements
where denotes a mapping from conservation-law index and subdomain index to decomposed degree of freedom. The decomposed state vector can be computed from the state vector as
where has elements , where is the indicator function, which evaluates to one if its argument is true, and zero if its argument is false. We note that this matrix can be decomposed as , where the elements of the volumetic matrices , and aggregation matrix comprise