Galerkin v. least-squares Petrov--Galerkin projection in nonlinear model reduction

Introduction

While modeling and simulation of parameterized systems has become an essential tool in many industries, the computational cost of executing high-fidelity simulations is infeasibly high for many time-critical applications. For example, real-time scenarios (e.g., model predictive control) require simulations to execute in seconds or minutes, while many-query scenarios (e.g., statistical inversion) can require thousands of simulations corresponding to different parameter instances of the system.

Reduced-order models (ROMs) have been developed to mitigate this computational bottleneck. First, they execute an offline stage during which computationally expensive training tasks (e.g., evaluating the high-fidelity model at several points in the parameter space) compute a representative low-dimensional ‘trial’ basis for the system state. Then, during the inexpensive online stage, these methods quickly compute approximate solutions for arbitrary points in the parameter space via projection: they compute solutions in the span of the trial basis while enforcing the high-fidelity-model residual to be orthogonal to a low-dimensional ‘test’ basis. They also introduce other approximations in the presence of general nonlinearities or non-affine parameter dependence. See Ref. and references within for a survey of current methods.

To address these shortcomings, alternative projection techniques have been developed, particularly in fluid dynamics. These include stabilizing inner products ; 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 ; and performing Petrov–Galerkin projection .

In spite of its promise, theoretical analysis has been limited to developing consistency conditions for snapshot collection and discrete-time error bounds for simple time integrators . In particular, major outstanding questions include: (1) What are time-continuous and time-discrete representations of the Galerkin and LSPG ROMs for broad classes of time integrators? (2) Are there conditions under which the two techniques are equivalent? (3) What discrete-time error bounds are available for the two techniques for broad classes time integrators? Related to the third issue is how parameters (e.g., time step or basis dimension) for the LSPG ROM should be chosen to optimize performance. This work aims to fill this gap by performing a number of detailed theoretical and computational studies that compare Galerkin and LSPG ROMs for the two most important classes of time integrators: linear multistep methods and Runge–Kutta schemes. We summarize the most important theoretical results (which map to the three questions above) as follows:

Galerkin projection and time discretization are commutative (Theorem 3.4).

LSPG ROMs can be derived for Runge–Kutta schemes (Section 4.1).

The LSPG ROM has a time-continuous (i.e., ODE) representation under certain conditions (Section 4.2, Figure 1). This ODE depends on the time step $\Delta t$.

Galerkin and LSPG ROMs are equivalent for explicit time integrators (Corollaries 5.1 and 5.2).

Galerkin and LSPG ROMs are equivalent in the limit of $\Delta t\rightarrow 0$ for implicit time integrators (Theorem 5.3).

Galerkin ROMs are discrete optimal for symmetric-positive-definite residual Jacobians (Theorems 5.4–5.6).

We provide local and global a posteriori and a priori error bounds for both Galerkin and LSPG ROMs for linear multistep schemes (Section 6.1).

For backward differentiation formulas, we show that the LSPG ROM can yield a lower local a posteriori error bound than the Galerkin ROM (Corollary 6.4) because it solves a time-global optimization problem (over a time window) rather than a time-local optimization problem (Remark 6.9).

For the backward Euler time integrator, we show that an intermediate time step should yield the lowest error bound (Corollary 6.17 and Remark 6.18).

For the backward Euler time integrator, we show that a larger basis size leads to a smaller optimal time step for the LSPG ROM (Corollary 6.17).

We provide global a posteriori and a priori error bounds for both Galerkin and LSPG ROMs for Runge–Kutta schemes (Section 6.3).

Figure 1 summarizes time-continuous and time-discrete representations of the two techniques.

In addition to the above theoretical results, we present numerical results for a large-scale compressible fluid-dynamics problem with turbulence modeling characterized by over one million degrees of freedom. These results illustrate the practical significance of the above theoretical results. Critically, we show that employing an intermediate time step for the LSPG ROM can decrease both the error and the simulation time by an order of magnitude, which is a highly non-intuitive result that is of immense practical significance.

The remainder of the paper is organized as follows. Section 2 formulates the full-order model, including its representation at the time-continuous and time-discrete levels. Section 3 presents the Galerkin ROM at the continuous and discrete levels, and Section 4 does so for the LSPG ROM. In particular, we provide conditions under which the LSPG ROM has a time-continuous representation. Section 5 provides conditions under which the Galerkin and LSPG ROMs are equivalent; in particular, equivalence holds for explicit integrators (Section 5.1), in the limit of the time step going to zero for implicit integrators (Section 5.2), and for symmetric-positive-definite residual Jacobians (Section 5.3). Section 6 provides error analysis for Galerkin and LSPG ROMs for linear multistep schemes (Section 6.1), Runge–Kutta schemes (Section 6.3), and a detailed analysis in the case of backward Euler (Section 6.2). Section 7 provides detailed numerical examples that illustrate the practical importance of the analysis. Finally, Section 8 provides conclusions.

Full-order model

We begin by formulating both the time-continuous (ODE) and time-discrete (O$\Delta$E) representations of the full-order model (FOM).

In this work, we consider the full-order model (FOM) to be an initial value problem characterized by a system of nonlinear ODEs

2 Discrete representation

A time-discretization method is required to solve Eq. (2.1) numerically. We now characterize the full-order-model O$\Delta$E, which is the time-discrete representation of the model, for two classes of time integrators: linear multistep schemes and Runge–Kutta schemes.

A linear $k$-step method applied to numerically solve Eq. (2.1) can be expressed as

Then, the state can be updated explicitly as

Hence, the unknown is equal to the state. These methods are implicit if $\beta_{0}\neq 0$.

2.2 Runge–Kutta schemes

For an $s$-stage Runge–Kutta scheme, the O$\Delta$E is characterized by the following system of algebraic equations to be solved at each time step:

Here, the Runge–Kutta residual is defined as

Galerkin ROM

This section provides the time-continuous and time-discrete representations of the Galerkin ROM, as well as key results related to optimality and commutativity of projection and time discretization.

In order to obtain computational efficiency, it is necessary to reduce the computational complexity of repeatedly computing matrix–vector products of the form $\bm{{\Phi}}^{T}\bm{{f}}$. This can be done using a variety of methods, e.g., collocation , gappy POD , the discrete empirical interpolation method (DEIM) , reduced-order quadrature , finite-element subassembly methods , or reduced-basis-sparsification techniques . However, in this work we limit ourselves to comparatively analyzing different projection techniques. For this reason, we do not perform additional analysis for such complexity-reduction mechanisms; this is the subject of follow-on work.

We now restate the well-known result that Galerkin projection leads to a notion of minimum-residual optimality at the continuous level. This is reflected in the top-right box of Figure 1, where the bolded outline indicates minimum-residual optimality.

where $g\left(\hat{\bm{{v}}}\right)\vcentcolon=\|\bm{{\Phi}}\hat{\bm{{v}}}-\bm{{f}}(\bm{{x}}_{0}+\bm{{\Phi}}{\hat{\bm{{x}}}},t)\|_{2}^{2}$. We now assess whether Eq. (3.4) holds, i.e., whether $\frac{d{{\hat{\bm{{x}}}}}}{dt}$ as defined by Eq. (3.2) is the minimizer of $g$.

The function $g$ can be expressed as $g\left(\hat{\bm{{v}}}\right)=\hat{\bm{{v}}}^{T}\bm{{\Phi}}^{T}\bm{{\Phi}}\hat{\bm{{v}}}-2\hat{\bm{{v}}}^{T}\bm{{\Phi}}^{T}\bm{{f}}(\bm{{x}}_{0}+\bm{{\Phi}}{\hat{\bm{{x}}}},t)+\bm{{f}}(\bm{{x}}_{0}+\bm{{\Phi}}{\hat{\bm{{x}}}},t)^{T}\bm{{f}}(\bm{{x}}_{0}+\bm{{\Phi}}{\hat{\bm{{x}}}},t)$. Due to the strict convexity of the function $g$, the global minimizer $\hat{\bm{{v}}}^{\star}$ is equal to the stationary point of $g$, i.e., $\hat{\bm{{v}}}^{\star}$ satisfies

where orthogonality of $\bm{{\Phi}}$ has been used. Comparing Eqs. (3.6) and (3.2) shows $\frac{d{{\hat{\bm{{x}}}}}}{dt}(\bm{{x}}_{0}+\bm{{\Phi}}{\hat{\bm{{x}}}},t)=\hat{\bm{{v}}}^{\star}$, which is the desired result. $\square$

Thus, Galerkin ROMs exhibit desirable properties (i.e., minimum-residual optimality) at the time-continuous level. We now derive a time-discrete representation for the Galerkin ROM, noting that these properties are lost at the time-discrete level.

2 Discrete representation

As before, a time-discretization method is needed to numerically solve Eq. (3.2). We now characterize the O$\Delta$E for the Galerkin ROM.

A linear $k$-step method applied to numerically solve Eq. (3.2) can be expressed as

Here, the O$\Delta$E is characterized by the following system of algebraic equations to be solved at each time step:

Here, the discrete residual corresponds to

and the generalized state is explicitly updated as

2.2 Runge–Kutta schemes

Applying an $s$-stage Runge–Kutta method to solve Eq. (3.2) leads to an O$\Delta$E characterized by the following system of algebraic equations to be solved at each time step:

Here, discrete the residual is defined as

and the generalized state is computed explicitly as

Note that the Galerkin-ROM solution satisfying Eqs. (3.7) or (3.9) does not in general associate with the solution to an optimization problem; therefore, the optimality property the method exhibits at the continuous level has been lost at the discrete level. We now show that Galerkin projection and time discretization are commutative; this implies that Galerkin ROMs can be analyzed, implemented, and interpreted equivalently at both the time-discrete and time-continuous levels. This corresponds to the rightmost part of Figure 1.

Performing a Galerkin projection on the governing ODE and subsequently applying time discretization yields the same model as first applying time discretization on the governing ODE and subsequently performing Galerkin projection.

Least-squares Petrov–Galerkin ROM

We note that other residual-minimizing approaches have been developed in the case of linear and nonlinear steady-state problems, and space–time solutions . In addition, a recently developed approach has suggested $L^{1}$ minimization of the residual arising at each time instance for hyperbolic problems.

We begin by developing the time-discrete representation for the LSPG ROM for both linear multistep schemes and Runge–Kutta schemes. The latter is a novel contribution, as previous work has derived discrete-optimal LSPG ROMs only for linear multistep schemes . Optimality of this approach corresponds to the bolded bottom-left box of Figure 1.

Note that necessary conditions for the solution to Eq. (4.2) corresponds to stationary of the objective function in Eq. (4.2), i.e., the solution satisfies

1.2 Runge–Kutta schemes

where $[\cdot]_{ij}$ denotes entry $(i,j)$ of the argument. This again leads to a least-squares Petrov–Galerkin interpretation for the discrete-optimal ROM.

We can then formulate the following sequence of optimization problems to compute discrete minimum-residual approximations:

Here, the associated Petrov–Galerkin projection is

Note that in the explicit case, the LSPG ROM generally requires solving a system of nonlinear equations at each Runge–Kutta stage. Because ${\partial\bm{{q}}^{n}_{i}}/\partial\bm{{w}}=\bm{{I}}$, the system of equations is linear if $\bm{{A}}_{i}$ are constant matrices, and only an explicit solution update is required if $\bm{{A}}_{i}=\bm{{I}}$ and $\bm{{\Phi}}^{T}\bm{{\Phi}}=\bm{{I}}$.

We now see the critical tradeoff between the Galerkin and LSPG ROMs: Galerkin ROMs exhibit continuous (minimum-residual) optimality, while LSPG ROMs exhibit discrete optimality. Without further analysis, it is unclear which of these attributes is preferable. Numerical experiments (Section 7) and supporting error analysis (Section 6) will highlight the benefits of discrete optimality over continuous optimality in practice.

2 Continuous representation

Because the LSPG ROM introduces approximations at the discrete level, it is unclear whether it can be interpreted at the continuous level. In fact, it has not previously been shown that a continuous representation of the LSPG ROM even exists. We now show that an ODE representation of the LSPG ROM does indeed exist for both linear multistep schemes and Runge–Kutta schemes under certain conditions; however, the ODE depends on the time step used to define the LSPG ROM. This associates with the top-left section of the relationship diagram in Figure 1.

The LSPG ROM for linear multistep integrators is equivalent to applying a Petrov–Galerkin projection to the ODE with test basis (in matrix form)

and subsequently applying time integration with a linear multistep scheme with time step $\Delta t$ if $\bm{{A}}$ is a constant matrix and (at least) one of the following conditions holds:

$\beta_{j}=0$, $j\geq 1$ (e.g., a single-step method),

the velocity $\bm{{f}}$ is linear in the state, or

Case 1 Applying a linear multistep time integrator with the stated assumption of $\beta_{j}=0$, $j\geq 1$ to numerically solve Eq. (4.16) results in the following discrete equations to be solved at each time instance:

Pre-multiplying by $\bm{{\Psi}}(\hat{\bm{{y}}}^{n},t^{n})^{T}\bm{{\Phi}}$ yields discrete equations $\hat{}{\bm{{r}}}^{n}\left(\hat{\bm{{y}}}^{n}\right)=0$ with residual

Comparing Eqs. (4.18) and (2.4) reveals $\hat{}{\bm{{r}}}^{n}\left(\hat{}\bm{{w}}\right)=\bm{{\Psi}}(\hat{}\bm{{w}},t^{n})^{T}{\bm{{r}}}^{n}\left(\bm{{x}}_{0}+\bm{{\Phi}}\hat{}\bm{{w}}\right)$ and so the solution $\hat{\bm{{y}}}^{n}$ satisfies

Under the stated assumptions, we have $\partial{\bm{{r}}}^{n}/\partial\bm{{w}}(\bm{{x}})=\alpha_{0}\bm{{I}}-\Delta t\beta_{0}\frac{\partial\bm{{f}}}{\partial\bm{{\xi}}}(\bm{{x}},t^{n})$ and so the LSPG test basis $\bm{{\Psi}}^{n}$ defined in Eq. (4.4) is equal to the test basis in Eq. (4.14) evaluated at time instance $n$, i.e., $\bm{{\Psi}}^{n}(\hat{}\bm{{w}})=\bm{{\Psi}}(\hat{}\bm{{w}},t^{n})$. Therefore, the solution $\hat{}\bm{{w}}^{n}$ to the LSPG O$\Delta$E (4.3) satisfies

This shows that $\hat{}\bm{{w}}^{n}=\hat{\bm{{y}}}^{n}$, i.e., the solutions to the LSPG O$\Delta$E and the O$\Delta$E obtained after applying Petrov–Galerkin projection with test basis $\bm{{\Psi}}(\bm{{x}},t)$ defined by Eq. (4.14) to the full-order model ODE and subsequently applying time integration are equivalent under the stated assumptions, which is the desired result. Case 2 In this case, the test basis is independent of the state, i.e.,

Applying a linear multistep time integrator to solve Eq. (4.16) and subsequently pre-multiplying by the constant matrix $\bm{{\Psi}}(t^{n})^{T}\bm{{\Phi}}$ yields the following discrete equations arising at each time step

Comparing Eqs. (4.23) and (2.4) reveals $\hat{}{\bm{{r}}}^{n}\left(\hat{}\bm{{w}}\right)=\bm{{\Psi}}(t^{n})^{T}{\bm{{r}}}^{n}\left(\bm{{x}}_{0}+\bm{{\Phi}}\hat{}\bm{{w}}\right)$ and so the solution $\hat{\bm{{y}}}^{n}$ satisfies

Under these assumptions, we have $\partial{\bm{{r}}}^{n}/\partial\bm{{w}}=\alpha_{0}\bm{{I}}-\Delta t\beta_{0}\partial\bm{{f}}/\partial\bm{{\xi}}(\cdot,t^{n})$ and so the LSPG test basis $\bm{{\Psi}}^{n}$ defined in Eq. (4.4) is equal to the test basis in Eq. (4.21) at time instance $n$, i.e., $\bm{{\Psi}}^{n}(\hat{}\bm{{w}})=\bm{{\Psi}}(t^{n})$. Therefore, the LSPG O$\Delta$E (4.3) can be expressed as

This shows that $\hat{}\bm{{w}}^{n}=\hat{\bm{{y}}}^{n}$, i.e., the solutions to the LSPG O$\Delta$E and the O$\Delta$E obtained after applying Petrov–Galerkin projection with test basis $\bm{{\Psi}}(t)$ defined by Eq. (4.14) to the full-order model ODE and subsequently applying time integration are equivalent under the stated assumptions. Case 3 The assumption $\beta_{0}=0$ results in a constant test basis

Applying a linear multistep time integrator to solve Eq. (4.16) and subsequently pre-multiplying by the constant matrix $\bm{{\Psi}}^{T}\bm{{\Phi}}$ yields

which is to be solved at each time step with a residual defined as

As in Case 2, this leads to $\hat{}{\bm{{r}}}^{n}\left(\hat{}\bm{{w}}\right)=\bm{{\Psi}}^{T}{\bm{{r}}}^{n}\left(\bm{{x}}_{0}+\bm{{\Phi}}\hat{}\bm{{w}}\right)$. Because $\frac{\partial{\bm{{r}}}^{n}}{\partial\bm{{w}}}(\bm{{x}})=\alpha_{0}\bm{{I}}$, we also again have $\bm{{\Psi}}^{n}(\hat{}\bm{{w}})=\bm{{\Psi}}$. This leads to the desired result, as the O$\Delta$Es for the LSPG ROM and the ROM obtained after applying Petrov–Galerkin projection with test basis $\bm{{\Psi}}$ to the full-order model ODE and subsequently applying time integration both satisfy $\bm{{\Psi}}^{T}{\bm{{r}}}^{n}(\bm{{x}}_{0}+\bm{{\Phi}}\hat{}\bm{{w}}^{n})=0$ under the stated assumptions. $\square$

We now provide conditions under which the LSPG ROM for Runge–Kutta schemes can be expressed as an ODE.

The LSPG ROM for Runge–Kutta integrators is equivalent to applying a Petrov–Galerkin projection to the ODE with test basis (in matrix form)

and subsequently applying time integration if $\bm{{A}}_{i}=\bm{{A}}$, $\forall i$ are constant matrices and the integrator is a singly diagonally implicit Runge–Kutta (SDIRK) scheme, i.e., $a_{ij}=0$, $\forall j>i$ and $a_{ii}=a$, $\forall i$.Note that summation on repeated indicies is not implied.

Applying an SDIRK time integrator to numerically solve Eq. (4.30) results in the following sequence of discrete equations to be solved at each time step:

Pre-multiplying by $\bm{{\Psi}}({\hat{\bm{{x}}}}^{n-1}+\Delta t\sum_{j=1}^{i}a_{ij}\bm{{\hat{\bm{{y}}}}}_{j},t^{n-1}+c_{i}\Delta t)^{T}\bm{{\Phi}}$ yields the following discrete equations

such that the solutions $\hat{\bm{{y}}}_{i}^{n}$ satisfy

Now, under the stated assumptions, we have

such that the LSPG test basis $\bm{{\Psi}}_{i}^{n}$ defined in Eq. (4.13) is related to the test basis in Eq. (4.29) as follows:

Therefore, the solutions $\hat{}\bm{{w}}^{n}_{i}$ to the LSPG O$\Delta$E (4.12) satisfy

We now show that the LSPG ROM has a time-continuous representation for all explicit and single-state Runge–Kutta schemes, which include the widely used forward Euler, backward Euler, and implicit midpoint schemes.

The LSPG ROM for Runge–Kutta integrators is equivalent to applying a Petrov–Galerkin projection to the ODE with test basis (in matrix form)

and subsequently applying time integration if $\bm{{A}}_{i}=\bm{{A}}$, $\forall i$ are constant matrices and an explicit Runge–Kutta scheme is employed.

Explicit Runge–Kutta schemes are characterized by $a_{ij}=0$, $j\geq i$ and so they satisfy the conditions Theorem 4.3 with $a=0$. $\square$

The LSPG ROM for Runge–Kutta integrators is equivalent to applying a Petrov–Galerkin projection to the ODE with test basis (in matrix form)

and subsequently applying time integration if $\bm{{A}}_{i}=\bm{{A}}$, $\forall i$ are constant matrices and a single-stage Runge–Kutta scheme is employed.

Single-stage Runge–Kutta schemes are characterized by $s=1$ and so they satisfy the conditions of Theorem 4.3 with $a=a_{11}$. $\square$

Equivalence conditions

This section performs theoretical analysis that highlights cases in which Galerkin and LSPG ROMs are equivalent, in which case the Galerkin ROM exhibits both continuous and discrete optimality. This suggests that time discretization and minimum-residual projection are commutative in these scenarios. Section 5.1 shows that equivalence holds for explicit time integrators, Section 5.2 demonstrates equivalence in the limit of $\Delta t\rightarrow 0$, and Section 5.3 shows equivalence in the case of symmetric-positive-definite residual Jacobians.

Galerkin projection is equivalent to LSPG projection with $\bm{{A}}=\frac{1}{\sqrt{\alpha_{0}}}\bm{{I}}$ for explicit linear multistep schemes.

In the case of explicit linear multistep schemes, $\beta_{0}=0$ and so Galerkin projection corresponds to Case 3 of Theorem 4.2 with $\bm{{A}}=\frac{1}{\sqrt{\alpha_{0}}}\bm{{I}}$, as $\bm{{\Psi}}=\bm{{\Phi}}$ in this case. $\square$

Galerkin projection is equivalent to LSPG projection with $\bm{{A}}=\bm{{I}}$ for explicit Runge–Kutta schemes.

In the case of explicit Runge–Kutta schemes, $a_{ij}=0$ $\forall j\geq i$ and so Galerkin projection corresponds to Theorem 4.3 with $\bm{{A}}=\bm{{I}}$ and $a=0$, as $\bm{{\Psi}}=\bm{{\Phi}}$ in this case. $\square$

2 Equivalence in the limit of Δ​t→0→Δ𝑡0\Delta t\rightarrow 0

Linear multistep schemes. Consider solving the LSPG O$\Delta$E (4.3) with $\bm{{A}}=\frac{1}{\sqrt{\alpha_{0}}}\bm{{I}}$. Then, the test basis defined in Eq. (4.4) is simply

From Eq. (2.4), we can write the residual Jacobian as

and so in the limit of $\Delta t\rightarrow 0$, the LSPG ROM solution satisfies

Now, from Eq. (2.6) the Jacobian can be expressed as

and so in the limit of $\Delta t\rightarrow 0$, the LSPG ROM solution satisfies

Because the Galerkin ROM solution also satisfies Eq. (5.2) (see Eq. (3.13) of Theorem 3.4), the two techniques are equivalent in this limit, which is the desired result. $\square$

3 Equivalence for symmetric-positive-definite residual Jacobians

In the case of linear multistep schemes, Galerkin projection is equivalent to LSPG projection with $\bm{{A}}\left(\bm{{z}}\right)=\bm{{U}}\left(\bm{{z}}\right)$, where $\bm{{U}}$ is the Cholesky factor Its derivative can be computed by solving the Lyapunov equation $\frac{\partial\bm{{U}}}{\partial w_{k}}^{T}\bm{{U}}+\bm{{U}}\frac{\partial\bm{{U}}}{\partial w_{k}}=-\left[\frac{\partial{\bm{{r}}}^{n}}{\partial\bm{{w}}}\right]^{-1}\frac{\partial^{2}{\bm{{r}}}^{n}}{\partial\bm{{w}}\partial w_{k}}\left[\frac{\partial{\bm{{r}}}^{n}}{\partial\bm{{w}}}\right]^{-1}$. of the residual-Jacobian inverse

if $\partial{\bm{{r}}}^{n}/\partial\bm{{w}}\left(\bm{{w}}^{n},t^{n}\right)=\alpha_{0}\bm{{I}}-\Delta t\beta_{0}\frac{\partial\bm{{f}}}{\partial\bm{{\xi}}}\left(\bm{{w}}^{n},t^{n}\right)$ is symmetric positive definite and if

Under the stated assumptions, the LSPG test basis defined in Eq. (4.4) is equal to the trial basis, i.e., $\bm{{\Psi}}^{n}(\hat{}\bm{{w}}^{n})=\bm{{\Phi}}$. By invoking Eq. (3.12), we can see that the O$\Delta$Es for the the LSPG ROM (4.3) and Galerkin ROM (3.7) both satisfy $\bm{{\Phi}}^{T}{\bm{{r}}}^{n}\left(\bm{{x}}_{0}+\bm{{\Phi}}\hat{}\bm{{w}}^{n}\right)=0$, which is the desired result. $\square$

In the case of diagonally implicit Runge–Kutta schemes, Galerkin projection is equivalent to LSPG projection with $\bm{{A}}_{i}\left(\bm{{z}}\right)=\underline{\bm{{U}}}_{i}\left(\bm{{z}}\right)$, where $\underline{\bm{{U}}}_{i}$ is the Cholesky factor of the residual-Jacobian inverse

if $\partial\bm{{q}}^{n}_{i}/\partial\bm{{w}}\left(\bm{{w}}^{n},t^{n}\right)=\bm{{I}}-\Delta ta_{ii}\frac{\partial\bm{{f}}}{\partial\bm{{\xi}}}\left(\bm{{x}}^{n-1}+\Delta ta_{ii}\bm{{w}}^{n}+\Delta t\sum_{j=1}^{i-1}a_{ij}\bm{{w}}^{n}_{j},t^{n-1}+c_{i}\Delta t\right)$ is symmetric positive definite and if

In the case of Runge–Kutta schemes, Galerkin projection exhibits discrete optimality if $\partial\bar{}{\bm{{r}}}^{n}/\partial\bar{}\bm{{w}}\left(\bar{}\bm{{w}}^{n},t^{n}\right)$ is symmetric positive definite and if

First, note that solution $(\hat{}\bm{{w}}^{n}_{1},\ldots,\hat{}\bm{{w}}^{n}_{s})$ to the Galerkin O$\Delta$E (3.13) equivalently satisfies

under the assumed conditions. This objective function can be written equivalently as

where $\bm{{A}}_{ij}=\bar{\bm{{U}}}_{ij}$. Comparing Eqs. (5.11) and (4.6) reveals that the Galerkin ROM satisfies a slightly more general notion of discrete optimality than the LSPG schemes considered in this work. $\square$

This analysis demonstrates that Galerkin projection exhibits discrete optimality when the residual Jacobian is symmetric positive definite. This is aligned with recent work that has shown Galerkin projection to be effective for Lagrangian dynamical systems —which are characterized by symmetric-positive-definite residual Jacobians—due to the fact that Galerkin projection preserves properties such as symplectic time evolution and energy conservation. For these reasons, using Galerkin projection is sensible for problems exhibiting these characteristics.

Error analysis

Ultimately, we are interested in assessing the state-space error between the (computed) time-discrete ROM solution and the (unknown) time-continuous FOM solution. This error comprises two contributions: the state-space error between (1) the time-continuous FOM and time-discrete FOM solutions (i.e., time-discretization error), and (2) the time-discrete ROM and time-discrete FOM solutions. This section focuses on the latter and performs time-discrete state-space error analyses for Galerkin and LSPG ROMs applied to different time integrators.

Section 6.1 derives error bounds for the Galerkin and LSPG ROMs for linear multistep schemes. Here, Theorem 6.5 provides a posteriori bounds that depend on the ROM solution, while Theorem 6.11 reports a priori bounds. Section 6.2 provides a posteriori error bounds for the backward Euler scheme, as well as additional analyses that highlight the important role of the time step in the LSPG ROM, which is discussed in Remark 6.18. Section 6.3 derives ROM error bounds for Runge–Kutta schemes. Here, Theorem 6.19 provides a posteriori bounds, Corollary 6.20 specializes these results to explicit Runge–Kutta and DIRK schemes, and Theorem 6.21 and Corollary 6.23 report a priori bounds for Runge–Kutta schemes.

and $\bm{{x}}_{\star}^{k}\vcentcolon=\bm{{x}}^{k}-\bm{{x}}_{0}$. We also define the FOM residuals at time instance $n$ associated with the trajectories associated with the FOM, Galerkin ROM, and LSPG ROM O$\Delta$Es as

We define the Galerkin and LSPG operators as

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

We proceed by deriving a posteriori error bounds for the Galerkin and LSPG ROMs for linear multistep schemes. We assume Lipschitz continuity of $\bm{{f}}$ in the first argument:

It is sufficient to show bound (6.26), as the arguments for (6.25) are similar. Let $n$ be fixed but arbitrary, then subtracting Eq. (6.3) from Eq. (6.1) yields

where $\delta{{\bm{{r}}}}_{P}^{n-1}\vcentcolon=\bar{}{\bm{{r}}}_{\star}^{n}\mathinner{\left[\bm{{x}}_{\star}^{n-k},\dots,\bm{{x}}_{\star}^{n-1}\right]}-\bm{{\Phi}}\bar{{\bm{{r}}}}_{P}^{n}\mathinner{\left[{\hat{\bm{{x}}}}_{P}^{n-k},\dots,{\hat{\bm{{x}}}}_{P}^{n-1}\right]}$. Adding and subtracting $\bm{{f}}\mathinner{\left(\bm{{x}}_{0}+\bm{{\Phi}}{\hat{\bm{{x}}}}_{P}^{n}\right)}$ and applying the triangle inequality leads to

Invoking $\mathinner{\left(\bf A_{1}\right)}$, and using $\Delta t<\mathinner{\!\left\lvert\alpha^{n}_{0}\right\rvert}/\mathinner{\!\left\lvert\beta^{n}_{0}\right\rvert}\kappa$, we deduce

Next, we will estimate $\mathinner{\!\left\lVert\delta{{\bm{{r}}}}_{P}^{n-1}\right\rVert}_{2}$. Using the definition of $\bar{}{\bm{{r}}}_{\star}^{n}$, $\bar{{\bm{{r}}}}_{P}^{n}$ from (6.4) we derive

We now establish a result that aids interpretability, as it provides a connection between the terms in the a posteriori error bounds and the optimality properties of LSPG and Galerkin ROMs. First, we note that

due to the optimality property associated with orthogonal projection; this term appears in local a posteriori error bound (6.8).

If $\beta_{j}^{n}=0$, $j\geq 1$ (e.g., backward differentiation formulas), then

If additionally the LSPG ROM employs $\bm{{A}}=\bm{{I}}$, then

Under the stated conditions, ${{\bm{{r}}}}_{G}^{n}={{\bm{{r}}}}_{P}^{n}$ and the optimality result Eq. (6.18) holds, yielding the desired result. $\square$

Corollary 6.3 shows that discrete-residual minimization (i.e., LSPG projection) rather than continuous-residual minimization (i.e., Galerkin projection) can produce a smaller value of a term that appears in the a posteriori error bounds; for example, this term appears on the right-hand side of Eqs. (6.8)–(6.9). This can be interpreted as arising from the fact that discrete-residual minimization computes the LSPG ROM solution that minimizes the entire normed quantity, while continuous-residual minimization performs orthogonal projection of the velocity given the Galerkin ROM solution.

Under the assumptions of Corollary 6.3 and Theorem 6.1, the local a posteriori error bound for the LSPG ROM (6.9) is smaller than that for the Galerkin ROM (6.8).

respectively. Under the conditions of Corollary 6.3, inequality (6.22) holds and the right-hand side of (6.24) will be smaller than the right-hand side of (6.23). $\square$

This is an interesting result, as it provides conditions under which the LSPG ROM produces a smaller local a posteriori error bound than the Galerkin ROM. It also provides some theoretical justification for the numerical experiments in Section 7, which use the three-point backward-differentiation formula, wherein the LSPG uniformly outperforms the Galerkin ROM. We now extend to results of Theorem 6.1 to obtain global a posteriori error bounds.

Under the assumptions of Theorem 6.1, we have

We now derive a simplified variant of these bounds.

Under the assumptions of Theorem 6.1, if additionally $\Delta t\leq|\alpha_{0}^{\star}|(1-\epsilon)/(\kappa|\beta_{0}^{\star}|)$ with $0<\epsilon<1$. Then,

We proceed by proving bound (6.28); the proof for bound (6.27) is similar. First, we define

as well as a ‘path’ from time instance $n$ backward to the initial time by defining $\mathcal{L}(0)=n$ with

with $\bar{n}\leq n$ and $\mathcal{L}(\bar{n})=0$. Then, from local a posteriori error bound (6.9), we have

the relation $(1+x)^{n}\leq\exp(nx)$, and the following result (with $x=1+\kappa\Delta t|\beta^{\star}|/|\alpha^{\star}|$ and $y=1-\kappa\Delta t|\beta_{0}^{\star}|/|\alpha_{0}^{\star}|$): if $x\geq y$, then $(x-y)/y\leq\epsilon^{-1}(x-y)$ if and only if $y\geq\epsilon>0$. $\square$

We note that due to the optimality property associated with orthogonal projection, we can write bound (6.27) equivalently as

We now prove conditions under which the a posteriori error bound is independent of the time step $\Delta t$ and total number of time instances $n$; the bound is fixed for a given time $t^{n}$.

Under the assumptions of Theorem 6.6, if additionally $k|\alpha^{\star}|=|\alpha_{0}^{\star}|$ (e.g., backward Euler, where $k=1$ and $|\alpha^{\star}|=|\alpha_{0}^{\star}|=1$) then the global a posteriori error bounds (6.27)–(6.28) are independent of the time step and simplify to

The result can be derived by substituting $k|\alpha^{\star}|=|\alpha_{0}^{\star}|$ into inequalities (6.27) and (6.28). $\square$

We now prove conditions under which a posteriori error bounds (6.27)–(6.28) can be written in ‘residual form’, i.e., in terms of the discrete residual arising at each time step. This will enable the respective optimality properties of the Galerkin and LSPG ROMs to be compared in Remark 6.9.

If additionally $\bm{{A}}=\bm{{I}}$, then

Comparing inequalities (6.41) and (6.46) highlights the differences in how optimality affects the Galerkin and LSPG ROM error bounds. Writing these expressions more compactly under the conditions of Corollary 6.8 yields

1.2 A priori error bounds

We now derive a priori error bounds by slightly modifying the steps in the above proofs. The most significant difference in the subsequent results is that the oblique projection associated with the LSPG ROM no longer associates with residual minimization, as the argument of the operators corresponds to the full-order-model solution.

Under the assumptions of Corollary 6.10, we have

The result can be derived by following the same steps as Theorem 6.5 based on the local bounds in Corollary 6.10. $\square$

The result can be derived by following the same steps as Theorem 6.6 based on the local bounds in Corollary 6.10. $\square$

We now demonstrate conditions under which the a priori error bound is independent of the time step $\Delta t$.

Under the assumptions of Corollary 6.12, if additionally $k|\alpha^{\star}|=|\alpha_{0}^{\star}|$ (e.g., backward Euler) for the Galerkin ROM, and $k|\bar{\alpha}^{\star}|=|\bar{\alpha}_{0}^{\star}|$ (e.g., backward Euler) for the LSPG ROM, then the a priori error bounds (6.53)–(6.54) are independent of the time step and simplify to

As in Corollary 6.7, the result can be shown by substitution of the appropriate assumptions (i.e., $k|\alpha^{\star}|=|\alpha_{0}^{\star}|$ for the Galerkin ROM, $k|\bar{\alpha}^{\star}|=|\bar{\alpha}_{0}^{\star}|$ for the LSPG ROM) into inequalities (6.53) and (6.54). $\square$

Because the argument of $\bm{{f}}$ in the a priori error bounds corresponds to the full-order-model solution, it is not possible to relate such quantities to ROM residuals as was done in the case of a posteriori error bounds in Lemma 6.2. As such, it is not possible to associate the oblique projection of the LSPG ROM with minimizing any component of the a priori error bounds, as was shown for a posteriori error bounds in Corollary 6.3. However, it is possible to associate Galerkin projection with minimizing terms in the a priori error bounds, i.e., the following optimality result holds under the conditions of Corollary 6.12:

This is analogous to inequality (6.47) for a posteriori error bounds, where

2 Backward Euler

We first note that the backward Euler scheme satisfies $\beta_{j}^{n}=0$, $j\geq 1$ in Corollary 6.4. Thus, that result provides conditions under which the LSPG ROM has a lower local a posteriori error bound than the Galerkin ROM. Next, we specialize the global a posteriori error bounds in Theorem 6.5 to the case of the backward Euler scheme.

Under the assumptions of Theorem 6.1, for the backward Euler scheme we obtain

where $h\vcentcolon=1-\kappa\Delta t$. Note that the time-step condition corresponds to $\Delta t<1/\kappa$ in this case.

Because it is a single-step method, the linear multistep coefficients do not vary between time instances; as a result, the constants appearing in Theorem 6.5 are also time-step independent and are $h=1-\kappa\Delta t$, $\varepsilon_{0}=\Delta t/h$, $\varepsilon_{1}=0$, $\gamma_{0}=(\kappa\Delta t+1)/h$, and $\gamma_{1}=1/h$. Substituting these values into error bound (6.26) and noting that

which simplifies to bound (6.59). Derivation of bound (6.58) is identical and is thus omitted. $\square$

We now specialize the results of the time-step-independent global a priori error bounds in Corollary 6.13 to the backward Euler scheme.

Under the assumptions of Corollary 6.13—which can be satisfied by backward Euler as $k|\alpha^{\star}|=|\alpha_{0}^{\star}|$ for this scheme—we obtain the following for the backward Euler scheme:

We now derive a result that highlights the (surprising) role the time step $\Delta t$ plays in the LSPG error bound. As will be shown in the numerical experiments in Section 7, this theoretical results can have an important effect on the performance of the LSPG ROM in practice.

If $\bar{\bm{{x}}}$ solves an auxiliary problem that computes the full-space solution increment centered on the LSPG ROM trajectory

Here, $\mu^{j}\vcentcolon=\mathinner{\!\left\lVert\bm{{\Phi}}\Delta{\hat{\bm{{x}}}}_{P}^{j}-\Delta\bar{\bm{{x}}}^{j}\right\rVert}_{2}$ denotes the difference in solution increments at time instance $j$, where $\Delta{\hat{\bm{{x}}}}_{P}^{j}\vcentcolon={\hat{\bm{{x}}}}_{P}^{j}-{\hat{\bm{{x}}}}_{P}^{j-1}$ and $\Delta\bar{\bm{{x}}}^{j}\vcentcolon=\bm{{\bar{x}}}^{j}-\bm{{\Phi}}{\hat{\bm{{x}}}}_{P}^{j-1}$. We denote the relative solution increment at time instance $j$ by $\bar{\mu}^{j}\vcentcolon=\mu^{j}/\|\Delta\bar{\bm{{x}}}^{j}\|_{2}$.

Eq. (6.59) in conjunction with (6.17) (with $\mathinner{\!\left\lvert\beta_{0}^{n}\right\rvert}=1$ and the appropriate form of the discrete residual for backward Euler) implies

Lipschitz continuity of $\bm{{f}}$ leads to the bound (6.64). To obtain Eq. (6.65), we multiply and divide by $\|\Delta\bar{\bm{{x}}}^{n-j}\|_{2}$ for each term in the summation and use $\Delta\bar{\bm{{x}}}^{n-j}=\Delta t\bm{{f}}\mathinner{\left(\bm{{\bar{x}}}^{n-j}\right)}$. $\square$

Corollary 6.17 is useful in that it expresses the LSPG ROM error in terms of the (time-local) single-step errors incurred by projection along the LSPG ROM trajectory. In addition, this result highlights the critical role of the time step $\Delta t$ in the performance of the LSPG ROM; the following remark provides this discussion.

The time step $\Delta t$ in the error bound (6.65) for the LSPG ROM solution plays an important role. In particular, decreasing the time step produces both beneficial effects (bound decrease) and deleterious effects (bound increase), which we denote by ‘+’ and ‘-’, respectively as follows:

The time-discretization error decreases (this does not appear in the time-discrete error analysis above).

The number of overall time instances $n$ increases, so there are more terms in the summation.

The terms $\Delta t(1+\kappa\Delta t)$ and $1/{(h)^{j+1}}$ decrease.

The term $\bar{\mu}^{n-j}$ may increase or decrease, depending on the spectral content of the basis $\bm{{\Phi}}$.

We now discuss this final ambiguous effect. The term $\bar{\mu}^{n}$ can be interpreted as the relative error in solution increment over $\left[(n-1)\Delta t,n\Delta t\right]$. Clearly, the ability of the LSPG ROM to make $\bar{\mu}^{n}$ small depends on the spectral content of the basis $\bm{{\Phi}}$: if the basis only captures modes that evolve over long time scales, then $\bar{\mu}^{n}$ will be large (i.e., close to one), as the basis does not contain the ‘fast evolving’ solution components that change over a single time step. This suggests that the time step should be ‘matched’ to the spectral content of the reduced basis $\bm{{\Phi}}$. In Section 7.5 of the experiments, we explore this issue numerically, and demonstrate that the error bound is minimized for an intermediate value of the time step $\Delta t$.

We note that the above arguments do not hold for the Galerkin ROM, which is simply an ODE that does not depend on the time step. Instead, decreasing the time step should increase accuracy, as it has the effect of reducing the time-discretization error.

3 Runge–Kutta schemes

We rewrite Eqs. (2.5), (3.9), and (4.7) as

and set $\bm{{x}}_{\star}^{0}=\bm{{x}}_{G}^{0}=\bm{{x}}_{P}^{0}=\bm{{0}}$.

We first derive a posteriori error bounds.

If $\mathinner{\left(\bf A_{1}\right)}$ holds and $\Delta t$ is such that

for every $\bm{{x}},\bm{{y}}\geq 0$, if $\bm{{D}}\bm{{x}}\leq\bm{{y}}$ then $\bm{{x}}\leq\bm{{D}}^{-1}\bm{{y}}$,

Galerkin ROM. First we will show bound (6.70). Subtracting Eq. (6.68) from Eq. (6.67) and applying the triangle inequality yields

where $\delta\bm{{w}}_{G,i}^{n}\vcentcolon=\bm{{w}}_{\star,i}^{n}-\bm{{\Phi}}\hat{}\bm{{w}}_{G,i}^{n}$. Adding and subtracting \bm{{f}}^{n}_{i}\Big{(}\bm{{x}}_{0}+\bm{{\Phi}}{{\hat{\bm{{x}}}}_{G}}^{n-1}+\Delta t\sum_{j=1}^{s}a_{ij}\bm{{\Phi}}\hat{}\bm{{w}}_{G,j}^{n}\Big{)} and invoking assumption $\mathinner{\left(\bf A_{1}\right)}$, we deduce

Selecting $\Delta t$ small enough such that $(a)$ and $(b)$ hold yields

where $[\cdot]_{ij}$ denotes entry $(i,j)$ of the argument. From explicit state updates (2.7) and (3.11), we obtain

Finally, an induction argument produces the desired result (6.70).

LSPG ROM. We now prove bound (6.71). Subtracting (6.69) from (6.67) and applying the triangle inequality yields

Again selecting $\Delta t$ small enough such that $(a)$ and $(b)$ hold yields

The explicit state updates from (2.7) and (3.11) and the bound (6.73) yield

An induction argument yields the bound (6.71). $\square$

Under the assumptions of Theorem 6.19 for explicit RK $(\theta=1)$ and DIRK $(\theta=0)$ schemes, we have

For explicit and diagonally implicit Runge–Kutta schemes (4.12), we have $\bm{{\Psi}}_{ij}^{n}=0$ when $i\neq j$. The proof is then an immediate consequence of (6.71). $\square$

Note that the LSPG error bound (6.74) resembles the Galerkin error bound (6.70) much more closely than the previous LSPG bound (6.71), as explicit Runge–Kutta and DIRK schemes remove the coupling of the test basis across stages.

3.2 A priori error bounds

We now state the a priori versions of the Galerkin Runge–Kutta schemes (6.70) and the LSPG Runge–Kutta schemes (6.71).

If $\mathinner{\left(\bf A_{1}\right)}$ holds and $\Delta t$ is such that

where we have used the convention that the empty product is equal to one.

Similarly to Corollary 6.16, we now derive a time-step-independent variant of the error bound (6.75). A similar result can be shown for the LSPG bound (6.76); however, we omit this result for simplicity and instead will provide (more readily interpretable) a priori LSPG ROM error bounds for explicit Runge–Kutta and DIRK schemes in Corollary 6.24.

Under the assumptions of Theorem 6.21, assume additionally $\Delta t\leq(1-\omega)/(\kappa a^{\star})$ with $a^{\star}\vcentcolon=\|[|a_{ij}|]_{ij}\|_{2}$ and $0<\omega<1$. Then,

As in Corollary 6.16, by using an upper bound for the right hand side in (6.75) we obtain

We now derive a bound for the term $\sum\limits_{k=1}^{s}|b_{k}|\sum\limits_{i=1}^{s}[\bm{{D}}^{-1}]_{ki}$ as

We can now compute an upper bound on the constant in inequality (6.78) as

Then, setting $t^{n}=n\Delta t$ we obtain

As in Theorem 6.1, we have used the relation $(1+x)^{n}\leq\exp(nx)$, and the result (with x=1+\kappa\Delta t\Big{(}b^{\star}s^{3/2}-a^{\star}\Big{)} and $y=1-\kappa\Delta ta^{\star}$) that states if $x\geq y$, then $(x-y)/y\leq\omega^{-1}(x-y)$ if and only if $y\geq\omega>0$. Finally, substituting in (6.82) in (6.3.2) and combining the resulting expression with (6.78) yields the desired result. $\square$

As with linear multistep schemes, the rightmost term in the Galerkin a priori bound will always be smaller than that for the LSPG bound, as the former associates with an orthogonal projection error of a fixed vector. In addition, the LSPG bound depends on the LSPG ROM solution; while this dependence could be removed, the bound in its current form facilitates comparison with the Galerkin bound. We again notice the complex structure of the estimator in (6.76) compared to (6.75). To better understand the behavior the LSPG estimator in (6.76) we consider two subcases: explicit Runge–Kutta and DIRK schemes.

Under the assumptions of Theorem 6.21 for explicit RK $(\theta=1)$ and DIRK $(\theta=0)$ schemes, we have

For explicit and diagonally implicit Runge–Kutta schemes (4.12), we have $\bm{{\Psi}}_{ij}^{n}=0$ when $i\neq j$. The proof is then an immediate consequence of (6.76). $\square$

Owing to the fact that the explicit Runge–Kutta and DIRK schemes removes the coupling of the basis, we again notice that the LSPG error bound (6.83) resembles the Galerkin ROM error bound (6.75).

The proof follows closely that of Corollary 6.22. We begin by deriving an upper bound for the right hand side in (6.75) as

Finally, we compute an upper bound on the constant in inequality (6.85) as

As this result is identical to upper bound (6.80) in the proof of Corollary 6.22 with $\bar{a}^{\star}$ replacing $a^{\star}$, we obtain the desired result by applying the remaining steps as Corollary 6.22. $\square$

Numerical experiments

This section compares the performance of Galerkin and LSPG ROMs on a computational-fluid-dynamics (CFD) application using a basis constructed by proper orthogonal decomposition. These experiments highlight the importance of the previous analyses, in particular the limiting equivalence of Galerkin and LSPG ROMs (Theorem 5.3), superior accuracy of the LSPG ROM compared with the Galerkin ROM (Corollary 6.4), and performance improvement of the LSPG ROM when an intermediate time step is selected (Corollary 6.17 and Remark 6.18).

Note that these experiments could be carried out on any dynamical system yielding a system of nonlinear ODEs (2.1); we have selected compressible turbulent fluid dynamics due to both its wide interest and challenging nature: limited progress has been made to date on developing robust, accurate ROMs for such problems. The numerical experiments highlight this fact, as standard Galerkin ROMs generate unstable responses in all cases.

The Galerkin and LSPG ROMs are implemented in AERO-F , a massively parallel compressible-flow solver. AERO-F solves the steady or unsteady compressible Navier–Stokes equations with various closure models available for turbulent flow, and employs a second-order node-centered finite-volume scheme. For model-reduction algorithms, all linear least-squares problems and singular value decompositions are computed in parallel using the ScaLAPACK library .

The full-order model corresponds to an unsteady Navier–Stokes simulation of a two-dimensional open cavity with a length-to-depth ratio of 4.5 using AERO-F’s DES turbulence model (based on the Spalart–Allmaras one-equation model ) and a wall-function boundary condition applied on solid surface boundaries. The fluid domain is discretized by a mesh with 192,816 nodes and 573,840 tetrahedra (Figure 2). The two-dimensional geometry is discretized in three dimensions by considering a slab of thin, but finite thickness, in the $z$-direction; the resulting grid is one element wide and is created by extruding a distance that is 1% of the cavity length. The viscosity is assumed to be constant, and the Reynolds number based on cavity length is $2.97\times 10^{6}$, while the free-stream Mach number is 0.6. Due to the turbulence model and three-dimensional domain, the number of conservation equations per node is $6$, and therefore the dimension of the CFD model is $N=1,156,896$. Roe’s scheme is employed to discretize the convective fluxes, and a linear variation of the solution is assumed within each control volume, which leads to a second-order space-accurate scheme on general unstructured, multi-dimensional meshes; however, we employ an extended stencil that gives fifth-order formal order of accuracy (with uniform mesh spacing) on inviscid, one-dimensional problems. The viscous flux is discretized using a centered Galerkin scheme.

Flow simulations are performed within a time interval $t\in\left[0,T\right]$ with $T=12.5$ time units. We employ the second-order accurate implicit three-point backward differentiation formula, which is a linear multistep scheme characterized by $k=2$, $\alpha_{0}=1$, $\alpha_{1}=-4/3$, $\alpha_{2}=1/3$, $\beta_{0}=2/3$, $\beta_{1}=\beta_{2}=0$, for time integration; future work will perform numerical experiments with Runge–Kutta schemes. The O$\Delta$E (2.3) arising at each time step is solved by a Newton–Krylov method, where GMRES is employed as the iterative linear solver with a restrictive additive Schwarz preconditioner (with no fill in) and the previous 50 Krylov vectors are employed for orthogonalization. Convergence is declared when the residual norm is reduced to a factor of $10^{-3}$ of its starting value. All flow computations are performed in a non-dimensional setting.

The initial condition $\bm{{x}}_{0}$ is provided by first computing a steady-state solution, and using that solution as an initial guess for an unsteady ‘transient’ simulation (which captures the initial transient before the flow reaches a quasi-periodic state) of 7.5 time units. The state at the end of the unsteady transient simulation is then used as the initial condition for the subsequent simulations. The steady-state calculation is characterized by the same parameters as above, except that it employs local time stepping with a maximum CFL number of 100, it uses the first-order implicit backward Euler time integration scheme, it assumes a linear variation of the solution within each control volume, it employs a Spalart–Allmaras turbulence model, and it employs only one Newton iteration per (pseudo) time step.

All computations are performed in double-precision arithmetic on a parallel Linux clusterThe cluster contains 8-core compute nodes that each contain a 2.93 GHz dual socket/quad core Nehalem X5570 processor with 12 GB of memory. The interconnect is a 3D torus InfiniBand. using 48 cores across 6 nodes.

2 Time-step verification

Because this paper considers the time step to be an important parameter in model reduction, we first perform a time-step verification study to ensure we employ an appropriate ‘nominal’ time step. Figure 3 reports these results using a time-step refinement factor of two. A time step of $\Delta t_{\star}=0.0015$ time units yields observed convergence rates in both the instantaneous drag force on the lower wall and instantaneous pressure at $t=T$ that are close to the asymptotic rate of convergence (2.0) of three-point BDF2 scheme. Further, this value also leads to sub-2% errors in both quantities, which we deem to be sufficient for this set of experiments.

Figure 4 shows several instantaneous snapshots of the vorticity field and corresponding pressure field generated by the high-fidelity CFD model. The flow within the cavity is quasi-periodic; during one cycle, vorticity is shed from the leading edge of the cavity, convects downstream, and impinges on the aft edge of the cavity. Upon impingement, an acoustic disturbance is generated which propagates upstream and scatters on the leading edge of the cavity, generating a new vortical disturbance to initiate the next oscillation cycle. The pressure fields in the bottom row of Figure 4 reveal regions of low pressure (blue contours) associated with vortices, as well as acoustic disturbances both within the cavity and radiating outside the cavity. This complex flow is governed by the interactions of several nonlinear processes, including roll-up of the shear layer vortices, impingement of the vortices on the aft wall resulting in sound generation, propagation of nonlinear acoustic waves, and interaction of these waves with the shear layer vorticity.

3 Reduced-order models

To construct both the Galerkin and LSPG ROMs, we employ the proper orthogonal decomposition (POD) technique; we employ a constant weighting matrix $\bm{{A}}=\bm{{I}}$ for the LSPG ROM. To construct the POD basis, we set $\bm{{\Phi}}\leftarrow\bm{{\Phi}}\left(\mathcal{X},\nu\right)$, where $\bm{{\Phi}}$ is computed via Algorithm 1 of the appendix with snapshots consisting of the initial-condition-centered full-order model states $\mathcal{X}=\{\bm{{x}}_{\star}(k\Delta t_{\star})-\bm{{x}}_{0}\}_{k=1}^{8334}$, where $\bm{{x}}_{\star}$ denotes the FOM response computed for a time step of $\Delta t_{\star}=0.0015$. Three values of the energy criterion $\nu\in$ are used during the experiments: $\nu=1-10^{-4}$ (${p}=204$), $\nu=1-10^{-5}$ (${p}=368$), and $\nu=1-10^{-6}$ (${p}=564$). Figure 5 shows a selection of the energy component of the computed POD modes. Note that as the mode number increases, the modes capture finer spatial-scale behavior, which we expect to be associated with finer time-scale behavior; this will be verified in Section 7.5.1.

We first repeat the time-step verification study, but we do so for the reduced-order models (again using the BDF2 scheme) in the time interval $0\leq t\leq 0.55$, as all Galerkin ROMs remain stable in this time interval. Figure 6 reports these results. First, we note that the Galerkin ROM converges an approximated rate of 2.0, which is what we expect given that the Galerkin ROM simply associates with a time-step-independent ODE (3.2). However, the LSPG ROM does not exhibit this behavior; in fact the error convergence is not even monotonic. This is likely due to the fact that the method does not associate with a time-step-independent ODE.

We next perform simulations for both reduced-order models for all tested basis dimensions and time steps; Figure 7 reports the time-dependent responses. When a response stops before the end of the time interval, this indicates that a negative pressure was encountered, which causes AERO-F to exit the simulation. We interpret this phenomenon as a non-physical instability.

First, note that the Galerkin ROMs become unstable (i.e., generate a negative pressure) for all time steps and all basis dimensions. This is consistent with previously reported results that indicate Galerkin projection almost always leads to inaccurate responses for compressible fluid-dynamics problems. In contrast, the LSPG ROM results in many stable, accurate responses for all basis dimensions. Further, LSPG responses exhibit a clear dependence on the time step $\Delta t$. Subsequent sections provide a deeper analysis of this dependence.

4 Limiting case: comparison

We next compare the responses of the Galerkin and LSPG ROMs for small time windows (when the Galerkin responses remain stable) and small time steps. Figure 8 reports $\varepsilon(p_{\text{discrete~{}opt.}},p_{\text{Gal}_{\star}})$—which is the difference between the pressure responses generated by the LSPG ROM with different time steps and the Galerkin ROM with a fixed time step $\Delta t=1.875\times 10^{-4}$ (the smallest tested time step)—for a time window $0\leq t\leq 1.1$. These responses support an important conclusion (see Theorem 5.3): the Galerkin and LSPG ROMs are equal in the limit of $\Delta t\rightarrow 0$ for $\bm{{A}}=1/\sqrt{\alpha_{0}}\bm{{I}}$, which is what we employ for the LSPG ROM (note that $\alpha_{0}=1$ for this time integrator).Note that in the ${p}=564$ case, it is not clear if the difference is converging to zero. This is likely due to the fact that the time steps are not sufficiently small to detect convergence to zero in this case. In fact, as the basis dimension ${p}$ increases, the basis captures finer temporal behavior (as will be shown in Figure 10) and so the time scale of the ROM response will be smaller; in turn, smaller time steps $\Delta t$ will be required to detect convergent behavior. This has significant consequences for the LSPG ROM, as decreasing the time step leads to the same unstable response as Galerkin; larger time steps are needed to ensure the LSPG ROM is stable for the entire time interval.

Figure 9 reports $\varepsilon(p_{\text{discrete~{}opt.}},p_{\text{FOM}_{\star}})$ and $\varepsilon(p_{\text{Gal.}},p_{\text{FOM}_{\star}})$—which are the differences between the two ROM-generated pressure responses and the full-order model pressure response for $\Delta t=1.875\times 10^{-4}$—as a function of the time step for all three basis dimensions and three time intervals. These results highlight a critical observation: the LSPG ROM is more accurate for an intermediate time step. This not only supports the result of Corollary 6.17, but provides an interesting insight: taking a larger time step not only leads to better speedups (i.e., the end of the time interval is reached in fewer time steps), but it also decreases the error, sometimes significantly. This is further explored in the next section.

5 Time-step selection

Recall from Corollary 6.17 and Remark 6.18 that decreasing the time step $\Delta t$ has a non-obvious effect on the error bound for the LSPG ROM. We now assess these effects for the current problem.

In our interpretation of the error bound (6.65) for the LSPG ROM applied to the backward Euler scheme, we noted that the time step should be ‘matched’ to the spectral content of the trial basis $\bm{{\Phi}}$. This is of practical importance, as selecting an appropriate time step for the ROM should take into account the relevant temporal dynamics associated with the basis. For example, a time step may be too small if the basis has filtered out modes with a time scale matching that of the time step. If we assume that the basis $\bm{{\Phi}}$ is computed via POD, then we would expect the vectors to be naturally ordered such that lower mode numbers are associated with lower temporal frequencies. Then, including additional modes has the effect of encoding information at higher frequencies. It follows that the time step should be decreased as additional modes are retained in construction of the ROM.

Thus, at least for the present application problem, we expect the optimal time step for the LSPG ROM to decrease as modes are added to the POD basis (this will be verified by Figure 12). Note that systematic calibration could be performed to attempt to automate selection of the ROM time step as a function of basis dimension. While this would be of clear practical interest, we do not pursue it here, as optimal-timestep computation would be complicated in practice by nonlinear interactions arising from the dynamical system, as well as effects from the spatial-discretization error and POD truncation error.

5.2 Error bound behavior

Having verified that higher POD mode numbers correspond to smaller wavelengths, we now numerically assess quantities related to the error bound (6.65). First, Figure 11(a) reports the dependence of the maximum relative projection error $\max_{k}\bar{\mu}_{\star}^{k}(\bm{{\Phi}},\Delta t)$ on the time step $\Delta t$ and the basis dimension, where

Note that $\bar{\mu}_{\star}^{k}$ is closely related to $\bar{\mu}^{k}$ from error bound (6.65), as they are equal if $\bm{{x}}_{0}+\bm{{\Phi}}{\hat{\bm{{x}}}}_{P}(t)=\bm{{x}}(t)$ and the LSPG ROM computes ${\hat{\bm{{x}}}}_{P}^{k}$ such that $\bar{\mu}^{k}$ is minimized.

These results confirm that adding basis vectors—which we know has the effect of encoding higher frequency content—significantly reduces the projection error for small time steps $\Delta t$, but has less of an effect on larger time steps, as retaining the first POD vectors already enables dynamics at that scale to be captured.

Next, Figure 11(b) plots the error bound (6.65) for a value of $\kappa=1$ and with $\bar{\mu}^{k}=\bar{\mu}_{\star}^{k}$. This highlights an important result: selecting an intermediate time step $\Delta t$ leads to the lowest error bound, regardless of the basis dimension. Even though this result corresponds to the backward Euler integrator, we expect a similar trend to hold for the present experiment, which uses the BDF2 scheme. The next section assesses the performance of the LSPG ROM, including its dependence on the time step.

6 LSPG ROM performance

We now compare the accuracy and walltime performance of the LSPG ROM as the dimension of the basis, time step, and time interval change. The most salient result from Figure 12 is that choosing an intermediate time step leads to both better accuracy and faster simulation times. This shows that our theoretical analysis of the error bound performed in Section 7.5.2 leads to an actual observed performance improvement. For example, consider the ${p}=564$ case over the time interval $0\leq t\leq 2.5$. In this case, a time step of $\Delta t=1.875\times 10^{-4}$ leads to a relative error of $0.0140$ and a simulation time of $289$ hours; increasing this value to $\Delta t=1.5\times 10^{-3}$ reduces the relative error to $9.46\times 10^{-4}$ and the simulation time to $35.8$ hours, which constitutes roughly an order of magnitude improvement in both quantities. Again, this supports the theoretical results of Corollary 6.17 and highlights the critical importance of the time step for LSPG reduced-order models.

In addition, Figure 12 shows that as the basis dimension increases, the optimal time step decreases; this was anticipated from the spectral analysis performed in Section 7.5.1. In addition, adding POD basis vectors does not improve accuracy for large time steps. We interpret this effect as follows: for larger time steps, the first few POD modes accurately capture ‘coarse’ phenomena on the scale of the time step. Therefore, accuracy improvement is not achieved by adding modes that encode dynamics that evolve on a time scale finer than the time step itself.

Further, Figure 12(g) highlights that as the basis dimension increases, the error generally decreases, which is an artifact of the monotonic decrease in the FOM O$\Delta$E residual achieved by the LSPG ROM (Remark 4.1). Finally, the figure shows that as the time interval grows, the optimal time step generally increases.

7 GNAT: ROM with complexity reduction

In this section, we perform a similar study, but equip the LSPG ROM with complexity reduction in order to achieve computational savings. In particular, we employ the GNAT method , which solves Eq. (4.1) with $\bm{{A}}=\left(\bm{{P}}\bm{{\Phi}}_{r}\right)^{+}\bm{{P}}$, where $\bm{{\Phi}}_{r}$ is a basis for the residual and $\bm{{P}}$ consisting of selected rows of the identity matrix.

The problem is identical to that described in Section 7.1 except that we take $T=5.5$ time units and employ a second-order space-accurate dissipation scheme wherein a linear variation of the solution is assumed within each control volume.This is done to ensure the sample mesh requires two layers of neighboring nodes for each sample node. For this simulation, the full-order model consumes 5.0 hours on 48 cores across six compute nodes.

The GNAT implementation in AERO-F is characterized by the sample-mesh concept . Figure 13 depicts the sample mesh for this problem, which was constructed using $n_{c}=2228$ working columns [21, Algorithm 3], and includes two layers of nodes around the sample nodes (to enable the residual to be computed at the sample nodes). It is characterized by 7,974 total nodes (4.1% of the original mesh) and 17,070 total volumes (3.0% of the original mesh). Due to the small footprint of the sample mesh, the GNAT simulations are run using only 2 cores on a single compute node.

Figure 14 reports the results obtained with the GNAT ROM using different time steps. Critically, note that the GNAT ROM also exhibits a ‘dip’ in the optimal time step, with a time step of $6.0\times 10^{-3}$ yielding the lowest error. In fact, increasing the time step from $1.5\times 10^{-3}$ to $6.0\times 10^{-3}$ decreases the error from 3.32% to 2.25% and also significantly increases the computational savings relative to the full-order model (as measured in core–hours) from 14.9 to 55.7. This highlights that the analysis is also relevant to ROMs equipped with complexity reduction.

8 Summary of experimental results

We now briefly summarize the main experimental results:

Galerkin ROMs are unstable for long time intervals (Figure 7).

LSPG ROMs are only unstable for small time steps (Figure 7).

Galerkin and LSPG ROMs are equivalent as $\Delta t\rightarrow 0$ (Figure 8).

LSPG ROMs are more accurate than Galerkin ROMs over small time windows where Galerkin is stable (Figure 9).

LSPG ROMs are most accurate for an intermediate time step (Figure 9).

Adding POD modes has the effect of including higher-frequency response components (Figure 10).

The theoretical error bound for the LSPG ROM exhibits the same time step ‘dip’ as the experimentally observed error (Figure 11).

The optimal time step for the LSPG ROM decreases as modes are added to the POD basis (Figure 12).

Adding modes to the POD basis has little effect on LSPG ROM accuracy for large time steps (Figure 12).

The optimal time step for the LSPG ROM tends to increase as the time interval increases (Figure 12(g)).

The GNAT ROM, which is discrete optimal and is equipped with complexity reduction, also produces minimal error for an intermediate time step (Figure 14).

Conclusions

This work has performed a comparative theoretical and experimental analysis of Galerkin and LSPG reduced-order models for linear multistep schemes and Runge–Kutta schemes. We have demonstrated a number of new findings that have important practical implications, including conditions under which the LSPG ROM has a time-continuous representation, conditions under which the two techniques are equivalent, and time-discrete error bounds for the two approaches.

Perhaps most surprisingly, we demonstrated that decreasing the time step does not necessarily decrease the error for the LSPG ROM. This phenomenon arose in both the theoretical analysis and in numerical experiments. In particular, our results suggest that the time step should be ‘matched’ to the spectral content of the reduced basis. In the experiments, we showed that increasing the time step to an intermediate value decreased both the error and the simulation time by an order of magnitude in certain cases. Alternatively, decreasing the time step cause the LSPG ROM to become unstable for longer time intervals. This highlights the critical importance of time-step selection for LSPG ROMs.

Acknowledgments

We thank Prof. Stephen Pope for insightful conversations related to comparing Galerkin and least-squares Petrov–Galerkin reduced-order models; these conversations inspired this work. We also thank the anonymous reviewers for their extremely helpful and insightful comments and suggestions. We also thank Prof. Charbel Farhat for permitting us the use of AERO-F, as well as Julien Cortial, Charbel Bou-Mosleh, and David Amsallem for their previous contributions in implementing nonlinear reduced-order models in AERO-F. K. Carlberg acknowleges an appointment to the Sandia National Laboratories Truman Fellowship in National Security Science and Engineering. The Truman Fellowship is sponsored by Sandia National Laboratories. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000. H. Antil acknowledges the support by the NSF-DMS-1521590. The content of this publication does not necessarily reflect the position or policy of any of these institutions, and no official endorsement should be inferred.

Appendix

Algorithm 1 reports the algorithm for computing a POD basis using normalized snapshots.

References