The Adjoint Petrov-Galerkin Method for Non-Linear Model Reduction

Introduction

High-fidelity numerical simulations play a critical role in modern-day engineering and scientific investigations. The computational cost of high-fidelity or full-order models (FOMs) is, however, often prohibitively expensive. This limitation has led to the emergence of reduced-order modeling techniques. Reduced-order models (ROMs) are formulated to approximate solutions to a FOM on a low-dimensional manifold. Common reduced-order modeling techniques include balanced truncation , Krylov subspace techniques , reduced-basis methods , and the proper orthogonal decomposition approach . Reduced-order models based on such techniques have been implemented in a wide variety of disciplines and have been effective in reducing the computational cost associated with high-fidelity numerical simulations .

Projection-based reduced-order models constructed from proper orthogonal decomposition (POD) have proved to be an effective tool for model order reduction of complex systems. In the POD-ROM approach, snapshots from a high-fidelity simulation (or experiment) are used to construct an orthonormal basis spanning the solution space. A small, truncated set of these basis vectors forms the trial basis. The POD-ROM then seeks a solution within the range of the trial basis via projection. Galerkin projection, in which the FOM equations are projected onto the same trial subspace, is the simplest type of projection. The Galerkin ROM (G ROM) has been used successfully in a variety of problems. When applied to general non-self-adjoint and non-linear problems, however, theoretical analysis and numerical experiments have shown that Galerkin ROM lacks a priori guarantees of stability, accuracy, and convergence . This last issue is particularly challenging as it demonstrates that enriching a ROM basis does not necessarily improve the solution . The development of stable and accurate reduced-order modeling techniques for complex non-linear systems is the motivation for the current work.

A significant body of research aimed at producing accurate and stable ROMs for complex non-linear problems exists in the literature. These efforts include, but are not limited to, “energy-based” inner products , symmetry transformations , basis adaptation , $L^{1}$-norm minimization , projection subspace rotations , and least-squares residual minimization approaches . The Least-Squares Petrov–Galerkin (LSPG) method comprises a particularly popular least-squares residual minimization approach and has been proven to be an effective tool for non-linear model reduction. Defined at the fully-discrete level (i.e., after spatial and temporal discretization), LSPG relies on least-squares minimization of the FOM residual at each time-step. While the method lacks a priori stability guarantees for general non-linear systems, it has been shown to be effective for complex problems of interest . Additionally, as it is formulated as a minimization problem, physical constraints such as conservation can be naturally incorporated into the ROM formulation . At the fully-discrete level, LSPG is sensitive to both the time integration scheme as well as the time-step. For example, in Ref. it was shown that LSPG produces optimal results at an intermediate time-step. Another example of this sensitivity is that, when applied to explicit time integration schemes, the LSPG approach reverts to a Galerkin approach. This limits the scope of LSPG to implicit time integration schemes, which can in turn increase the cost of the ROM It is possible to use LSPG with an explicit time integrator by formulating the ROM for an implicit time integration scheme, and then time integrating the resulting system with an explicit integrator.. This is particularly relevant in the case where the optimal time-step of LSPG is small, thus requiring many time-steps of an implicit solver. Despite these challenges, the LSPG approach is arguably the most robust technique that is used for ROMs of non-linear dynamical systems.

Research has examined the application of both phenomenological and residual-based subgrid-scale models to POD-ROMs. In Refs. , Iliescu and co-workers examine the construction of eddy-viscosity-based ROM closures via the VMS method. These eddy-viscosity methods are directly analogous to the eddy-viscocity philosophy used in turbulence modeling. While they do not guarantee stability a priori, these ROMs have been shown to enhance accuracy on a variety of problems in fluid dynamics. However, as eddy-viscosity methods are based on phenomenological assumptions specific to three-dimensional turbulent flows, their scope may be limited to specific types of problems. Residual-based methods, which can also be derived from VMS, constitute a more general modeling strategy. The subgrid-scale model emerging from a residual-based method typically appears as a term that is proportional to the residual of the full-order model; if the governing equations are exactly satisfied by the ROM, then the model is inactive. While residual-based methods in ROMs are not as well-developed as they are in finite element methods, they have been explored in several contexts. In Ref. , ROMs of the Navier-Stokes equations are stabilized using residual-based methods. This stabilization is performed by solving a ROM stabilized with a method such as streamline upwind Petrov–Galerkin (SUPG) and augmenting the POD basis with additional modes computed from the residual of the Navier-Stokes equations. In Ref. , residual-based stabilization is developed for velocity-pressure ROMs of the incompressible Navier-Stokes equations. Both eddy-viscosity and residual-based methods have been shown to improve ROM stability and performance. The majority of existing work on residual-based stabilization (and eddy-viscosity methods) is focused on ROMs formulated from continuous projection (i.e., projecting a continuous PDE using a continuous basis). In this instance, the ROM residual is defined at the continuous level and is directly linked to the governing partial differential equation. In many applications (arguably the majority ), however, the ROM is constructed through discrete projection (i.e., projecting the spatially discretized PDE using a discrete basis). In this instance, the ROM residual is defined at the semi-discrete level and is tied to the spatially discretized governing equations. Residual-based methods for ROMs developed through discrete projections have, to the best of the authors’ knowledge, not been investigated.

Another approach that displays similarities to the variational multiscale method is the Mori-Zwanzig (MZ) formalism. Originally developed by Mori and Zwanzig and reformulated by Chorin and co-workers , the MZ formalism is a type of model order reduction framework. The framework consists of decomposing the state variables in a dynamical system into a resolved (coarse-scale) set and an unresolved (fine-scale) set. An exact reduced-order model for the resolved scales is then derived in which the impact of the unresolved scales on the resolved scales appears as a memory term. This memory term depends on the temporal history of the resolved variables. In practice, the evaluation of this memory term is not tractable. It does, however, serve as a starting point to develop closure models. As MZ is formulated systematically in a dynamical system setting, it promises to be an effective technique for developing stable and accurate ROMs of non-linear dynamical systems. A range of research examining the MZ formalism as a multiscale modeling tool exists in the community. Most notably, Stinis and co-workers have developed several models for approximating the memory, including finite memory and renormalized models, and examined their application to the semi-discrete systems emerging from Fourier-Galerkin and Polynomial Chaos Expansions of Burgers’ equation and the Euler equations. Application of MZ-based techniques to the classic POD-ROM approach has not been undertaken.

This manuscript leverages work that the authors have performed on the use of the MZ formalism to develop closure models of partial differential equations . In addition to focusing on the development and analysis of MZ models, the authors have examined the formulation of the MZ formalism within the context of the VMS method . By expressing MZ models within a VMS framework, similarities were discovered between MZ and VMS models. In particular, it was discovered that several existing MZ models are residual-based methods.

The development of a novel projection-based reduced-order modeling technique, termed the Adjoint Petrov–Galerkin (APG) method. The method leads to a ROM equation that is driven by the residual of the discretized governing equations. The approach is equivalent to a Petrov–Galerkin ROM and displays similarities to the LSPG approach. The method can be evolved in time with explicit integrators (in contrast to LSPG). This potentially lowers the cost of the ROM.

Theoretical error analysis examining conditions under which the a priori error bounds in APG may be smaller than in the Galerkin method.

Computational cost analysis (in FLOPS) of the proposed APG method as compared to the Galerkin and LSPG methods. This analysis shows that the APG ROM is twice as expensive as the G ROM for a given time step, for both explicit and implicit time integrators. In the implicit case, the ability of the APG ROM to make use of Jacobian-Free Newton-Krylov methods suggests that it may be more efficient than the LSPG ROM.

Numerical evidence on ROMs of compressible flow problems demonstrating that the proposed method is more accurate and stable than the G ROM on problems of interest. Improvements over the LSPG ROM are observed in most cases. An analysis of the computational cost shows that the APG method can lead to lower errors than the LSPG and G ROMs for the same computational cost.

Theoretical results and numerical evidence that provides a relationship between the time-scale in the APG ROM and the spectral radius of the right-hand side Jacobian. Numerical evidence suggests that this relationship also applies to the selection of the optimal time-step in LSPG.

The structure of this paper is as follows: Section 2 outlines the full-order model of interest and its formulation in generalized coordinates. Section 3 outlines the reduced-order modeling approach applied at the semi-discrete level. Galerkin, Petrov–Galerkin, and VMS ROMs will be discussed. Section 3.3 details the Mori-Zwanzig formalism and the construction of the Adjoint Petrov–Galerkin ROM. Section 4 provides theoretical error analysis. Section 5 discusses the implementation and computational cost of the Adjoint Petrov–Galerkin method. Numerical results and comparisons with Galerkin and LSPG ROMs are presented in Section 6. Conclusions are provided in Section 7.

Mathematical notation in this manuscript is as follows: matrices are written as bold uppercase letters (e.g. ${\mathbf{V}}$), vectors as lowercase bold letters (e.g. ${\mathbf{u}}$), and scalars as italicized lowercase letters (e.g. $a_{i}$). Calligraphic script may denote vector spaces or special operators (e.g. $\mathcal{V}$, $\mathcal{L}$). Bold letters followed by parentheses indicate a matrix or vector function (e.g. $\mathbf{R}(\cdot)$, ${\mathbf{u}}(\cdot)$), and those followed by brackets indicate a linearization about the bracketed argument (e.g. $\mathbf{J}[\cdot]$).

Full-Order Model and Generalized Coordinates

Consider a full-order model that is described by the dynamical system,

In many practical applications, the computational cost associated with solving Eq. 1 is prohibitively expensive due to the high dimension of the state. The goal of a ROM is to transform the $N$-dimensional dynamical system presented in Eq. 1 into a $K$ dimensional dynamical system, with $K\ll N$. To achieve this goal, we pursue the following agenda:

Develop a weak form of the FOM in generalized coordinates.

Decompose the generalized coordinates into a $K$-dimensional resolved coarse-scale set and an $N-K$ dimensional unresolved fine-scale set.

Develop a $K$-dimensional ROM for the coarse-scales by making approximations to the fine-scale coordinates.

The remainder of this section will address task 1 in the above agenda.

To develop the weak form of Eq. 1, we start by defining a trial basis matrix comprising $N$ orthonormal basis vectors,

Equation 1 can be expressed in terms of the generalized coordinates by inserting Eq. 2 into Eq. 1,

The weak form of Eq. 3 is obtained by taking the $L^{2}$ inner product with ${\mathbf{W}}$,The authors recognize that many types of inner products are possible in formulating a ROM. To avoid unnecessary abstraction, we focus here on the simplest case.

Manipulation of Eq. 4 yields the following dynamical system,

Reduced-Order Models

This subsection addresses task 2 in the aforementioned agenda. Reduced-order models seek a low-dimensional representation of the original high-fidelity model. To achieve this, we examine a multiscale formulation of Eq. 5. Consider sum decompositions of the trial and test space,

The fine-scale space is a subspace of $\mathcal{V}$, i.e., ${\mathcal{V}^{\prime}}\subset{\mathcal{V}}$.

For notational purposes, we make the following definitions for the trial and test spaces:

where $[\cdot\hskip 3.61371pt;\hskip 3.61371pt\cdot]$ denotes the concatenation of two matrices and,

The coarse and fine-scale states are defined as,

Equation 7 is referred to as the coarse-scale equation, while Eq. 8 is referred to as the fine-scale equation. It is important to emphasize that the system formed by Eqs. 7 and 8 is still an exact representation of the original FOM.

The objective of ROMs is to solve the coarse-scale equation. The challenge encountered in this objective is that the evolution of the coarse-scales depends on the fine-scales. This is a type of “closure problem” and must be addressed to develop a closed ROM.

2 Reduced-Order Models

As noted above, the objective of a ROM is to solve the (unclosed) coarse-scale equation. We now develop ROMs of Eq. 1 by leveraging the multiscale decomposition presented above. This section addresses task 3 in the mathematical agenda.

The most straightforward technique to develop a ROM is to make the approximation,

This allows for the coarse-scale equation to be expressed as,

Equation 9 forms a $K$-dimensional reduced-order system (with $K\ll N$) and provides the starting point for formulating several standard ROM techniques. The Galerkin and Least-Squares Petrov–Galerkin ROMs are outlined in the subsequent subsections.

When applied to unsteady non-linear problems, the Galerkin ROM is often inaccurate and, at times, unstable. Examples of this are seen in Ref. . These issues motivate the development of more sophisticated reduced-order modeling techniques.

2.2 Petrov–Galerkin and Least-Squares Petrov–Galerkin Reduced-Order Models

In writing Eq. 14, we have coupled all of the terms from the standard Galerkin ROM on the left-hand side, and have similarly coupled the terms introduced by the Petrov–Galerkin projection on the right-hand side. One immediately observes that the LSPG approach is a residual-based method, meaning that the stabilization added by LSPG is proportional to the residual. The LSPG method is similar to the Galerkin/Least-Squares (GLS) approach commonly employed in the finite element community . This can be made apparent by writing Eq. 14 as,

where $(\mathbf{a},\mathbf{b})=\mathbf{a}^{T}\mathbf{b}$ and $\tau=\alpha\Delta t$ is the stabilization parameter. Compare the above to, say, Eq. 70 and 71 in Ref. . A rich body of literature exists on residual-based methods, and viewing the LSPG approach in this light helps establish connections with other methods. We highlight several important aspects of LSPG. Remarks 1 through 3 are derived by Carlberg et al. in Ref. :

The LSPG approach is inherently tied to the temporal discretization. For different time integration schemes, the “stabilization” added by the LSPG method will vary. For optimal accuracy, the LSPG method requires an intermediary time-step size.

In the limit of $\Delta t\rightarrow 0$, the LSPG approach recovers a Galerkin approach.

For explicit time integration schemes, the LSPG and Galerkin approach are equivalent.

For backwards differentiation schemes, the LSPG approach is a type of GLS stabilization for non-linear problems.

While commonalities exist between LSPG and multiscale approaches, the authors believe that the LSPG method should not be viewed as a subgrid-scale model. The reason for this is that it is unclear how Eq. 14 can be derived from Eq. 7. This is similar to the fact that, in Ref. , adjoint stabilization is viewed as a subgrid-scale model while GLS stabilization is not. The challenge in deriving Eq. 14 from Eq. 7 lies primarily in the fact that the Jacobian in Eq. 14 contains a transpose operator. We thus view LSPG as mathematical stabilization rather than a subgrid-scale model.

While the LSPG approach has enjoyed much success for constructing ROMs of non-linear problems, remarks 1, 2, 3, and 5 suggest that improvements over the LSPG method are possible. Remark 1 suggests improvements in computational speed and accuracy are possible by removing sensitivity to the time-step size. Remark 3 suggests that improvements in computational speed and flexibility are possible by formulating a method that can be used with explicit time-stepping schemes. Lastly, remark 5 suggests that improvements in accuracy are possible by formulating a method that accounts for subgrid effects.

3 Mori-Zwanzig Reduced-Order Models

The optimal prediction framework formulated by Chorin et al. , which is a (significant) reformulation of the Mori-Zwanzig (MZ) formalism of statistical mechanics, is a model order reduction tool that can be used to develop representations of the impact of the fine-scales on the coarse-scale dynamics. In this section, the optimal prediction framework is used to derive a compact approximation to the impact of the fine-scale POD modes on the evolution of the coarse-scale POD modes. For completeness, the optimal prediction framework is first derived in the context of the Galerkin POD ROM. It is emphasized that the content presented in Sections 3.3.1 and 3.3.2 is simply a formulation of Chorin’s framework, with a specific projection operator, in the context of the Galerkin POD ROM.

We pursue the MZ approach on a Galerkin formulation of Eq. 5. Before describing the formalism, it is beneficial to re-write the original FOM in terms of the generalized coordinates with the solution being defined implicitly as a function of the initial conditions,

for arbitrary q. Equation 16 is referred to as the Liouville equation and is an exact statement of the original dynamics. The Liouville equation describes the solution to Eq. 15 for all possible initial conditions. The advantage of reformulating the system in this way is that the Liouville equation is linear, allowing for the use of superposition and aiding in the removal of the fine-scales.

The solution to Eq. 16 can be written as,

The operator $e^{t\mathcal{L}}$, which has been referred to as a “propagator”, evolves the solution along its trajectory in phase-space . The operator $e^{t\mathcal{L}}$ has several interesting properties. Most notably, the operator can be “pulled” inside of a non-linear functional ,

This is similar to the composition property inherent to Koopman operators . With this property, the solution to Eq. 16 may be written as,

The implications of $e^{t\mathcal{L}}$ are significant. It demonstrates that, given trajectories $\mathbf{a}(\mathbf{a_{0}},t)$, the solution $v$ is known for any observable $g$.

Noting that $\mathcal{L}$ and $e^{t\mathcal{L}}$ commute, Eq. 16 may be written as,

3.2 Projection Operators and the Generalized Langevin Equation

The objective now is to remove the dependence of Eq. 17 on the fine-scale variables. Similar to the VMS decomposition, $\mathcal{R}^{N}$ can be decomposed into resolved and unresolved subspaces,

The projection operators can be used to split the Liouville equation,

Inserting Eq. 19 into Eq. 18, the generalized Langevin equation is obtained,

By the definition of the initial conditions (Eq. 11), the noise-term is zero and we obtain,

The system described in Eq. 20 is precise and not an approximation to the original ODE system. For notational purposes, define,

Note that the time derivative is represented as a partial derivative due to the Liouville operators embedded in the memory.

The derivation up to this point has cast the original full-order model in generalized coordinates (Eq. 15) as a linear PDE. Through the use of projection operators and Duhamel’s principle, an exact equation (Eq. 23) for the coarse-scale dynamics only in terms of the coarse-scale variables was then derived. The effect of the fine-scales on the coarse-scales appeared as a memory integral. This memory integral may be thought of as the closure term that is required to exactly account for the unresolved dynamics.

3.3 The τ𝜏\tau-model and the Adjoint Petrov–Galerkin Method

The direct evaluation of the memory term in Eq. 23 is, in general, computationally intractable. To gain a reduction in computational cost, an approximation to the memory must be devised. A variety of such approximations exist, and here we outline the $\tau$-model . The $\tau$-model can be interpreted as the result of assuming that the memory is driven to zero in finite time and approximating the integral with a quadrature rule. This can be written as a two-step approximation,

Equation 24 provides a closed equation for the evolution of the coarse-scales. The left-hand side of Eq. 24 is the standard Galerkin ROM, and the right-hand side can be viewed as a subgrid-scale model.

When compared to existing methods, the inclusion of the $\tau$-model leads to a method that is analogous to a non-linear formulation of the adjoint stabilization technique developed in the finite element community. The “adjoint” terminology arises from writing Eq. 24 in a Petrov–Galerkin form,

It is seen that Eq. 25 involves taking the inner product of the coarse-scale ODE with a test-basis that contains the adjoint of the coarse-scale Jacobian. Unlike GLS stabilization, adjoint stabilization can be derived from the multiscale equations . Due to the similarity of the proposed method with adjoint stabilization techniques, as well as the LSPG terminology, the complete ROM formulation will be referred to as the Adjoint Petrov–Galerkin (APG) method.

3.4 Comparison of APG and LSPG

The APG method displays similarities to LSPG. From Eq. 25, it is seen that the test basis for the APG ROM is given by,

Recall the LSPG test basis for backward differentiation schemes,

Analysis

This section presents theoretical analyses of the Adjoint Petrov–Galerkin method. Specifically, error and eigenvalue analyses are undertaken for linear time-invariant (LTI) systems. Section 4.1 derives a priori error bounds for the Galerkin and Adjoint Petrov–Galerkin ROMs. Conditions under which the APG ROM may be more accurate than the Galerkin ROM are discussed. Section 4.2 outlines the selection of the parameter $\tau$ that appears in APG.

where the Galerkin and Adjoint Petrov–Galerkin projections are, respectively,

The residual of the full-order model is defined as,

We define the error in the Galerkin and Adjoint Petrov–Galerkin method as,

Similarly, the coarse-scale error is defined as,

To simplify the analysis, the Adjoint Petrov–Galerkin projection is approximated to be stationary in time. Note that the Galerkin projection is stationary in time. For clarity, we suppress the temporal argument on the states when possible in the proofs.

A priori error bounds for the Galerkin and Adjoint Petrov–Galerkin ROMs are, respectively,

Invoking the assumption of Lipschitz continuity,

Noting that $\frac{d\left\lVert\mathbf{e}_{A}\right\rVert_{2}}{dt}\leq\left\lVert\frac{d\mathbf{e}_{A}}{dt}\right\rVert_{2}$ we have $\frac{d\left\lVert\mathbf{e}_{A}\right\rVert_{2}}{dt}=\frac{1}{\left\lVert\mathbf{e}_{A}\right\rVert_{2}}\mathbf{e}_{A}^{T}\frac{d\mathbf{e}_{A}}{dt}\leq\left\lVert\frac{1}{\left\lVert\mathbf{e}_{A}\right\rVert_{2}}\mathbf{e}_{A}\right\rVert_{2}\left\lVert\frac{d\mathbf{e}_{A}}{dt}\right\rVert_{2}\leq\left\lVert\frac{d\mathbf{e}_{A}}{dt}\right\rVert_{2}$ ,

An upper bound on the error for the Adjoint Petrov–Galerkin method is then obtained by solving Eq. 33 for $\left\lVert\mathbf{e}_{A}\right\rVert_{2}$, which yields,

Unfortunately, for general non-linear systems, a priori error analysis provides minimal insight beyond what was just mentioned. To obtain a more intuitive understanding of the APG method, error analysis in the case that $\mathbf{R}(\cdot)$ is a linear time-invariant operator is now considered.

Let $\mathbf{R}({\mathbf{u}})=\mathbf{A}{\mathbf{u}}$ be a linear time-invariant operator. Error bounds for the coarse-scales in the Galerkin and APG ROMs, are, respectively,

Equation 36 is a first-order linear inhomogeneous differential equation that can be solved analytically. A standard way to do this is by performing the eigendecomposition,

Using the definition of the full-order residual,

Let $\mathbf{R}({\mathbf{u}})=\mathbf{A}{\mathbf{u}}$ be a linear time-invariant operator. Upper bounds of the Galerkin and Adjoint Petrov-Galerkin projections of the full-order coarse-scale residual, are, respectively,

The result will be proved first for the Galerkin method and then for the Adjoint Petrov-Galerkin method. By definition of the residual and the Galerkin projector,

The evolution equation of the fine-scales is defined by the FOM projected onto the fine-scale space,

The above equation is a first-order nonhomogeneous linear system of differential equations. By the definition of the initial conditions, the solution to the fine-scales reduces to,

where $e^{{\Pi^{\prime}}\mathbf{A}}$ is the matrix exponential. Introducing a change of variables, $\zeta=t-s$, this becomes,

Breaking the integral up into two intervals and applying the triangle inequality, the desired result for Galerkin is obtained,

The first two terms on the right-hand side are the Galerkin term from Eq. 41. Substituting this into the above equation results in,

Theorem 3 provides an upper bound on the error introduced due to the closure problem for both the Galerkin and APG ROMs. These error bounds were obtained by analytically solving for the fine-scales as a function of coarse-scales. This led to a memory integral involving the past history of the coarse-scales. The memory integral expressing the solution to the fine-scales was then split into two intervals: one from $[0,\tau]$ and one from $[\tau,t]$. The APG method can be interpreted as approximating the integral over the first interval with a quadrature rule and ignoring the second interval. The Galerkin ROM ignores both intervals. Next, Theorem 4 shows that, for sufficiently small $\tau$, the APG approximation is more accurate than the Galerkin approximation. This result is intuitive since, for sufficiently small $\tau$, the quadrature rule used in deriving APG is accurate.

In the limit of $\tau\rightarrow 0^{+}$ the error bound provided in Theorem 3 is less per unit $\tau$ in APG than in Galerkin,

Define the difference between the APG and Galerkin error incurred per unit $\tau$ as,

After canceling terms, we equivalently have

Apply the rectangle rule to alternatively write the integrals in the above as,

for some $c$ in the interval $[0,\tau$]. Inserting Eq. 45 into Eq. 44,

Factoring out $\tau$, which is non-negative by definition,

Taking the limit as $\tau\rightarrow 0^{+}$ and noting the norm is always non-negative,

From Eq. 43 and the algebraic limit theorem, this implies that,

Theorem 4 shows that, for sufficiently small $\tau$, an upper bound on the a priori error introduced due to the closure problem in APG is less than in Galerkin. This suggests that APG may be more accurate than Galerkin. The results of Theorem 4 will be discussed more in Section 4.1.1. While Theorem 4 shows that, for sufficiently small $\tau$, APG may be more accurate than Galerkin, it does not give insight into appropriate values of $\tau$. Theorem 5 and Corollary 5.1 seek to provide insight into this issue.

Expanding the Adjoint Petrov-Galerkin operator,

Noting that the first term on the right-hand side is just the Galerkin term,

We again expand the APG projection and simplify,

To ensure that $\lambda_{A_{1}}\leq 0$, we need only show under what conditions $\kappa_{1}\leq 0$. Thus, we seek bounds on $\tau$ such that $\kappa_{1}\leq 0$.

As $\mathbf{A}$ and $\tau\mathbf{A}{\Pi^{\prime}}\mathbf{A}$ are both self-adjoint, we can use the Weyl inequalities to bound the smallest eigenvalue of $\mathbf{A}+\tau\mathbf{A}{\Pi^{\prime}}\mathbf{A}$. Let $\chi_{1}\geq\chi_{2}\geq\ldots\geq\chi_{N}$ be the $N$ eigenvalues of $\mathbf{A}{\Pi^{\prime}}\mathbf{A}$, thus,

To ensure that $\kappa_{1}\leq 0$, we specify that $\lambda_{1}+\tau\chi_{1}\leq 0$. Noting that $\gamma_{1}\leq 0$ leads to the inequality,

To write the bounds on $\tau$ in Eq. 49 into a more intuitive form, the $\chi_{i}$ can be related to $\lambda_{i}$ through the the Weyl inequalities. Recall that $\chi_{i}$ are the eigenvalues of,

Inserting the above inequality into Eq. 49 obtains the desired result,

This upper bound on $\tau$ is more conservative than Eq. 49, though both are equally valid.

Theorems 4 and 5 contain three interesting results that are worth discussing. First, as discussed in Corollary 2.1, Theorem 4 shows that, in the limit $\tau\rightarrow 0^{+}$, the upper-bound provided in Theorem 3 on the error introduced at time $t$ in the APG ROM due to the fine-scales is less than that introduced in the Galerkin ROM (per unit $\tau$). This is an appealing result as APG is derived as a subgrid-scale model. Two remarks are worth making regarding this result. First, while the result was demonstrated in the limit $\tau\rightarrow 0^{+}$, it will hold so long as the norm of the truncation error in the quadrature approximation is less than the norm of the integral it is approximating; i.e. the approximation is doing a better job than neglecting the integral entirely. Second, the result derived in Theorem 4 does not directly translate to showing that,

This is a consequence of Theorem 3, in which the integral that defines the fine-scale solution was split into two intervals. The APG ROM attempts to approximate the first integral, while the Galerkin ROM ignores both terms. Theorem 4 showed that, in the limit $\tau\rightarrow 0^{+}$, the APG approximation to the first integral is better than in the case of Galerkin (i.e., the APG approximation is better than no approximation). The only time that the result provided in Theorem 4 will not translate to APG providing a better approximation to the entire integral (and thus proving Eq. 51), is when integration over the second interval “cancels out” the integration over the first interval. To make this idea concrete, consider the integral,

Clearly, $\int_{0}^{2\pi}\sin(x)dx=0$. If the entire integral was approximated using just the interval to $\pi$ (which is analogous to APG), then one would end up with the approximation $\int_{0}^{2\pi}\sin(x)dx\approx\int_{0}^{\pi}\sin(x)dx=2$. Alternatively, if one were to ignore the integral entirely (which is analogous to Galerkin) and make the approximation $\int_{0}^{2\pi}\sin(x)dx\approx 0$ (which in this example is exact), a better approximation would be obtained.

The next interesting result is presented in Theorem 5, where it is shown that for a self-adjoint system with negative eigenvalues, the eigenvalues associated with the APG ROM error equation are greater than the Galerkin ROM. This implies that APG is less dissipative than Galerkin, and means that errors may be slower to decay in time. Thus, while Theorem 4 shows that the a priori contributions to the error due to the closure problem may be smaller in the APG ROM than in the Galerkin ROM, the errors that are incurred may be slower to decay. Finally, Corollary 5.1 shows that, for self-adjoint systems with negative eigenvalues, the bounds on the parameter $\tau$ such that all eigenvalues associated with the evolution of the error in APG remain negative depends on the spectral content of the Jacobian of $\mathbf{A}$. This observation has been made heuristically in Ref . Although the upper bound on $\tau$ in Eq. 50 is very conservative due to the repeated use of inequalities, it provides insight into the selection and behaviour of $\tau$.

2 Selection of Memory Length τ𝜏\tau

The APG method requires the specification of the parameter $\tau$. Theorem 5 showed that, for a self-adjoint linear system, bounds on the value of $\tau$ are related to the eigenvalues of the Jacobian of the full-dimensional right-hand side operator. While such bounds provide intuition into the behavior of $\tau$, they are not particularly useful in the selection of an optimal value of $\tau$ as they 1.) are conservative due to repeated use of inequalities and 2.) require the eigenvalues of the full right-hand side operator, which one does not have access to in a ROM. Further, the bounds were derived for a self-adjoint linear system, and the extension to non-linear systems is unclear.

where $C$ is a model parameter and $\rho(\cdot)$ indicates the spectral radius. In Ref. , $C$ was reported to be $0.2.$ In the numerical experiments presented later in this manuscript, the sensitivity of APG to the value of $\tau$ and the validity of Eq. 52 are examined.

Similar to the selection of $\tau$ in the APG method, the LSPG method requires the selection of an appropriate time-step . In practice, this fact can be problematic as finding an optimal time-step for LSPG which minimizes error may result in a small time-step and, hence, an expensive simulation. The selection of the parameter $\tau$, on the other hand, does not impact the computational cost of the APG ROM.

Implementation and Computational Cost of the Adjoint Petrov–Galerkin Method

This section explores the cost of the APG method within the scope of explicit time integration schemes. For simplicity, the analysis is carried out only for the explicit Euler scheme. The computational cost of more sophisticated time integration methods, such as Runge-Kutta and multistep schemes, is generally a proportional scaling of the cost of the explicit Euler scheme. Algorithm 1 provides the step-by-step procedure for performing an explicit Euler update to the Adjoint Petrov–Galerkin ROM. Table 1 provides the approximate floating-point operations for the steps reported in Algorithm 1. The algorithm for an explicit update to the Galerkin ROM, along with the associated FLOP counts, is provided in Algorithm 7 and Table 8 in Appendix C. As noted previously, LSPG reverts to the Galerkin method for explicit schemes, and so is not detailed in this section. Table 1 shows that, in the case that $K\ll N$ (standard for a ROM) and $\omega\gg 1$ (sufficiently complex right-hand side), the Adjoint Petrov–Galerkin ROM is approximately twice as expensive as the Galerkin ROM.

2 Implicit Time Integration Schemes

This section evaluates the computational cost of the Galerkin, Adjoint Petrov–Galerkin, and Least-Squares Petrov–Galerkin methods for implicit time integration schemes. For non-linear systems, implicit time integration schemes require the solution of a non-linear algebraic system at each time-step. Newton’s method, along with a preferred linear solver, is typically employed to solve the system. For simplicity, the analysis provided in this section is carried out for the implicit Euler time integration scheme along with Newton’s method to solve the non-linear system. Before proceeding, the full-order residual, Galerkin residual, and APG residual at time-step $(n+1)$ are denoted as,

Two methods are considered for the solution to the non-linear algebraic system arising from implicit time discretizations of the G and APG ROMs: Newton’s method with direct Gaussian elimination and Jacobian-Free Newton-Krylov GMRES. The Gauss-Newton method with Gaussian elimination is considered for the solution to the least-squares problem arising in LSPG.

Algorithm 2 provides the step-by-step procedures for performing an implicit Euler update to the Adjoint Petrov–Galerkin ROM with the use of Newton’s method and Gaussian elimination. Table 2 provides the approximate floating-point operations for the steps reported in these algorithms. Analogous results for the Galerkin and LSPG ROMs are reported in Algorithms 8 and 9, and Tables 9 and 10 in Appendix C. In the limit that $K\ll N$ and $\omega\gg 1$, the total FLOP counts reported show that APG is twice as expensive as both the LSPG and Galerkin ROMs. It is observed that the dominant cost for all three methods lies in the computation of the low-dimensional residual Jacobian. Computation of the low-dimensional Jacobian requires $K$ evaluations of the unsteady residual. Depending on values of $\omega$, $N$, and $K$, this step can consist of over $50\%$ $(K\ll N)$ of the CPU time.It is noted that the low-dimensional Jacobian can be computed in parallel.

The LSPG method is formulated as a non-linear least-squares problem. The use of Jacobian-free methods to solve non-linear least-squares problems is significantly more challenging. The principle issue encountered in attempting to use Jacobian-free methods for such applications it that one requires the action of the transpose of the residual Jacobian on a vector. This quantity cannot be computed via a standard finite difference approximation or linearization. It is only recently that true Jacobian-free methods have been utilized for solving non-linear least-squares problems. In Ref , for example, automatic differentiation is utilized to compute the action of the transposed Jacobian on a vector. Due to the challenges associated with Jacobian-free methods for non-linear least-squares problems, this method is not considered here as a solution technique for LSPG.

Algorithm 3 and Table 3 report the algorithm and FLOPs required for an implicit Euler update to APG using JFNK GMRES. The term $\eta\leq K$ is the number of iterations needed for convergence of the GMRES solver at each Newton iteration. For a concise presentation, the same update for the Galerkin ROM is not presented. Figure 1 shows the ratio of the cost of the various implicit ROMs as compared to the Galerkin ROM solved with Gaussian elimination. The standard LSPG method is seen to be approximately the same cost of Galerkin, while APG is seen to be approximately $2$x the cost of Galerkin. The success of the JFNK methods depends on the number of GMRES iterations required for convergence. If $\eta=K$, which is the maximum number of iterations required for GMRES, the cost of JFNK methods is seen to be the same as their direct-solve counterparts. For cases where JFNK converges at a rate of $\eta<K,$ the iterative methods out-perform their direct-solve counterparts.

The analysis presented here shows that, for a given basis dimension, the Adjoint Petrov–Galerkin ROM is approximately twice the cost of the Galerkin ROM for both implicit and explicit solvers. In the implicit case, the APG ROM utilizing a direct linear solver is approximately 2x the cost of LSPG. It was highlighted, however, that APG can be solved via JFNK methods. For cases where one either doesn’t have access to the full Jacobian, or the full Jacobian can’t be stored, JFNK methods can significantly decrease the ROM cost. The use of JFNK methods within the LSPG approach is more challenging due to the presence of the transpose of the residual Jacobian. Lastly it is noted that, although hyper-reduction can decrease the cost of a residual evaluation, it does not entirely alleviate the cost of forming the Jacobian.

Numerical Examples

Applications of the APG method are presented for ROMs of compressible flows: the 1D Sod shock tube problem and 2D viscous flow over a cylinder. In both problems, the test bases are chosen via POD. The shock tube problem highlights the improved stability and accuracy of the APG method over the standard Galerkin ROM, as well as improved performance over the LSPG method. The impact of the choice of $\tau$ (APG) and $\Delta t$ (LSPG) time-scales are also explored. The cylinder flow experiment examines a more complex problem and assesses the predictive capability of APG in comparison with Galerkin and LSPG ROMs. The effect of the choice of $\tau$ on simulation accuracy is further explored.

The first case considered is the Sod shock tube, described in more detail in . The experiment simulates the instantaneous bursting of a diaphragm separating a closed chamber of high-density, high-pressure gas from a closed chamber of low-density, low pressure gas. This generates a strong shock, a contact discontinuity, and an expansion wave, which reflect off the shock tube walls at either end and interact with each other in complex ways. The system is described by the one-dimensional compressible Euler equations with the initial conditions,

with $x\in$. Impermeable wall boundary conditions are enforced at x = 0 and x = 1.

The 1D compressible Euler equations are solved using a finite volume method and explicit time integration. The domain is partitioned into 1,000 cells of uniform width. The finite volume method uses the first-order Roe flux at the cell interfaces. A strong stability-preserving RK3 scheme is used for time integration. The solution is evolved for $t\in[0.0,1.0]$ with a time-step of $\Delta t=0.0005$, ensuring CFL$\leq 0.75$ for the duration of the simulation. The solution is saved every other time-step, resulting in 1,000 solutions snapshots for each conserved variable.

1.2 Solution of the Reduced-Order Model

Remark: The Adjoint Petrov–Galerkin ROM requires specification of $\tau$.

Least-Squares Petrov–Galerkin ROM (Implicit Euler Time Integration):

Remark: The LSPG approach is strictly coupled to the time integration scheme and time-step.

1.3 Numerical Results

The first case considered uses $50$ basis vectors each for the conserved variables $\rho,\rho u,$ and $\rho E$. The total dimension of the reduced model is thus $K=150$. Roughly 99.9–99.99% of the POD energy is captured by this 150-mode basis. In fact, 99% of the energy is contained in the first 5-12 modes of each conserved variable.

The Adjoint Petrov–Galerkin ROM requires specification of the memory length $\tau$. Similarly, LSPG requires the selection of an appropriate time-step. The sensitivity of both methods to this selection will be discussed later in this section. The simulation parameters are provided in Table 4.

Density profiles at $t=0.25$ and $t=1.0$ for explicit Galerkin and APG ROMs, along with an implicit LSPG ROM, are displayed in Fig. 2. All three ROMs are capable of reproducing the shock tube density profile in Fig. 2(a); a normal shock propagates to the right and is followed closely behind by a contact discontinuity, while an expansion wave propagates to the left. All three methods exhibit oscillations at $x=0.5$, the location of the imaginary burst diaphragm, and near the shock at $x=0.95$. At $t=1.0$, when the shock has reflected from the right wall and interacted with the contact discontinuity, much stronger oscillations are present, particularly near the reflected shock at $x=0.45$. These oscillations are reminiscent of Gibbs phenomenon, and are an indicator of the inability to accurately reconstruct sharp gradients. The Galerkin ROM exhibits the largest oscillations of the ROMs considered, while LSPG exhibits the smallest.

Figure 3 shows the evolution of the error for all of the ROMs listed in Table 4. The $L^{2}$-norm of the error is computed as,

In Figure 3, it is seen that the APG ROM exhibits improved accuracy over the Galerkin ROM. The LSPG ROM for $\Delta t=0.0005$ performs slightly better than the explicit Galerkin ROM, and worse than the implicit Galerkin ROM. Increasing the time-step to $\Delta t=0.001$ results in a significant increase in error for the LSPG ROM. This is due to the fact that the performance of LSPG is influenced by the time-step. For a trial basis containing much of the residual POD energy, LSPG will generally require a very small time-step to improve accuracy; this sensitivity will be explored later. Lastly, it is observed that the APG ROM is not significantly affected by the time-step. The APG ROM with $\Delta t=0.001$ shows moderately increased error prior to $t=0.3$ and similar error afterwards when compared against the $\Delta t=0.0005$ APG ROM case.

Figure 4 studies the effect of the number of modes retained in the trial basis on the stability and accuracy, over the range $K=60,75,90,\ldots,180$. Missing data points indicate an unstable solution. Values of $\tau$ for the APG ROMs are again selected by user choice. The most striking feature of these plots is the fact that even though the explicit Galerkin ROM is unstable for $K\leq 135$ and the implicit Galerkin ROM is unstable for $K\leq 75$, the APG and LSPG ROMs are stable for all cases. Furthermore, the APG and LSPG ROMs are capable of achieving stability with a time-step twice as large as that of the Galerkin ROM. The cost of the APG and LSPG ROMs are effectively halved, but they are still able to stabilize the simulation. Interestingly, the Galerkin and APG ROMs both exhibit abrupt peaks in error at $K=120$, while the LSPG ROMs do not. The exact cause of this is unknown, but displays that a monotonic decrease in error with enrichment of the trial space is not guaranteed.

Several interesting comparisons between APG and LSPG arise from Figure 4. First, with the exception of the $K=120$ case, Fig. 3(a) shows that the APG ROM with explicit time integration exhibits accuracy comparable to that of the LSPG ROM with implicit time integration. As can be seen in comparing Tables 1 and 10, the cost of APG with explicit time integration is significantly lower than the cost of LSPG. This is an attractive feature of APG, as it is able to use inexpensive explicit time integration while LSPG is restricted to implicit methods. Additionally, we draw attention to the poor performance of LSPG at high $K$ for a moderate time-step in Fig. 4(b). Increasing the time-step to $\Delta t=0.001$ to decrease simulation cost only exacerbates this issue; as the trial space is enriched, LSPG requires a smaller time-step to yield accurate results. If we wish to improve the LSPG solution for $K=150$, we must decrease the time-step below that of the FOM. The accuracy of the APG ROM does not change when the time-step is doubled from $\Delta t=0.0005$ to $\Delta t=0.001$. This halves the cost of the APG ROM with no significant drawbacks.

1.4 Optimal Memory Length Investigations

As mentioned previously, the success of LSPG is tied to the physical time-step and the time integration scheme, while the parameter $\tau$ in the APG method may be chosen independently from these factors. In minimizing ROM error, finding an optimal value of $\tau$ for the APG ROM may permit the choice of a much larger time-step than the optimal LSPG time-step. Further, the APG method may be applied with explicit time integration schemes, which are generally much less expensive than the implicit methods which LSPG is restricted to. To demonstrate this, the APG ROM and LSPG ROM with $K=150$ are simulated for a variety of time scales ($\tau$ for APG and $\Delta t$ for LSPG).

Figure 6 shows the integrated error of the ROMs versus the relevant time scale. For this case, the optimal value of $\Delta t$ for LSPG is less than $0.0001$ and is not shown. The optimal value of $\tau$ is not greatly affected by the choice of time integration scheme (implicit or explicit) or time-step. Furthermore, because $\tau$ can be chosen independently from $\Delta t$ for APG, the APG ROM can produce low error at a much larger time-step ($\Delta t=0.001$) than the optimal time-step for the LSPG ROM. This highlights the fact that the “optimal” LSPG ROM may be computationally expensive due to a small time-step, whereas the “optimal” APG ROM requires only the specification of $\tau$ and can use, potentially, much larger time-steps than LSPG. It has to be mentioned, however, that the choice of $\tau$ has an impact on performance — selections of $\tau$ larger than those plotted in Fig. 6 caused the ROM to lose stability.

The spectral radius plays an important role in both implicit and explicit time integrators and is often the determining factor in the choice of the time-step. Theoretical analysis on the stability of explicit methods (and convergence of implicit methods) shows a similar dependence to the spectral radius. Choosing the memory length to be $\tau=\Delta t$ is one simple heuristic that may be used.

While a linear relationship between $\tau$ and the spectral radius of the coarse-scale Jacobian has been observed in every problem the authors have examined, the slope of the fit is somewhat problem dependent. For the purpose of reduced-order modeling, however, this is only a minor inconvenience as an appropriate value of $\tau$ can be selected by assessing the performance of the ROM on the training set, i.e., on the simulation used to construct the POD basis.

Finally, more complex methods may be used to compute $\tau$. A method to dynamically compute $\tau$ based on Germano’s identity, for instance, was proposed in in the context of the simulation of turbulent flows with Fourier-Galerkin methods. Extension of this technique to projection-based ROMs and the development of additional techniques to select $\tau$ will be the subject of future work.

2 Example 2: Flow Over Cylinder

The second case considered is viscous compressible flow over a circular cylinder. The flow is described by the two-dimensional compressible Navier-Stokes equations. A Newtonian fluid and a calorically perfect gas are assumed.

The compressible Navier-Stokes equations are solved using a discontinuous Galerkin (DG) method and explicit time integration. Spatial discretization with the discontinuous Galerkin method leads to a semi-discrete system of the form,

For the flow over cylinder problem considered in this section, a single block domain is constructed in polar coordinates by uniformly discretizing in $\theta$ and by discretizing in the radial direction by,

where $R_{g}$ is a stretching factor and is defined by,

The DG method utilizes the Roe flux at the cell interfaces and uses the first form of Bassi and Rebay for the viscous fluxes. Temporal integration is again performed using a strong stability preserving RK3 method. Far-field boundary conditions and linear elements are used. Details of the FOM are presented in Table 5.

2.2 Solution of the Full-Order Model and Construction of the ROM Trial Space

Flow over a cylinder at Re=$100,200,$ and $300$, where Re=$\rho_{\infty}U_{\infty}D/\mu$ is the Reynolds number, are considered. These Reynolds numbers give rise to the well studied von Kármán vortex street. Figure 7 shows the FOM solution at Re=100 for several time instances to illustrate the vortex street.

The FOM is used to construct the trial spaces used in the ROM simulations. The process used to construct these trial spaces is as follows:

Initialize FOM simulations at Reynold’s numbers of Re=$100,200,$ and $300.$ The Reynold’s number is controlled by raising or lowering the viscosity.

Time-integrate the FOM at each Reynolds number until a statistically steady-state is reached.

Once the flow has statistically converged to a steady state, reset the time coordinate to be $t=0$, and solve the FOM for $t\in$.

Take snapshots of the FOM solution obtained from Step 3 at every $t=0.5$ time units over a time-window of $t\in$ time units, for a total of $200$ snapshots at each Reynolds number. This time window corresponds to roughly two cycles of the vortex street, with $100$ snapshots per cycle.

Assemble the snapshots from each case into one global snapshot matrix of dimension $N\times 600$. This snapshot matrix is used to construct the trial subspace through POD. Note that only one set of basis functions for all conserved variables is constructed.

Construct trial spaces of dimension $N\times 11$, $N\times 43$, and $N\times 87$. These subspace dimensions correspond to an energy criterion of $99,99.9,$ and $99.99\%$. The different trial spaces are summarized in Table 6.

2.3 Solution of the Reduced-Order Models

The G ROM, APG ROM, and LSPG ROMs are considered. Details on their implementation are as follows:

Galerkin ROM: The Galerkin ROM is evolved in time using both explicit and implicit time integrators. In the explicit case, a strong stability RK3 method is used. In the implicit case, Crank-Nicolson time integration is used. The non-linear algebraic system is solved using SciPy’s Jacobian-Free Netwon-Krylov solver. LGMRES is employed as the linear solver. The convergence tolerance for the max-norm of the residual is set at the default ftol=6e-6.

LSPG ROM: The LSPG ROM is formulated from an implicit Crank-Nicolson temporal discretization. The resulting non-linear least-squares problem is solved using SciPy’s least-squares solver with the ‘dogbox’ method . The tolerance on the change to the cost function is set at ftol=1e-8. The tolerance on the change to the generalized coordinates is set at xtol=1e-8. The SciPy least-squares solver is comparable in speed to our own least-squares solver that utilizes the Gauss-Newton method with a thin QR factorization to solve the least-squares problem. The SciPy solver, however, was observed to be more robust in driving down the residual than the basic Gauss-Newton method with QR factorization, presumably due to SciPy’s inclusion of trust-regions, and hence results are reported with the SciPy solver.

All ROMs are initialized with the solution of the Re=$100$ FOM at time $t=0$, the $x$-velocity of which is shown in Figure 7(a).

2.4 Reconstruction of Re=100 Case

Reduced-order models of the Re=100 case are first considered. This case was explicitly used in the construction of the POD basis and tests the ability of the ROM to reconstruct previously “seen” dynamics. Unless otherwise noted, the default time-step for all ROMs is taken to be $\Delta t=0.5$. The values of $\tau$ used in the APG ROMs are selected from the spectral radius heuristic and are given in Table 6. Figures 8(a) and 8(b) show the lift coefficient as well as the mean squared error (MSE) of the full-field ROM solutions for the G ROM, APG ROM, and LSPG ROMs for Basis #2, while Figure 8(c) shows the integrated MSE for $t\in$ for Basis #1, 2, and 3. Figure 8(d) shows the integrated error as a function of relative CPU time for the various ROMs. The relative CPU time is defined with respect to the FOM, which is integrated with an explicit time-step 100 times lower than the ROMs. The lift coefficients predicted by all three ROMs are seen to overlay the FOM results. The mean squared error shows that, for a given trial basis dimension, the APG ROM is more accurate than both the Galerkin and LSPG ROMs. This is the case for both explicit and implicit time integrators. As shown in Figure 8(c), the APG ROM converges at a similar rate to the G ROM as the dimension of the trial space grows. For Basis #2 and #3, the implicit time-marching schemes are slightly less accurate than the explicit time-marching schemes. Finally, Figure 8(d) shows that, for a given CPU time, the G ROM with explicit time-marching produces the least error. The APG ROM with explicit time-marching is the second-best performing method. In the implicit case, both the G and APG ROMs lead to lower error at a given CPU time than LSPG. This decrease in cost is due to the fact that the G and APG ROMs utilize Jacobian-Free Netwon-Krylov solvers. As discussed in Section 5, it is much more challenging for LSPG to utilize Jacobian-Free methods. Due to the increased cost associated with implicit solvers, only explicit time integration is used for the G ROM and APG ROM beyond this point.

Next, we investigate the sensitivity of the different ROMs to the time-step size. Reduced-order models of the Re=$100$ case using Basis #2 are solved using time-steps of \Delta t=\big{[}0.1,0.2,0.5,1]. Note that the largest time-step considered is $200$ times larger than the FOM time-step, thus reducing the temporal dimensionality of the problem by 200 times. The mean-squared error of each ROM solution is shown in Figure 9. The G and APG ROMs are stable for all time-steps considered. Further, it is seen that varying the time-step has a minimal effect on the accuracy of the G and APG ROMs. In contrast, the accuracy of LSPG deteriorates if the time-step grows too large. This is due to the fact that, as shown in Ref. , the stabilization added by LSPG depends on the time-step size. Optimal accuracy of the LSPG method requires an intermediate time-step. The ability of the APG and G ROMs to take large time-steps without a significant degradation in accuracy allows for significant computational savings. This advantage is further amplified when large time-steps can be taken with an explicit solver, as is the case here. This will be discussed in more detail in Section 6.2.6.

Lastly, we numerically investigate the sensitivity of APG to the parameter $\tau$ by running simulations for $\tau=[0.001,0.01,0.1,0.3,0.5,1.]$. All simulations are run at $\Delta t=0.5$. The results of the simulations are shown in Figure 10. It is seen that, for all values of $\tau$, the APG ROM produces a better solution than the G ROM. The lowest error is observed for an intermediate value of $\tau$, in which case the APG ROM leads to over a $50\%$ reduction in error from the G ROM. As $\tau$ approaches zero, the APG ROM solution approaches the Galerkin ROM solution. Convergence plots for LSPG as a function of $\Delta t$ are additionally shown in Figure 10(b). It is seen that the optimal time-step in LSPG is similar to the optimal value of $\tau$ in APG.

2.5 Parametric Study of Reynolds Number Dependence

Next, the ability of the different ROMs to interpolate between between different Reynolds numbers is studied. Simulations at Reynolds numbers of Re=100,150,200,250, and 300 with Basis $\#2$ and $\#3$ are performed. All cases are initialized from the Re=100 simulation. Note that the trial spaces in the ROMs were constructed from statistically steady-state FOM simulations of Re=100,200,300. The Reynolds number is modified by changing the viscosity.

Figure 11 summarizes the amplitude of the lift coefficient signal as well as the shedding frequency for the various methods. The values reported in Figure 11 are computed from the last 150s of the simulationsNot all G ROMs reached a statistically steady state over the time window considered. The Galerkin ROM is seen to do poorly in predicting the lift coefficient amplitude for both Basis $\#2$ and Basis $\#3$. Unlike in the Re=100 case, enhancing the basis dimension does not improve the performance of the ROMs. Both the LSPG and APG ROMs are seen to offer much improved predictions over the Galerkin and ROM. This result is promising, as the ultimate goal of reduced-order modeling is to provide predictions in new regimes.

The results presented in this example highlight the shortcomings of the Galerkin ROM. To obtain results that are even qualitatively correct for the Re={150,200,250,300} cases, either APG or LSPG must be used. As reported in Figure 8(d), explicit APG is over an order of magnitude faster than LSPG, and implicit APG with a JFNK solver is anywhere from 2x to 5x faster than LSPG. Therefore, APG is the best-performing method for this example.

2.6 Flow Over Cylinder at Re=100100100 with Hyper-Reduction

The last example considered is again flow over a cylinder, but this time the reduced-order models are augmented with hyper-reduction. The purpose of this example is to examine the performance of the different reduced-order models when fully equipped with state-of-the-art reduction techniques. For hyper-reduction, an additional snapshot matrix of the right-hand side is generated. This additional snapshot matrix is generated by following steps one through five provided in Section 6.2.2. Hyper-reduction for the G ROM and the APG ROM is achieved through the Gappy POD method . Hyper-reduction for LSPG is achieved through collocation using the same sampling points.It is noted that collocated LSPG out-performed the GNAT method for this example, and thus GNAT is not considered. When augmented with hyper-reduction, the trial basis dimension ($K$), right-hand side basis dimension ($r$), and number of sample points $(N_{s})$ can impact the performance of the ROMs. Table 7 summarizes the various permutations of $K$, $r$, and $N_{s}$ considered in this example. The sample points are selected through a QR factorization of the right-hand side snapshot matrix . These sample points are then augmented such that they contain every conserved variable and quadrature point at the selected cells. The sample mesh corresponding to Basis numbers 4,5, and 6 in Table 7 is shown in Figure 12. Details on hyper-reduction and its implementation in our discontinuous Galerkin code are provided in Appendix A.

Conclusion

This work introduced the Adjoint Petrov–Galerkin method for non-linear model reduction. Derived from the variational multiscale method and Mori-Zwanzig formalism, the Adjoint Petrov–Galerkin method is a Petrov–Galerkin projection technique with a non-linear time-varying test basis. The method is designed to be applied at the semi-discrete level, i.e., after spatial discretization of a partial differential equation, and is compatible with both implicit and explicit time integration schemes. The method displays commonalities with the adjoint-stabilization method used in the finite element community as well as the Least-Squares Petrov–Galerkin approach used in non-linear model-order reduction. Theoretical error analysis was presented that showed conditions under which the Adjoint Petrov–Galerkin ROM may have lower a priori error bounds than the Galerkin ROM. The theoretical cost of the Adjoint Petrov–Galerkin method was considered for both explicit and implicit schemes, where it was shown to be approximately twice that of the Galerkin method. In the case of implicit time integration schemes, the Adjoint Petrov–Galerkin ROM was shown to be capable of being more efficient than Least-Squares Petrov–Galerkin when the non-linear system is solved via Jacobian-Free Newton-Krylov methods.

Numerical experiments with the Adjoint Petrov–Galerkin, Galerkin, and Least-Squares Petrov–Galerkin method were presented for the Sod shock tube problem and viscous compressible flow over a cylinder parameterized by the Reynolds number. In all examples, the Adjoint Petrov–Galerkin method provided more accurate predictions than the Galerkin method for a fixed basis dimension. Improvements over the Least-Squares Petrov–Galerkin method were observed in most cases. In particular, the Adjoint Petrov–Galerkin method was shown to provide relatively accurate predictions for the cylinder flow at Reynolds numbers outside of the training set used to construct the POD basis. The Galerkin method, with both an equivalent and an enriched trial space, failed to produce accurate results in these cases. Additionally, numerical evidence showed a correlation between the spectral radius of the reduced Jacobian and the optimal value of the stabilization parameter appearing in the Adjoint Petrov–Galerkin method.

When augmented with hyper-reduction, the Adjoint Petrov–Galerkin ROM was shown to be capable of producing accurate predictions within the POD training set with computational speedups up to 5000 times compared to the full-order models. This speed-up is a result of hyper-reduction of the right-hand side, as well as the ability to use explicit time integration schemes at large time-steps. A study of the Pareto front for simulation error versus relative wall time showed that, for the compressible cylinder problem, the Adjoint Petrov–Galerkin ROM is competitive with the Galerkin ROM, and more efficient than the LSPG ROM for the problems considered.

Acknowledgements

The authors acknowledge support from the US Air Force Office of Scientific Research through the Center of Excellence Grant FA9550-17-1-0195 (Tech. Monitors: Mitat Birkan & Fariba Fahroo) and the project LES Modeling of Non-local effects using Statistical Coarse-graining (Tech. Monitors: Jean-Luc Cambier & Fariba Fahroo). E. Parish acknowledges an appointment to the Sandia National Laboratories John von Neumann fellowship. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525.

Appendix A Hyper-reduction for the Adjoint Petrov–Galerkin Reduced-Order Model

The gappy POD method seeks to find an approximation for $\mathbf{R}(\cdot)$ that evaluates the right-hand side term at a reduced number of spatial points $r\ll N$. This is achieved through the construction of a trial space for the right-hand side and least-squares reconstruction of a sampled signal. The offline steps required in the Gappy POD method are given in Algorithm 4.

A.2 Selection of Sampling Points

Step 3 in Algorithm 4 requires the construction of the sampling point matrix, for which several methods exist. In the discrete interpolation method proposed by Chaturantabut and Sorensen , the sample points are selected inductively from the basis $\mathbf{U}$, based on an error measure between the basis vectors and approximations of the basis vectors via interpolation. The method proposed by Drmac and Gugercin leverages the rank-revealing QR factorization to compute $\mathbf{P}$. Dynamic updates of the $\mathbf{U}$ and $\mathbf{P}$ via periodic sampling of the full-order right-hand side term is even possible via the methods developed by Peherstorfer and Willcox .

The 2D compressible cylinder simulations presented in this manuscript uses a modified version of the rank-revealing QR factorization proposed in to obtain the sampling points. The modifications are added to enhance the stability and accuracy of the hyper-reduced ROM within the discontinuous Galerkin method. Algorithm 5 outlines the steps used in this manuscript to compute the sampling points.

A.3 Hyper-Reduction of the Adjoint Petrov–Galerkin ROM

Lastly, the online steps required for an explicit Euler update to the Adjoint Petrov–Galerkin ROM with Gappy POD hyper-reduction is provided in Algorithm 6. It is worth noting that Step 3 in Algorithm 6 requires one to reconstruct the right-hand side at the stencil points. Hyper-reduction via a standard collocation method, which provides no means to reconstruct the right-hand side, is thus not compatable with the Adjoint Petrov–Galerkin ROM.

Appendix B POD Basis Construction

For the 1D Euler case detailed in this manuscript, the procedure for constructing separate POD bases for each conserved variables is as follows:

Run the full-order model for $t\in(0,1)$ at a time-step of $\Delta t=0.0005$. The state vector is saved at every other time step to create 1000 state snapshots.

Collect the snapshots for each state into three state snapshot matrices:

Compute the singular-value decomposition (SVD) of each snapshot matrix, e.g. for $\mathbf{S}_{\rho}$,

The columns of $\mathbf{V}_{\rho}$ and $\mathbf{U}_{\rho}$ are the left and right singular vectors of $\mathbf{S}_{\rho}$, respectively. ${\Sigma}_{\rho}$ is a diagonal matrix of the singular values of $\mathbf{S}_{\rho}$. The columns of $\mathbf{V}_{\rho}$ form a basis for the solution space of $\rho_{i}$.

In this example, a separate basis is computed for each conserved quantity. It is also possible to construct a global basis by stacking $\mathbf{S}_{\rho}$, $\mathbf{S}_{\rho\mathbf{u}}$, and $\mathbf{S}_{\rho\mathbf{E}}$ into one snapshot matrix and computing one “global” SVD.

Decompose each basis into bases for the resolved and unresolved scales by selecting the first $K$ columns and last $1000-K$ columns, respectively, e.g.

Each basis vector is orthogonal to the others, hence the coarse and fine-scales are orthogonal.

In this example, we have selected 1000 snapshots such that the column space of ${\mathbf{V}}$ spans ${\mathcal{V}}$. In general, this is not the case. As the APG method requires no processing of the fine-scale basis functions, however, this is not an issue.

Appendix C Algorithms for the Galerkin and LSPG ROMs

Section 5 presented an analysis on the computational cost of the Adjoint Petrov–Galerkin ROM. This appendix presents similar algorithms and FLOP counts for the Galerkin and LSPG ROMs. The following algorithms and FLOP counts are reported:

An explicit Euler update to the Galerkin ROM (Algorithm 7, Table 8).

An implicit Euler update to the Galerkin ROM using Newton’s method with Gaussian elimination (Algorithm 8, Table 9).

An implicit Euler update to the Least-Squares Petrov–Galerkin ROM using the Gauss-Newton method with Gaussian elimination (Algorithm 9, Table 10).

References