Nonlinear model reduction: a comparison between POD-Galerkin and POD-DEIM methods

05/06/2020 ∙ by Denis Sipp, et al. ∙ Universidad de Cádiz Imperial College London ONERA 0

Several nonlinear model reduction techniques are compared for the three cases of the non-parallel version of the Kuramoto-Sivashinsky equation, the transient regime of flow past a cylinder at Re=100 and fully developed flow past a cylinder at the same Reynolds number. The linear terms of the governing equations are reduced by Galerkin projection onto a POD basis of the flow state, while the reduced nonlinear convection terms are obtained either by a Galerkin projection onto the same state basis, by a Galerkin projection onto a POD basis representing the nonlinearities or by applying the Discrete Empirical Interpolation Method (DEIM) to a POD basis of the nonlinearities. The quality of the reduced order models is assessed as to their stability, accuracy and robustness, and appropriate quantitative measures are introduced and compared. In particular, the properties of the reduced linear terms are compared to those of the full-scale terms, and the structure of the nonlinear quadratic terms is analyzed as to the conservation of kinetic energy. It is shown that all three reduction techniques provide excellent and similar results for the cases of the Kuramoto-Sivashinsky equation and the limit-cycle cylinder flow. For the case of the transient regime of flow past a cylinder, only the pure Galerkin techniques are successful, while the DEIM technique produces reduced-order models that diverge in finite time.



There are no comments yet.


page 6

page 9

page 24

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Reduced-order models play an important role in many fluid applications. Rapid evaluation of fluid systems in multi-query situations (such as Monte-Carlo techniques or parameter sweeps) or the low-dimensional representation of input-output behavior (such as in control design) are but two of many applications where a more compact, yet accurate description of fluid systems is a critical component in the overall analysis or design process. Model reduction tries to capture the essential features of the full-size system with far fewer degrees of freedom at a fraction of the computational cost. While the reduction of linear systems can be firmly based on a mathematical framework involving linear-algebra techniques, the reduction of nonlinear models is far less developed and understood.

Most commonly, model reduction for high-dimensional fluid systems relies on a projection of the governing equations onto a proper, low-dimensional basis. This basis is often extracted from a series of snapshots generated by the governing equations, and thus contains coherent structures that are deemed pertinent to the transport process one wishes to capture or model. The proper orthogonal decomposition (POD) Berkooz et al. (1993), producing a hierarchy of flow fields that optimally (in the

-sense) express the variance of fluctuations about a mean state, is quite popular among numericists and experimentalists. Once the basis has been determined, the governing equations are then recast, via a Galerkin projection, into a nonlinear dynamical system for the coefficients of the expansion basis: only the most energetic POD-modes are retained in the projection; modes with negligible energy content are dismissed by truncating the Galerkin expansion.

First efforts in this direction Berkooz et al. (1993); Aubry et al. (1988); Gillies (1998) revealed that the use of a small number of POD modes required a modification of the projection procedure, which accounts for the collective effect of the neglected POD modes, in order to ensure stability of the resulting dynamical system. In particular, the introduction of a shift-mode Noack et al. (2003) notably improved the stability and robustness of the model-reduction method. Prompted by this observation, various calibrations of the nonlinear Galerkin models were explored and investigated Couplet et al. (2005); Bergmann et al. (2009); Cordier et al. (2013), while the use of these models in flow control applications became increasingly common King et al. (2005); Bergmann et al. (2005); Bergmann and Cordier (2008).

The computational savings that are obtained in POD-Galerkin models greatly depend on the underlying structure of the governing equations of the full-order model. For instance, in the case of linear systems, the projection of the dynamics onto the subspace spanned by the selected POD modes leads to a reduced system that can be evaluated at a negligible computational cost. In the case of general nonlinear dynamics, however, the evaluation cost of the projected nonlinear model is still comparable to that of the full-order model.

To circumvent this limitation, an alternative technique known as POD-DEIM Chaturantabut and Sorensen (2010) advocates a different treatment of linear and nonlinear terms. While linear terms are treated in exactly the same fashion as in POD-Galerkin models, an additional POD basis is introduced to represent the nonlinear terms in the reduced-order model. Then, nonlinear terms are incorporated into the reduced system according to their values at selected interpolation points. The location of these interpolation

points is determined according to a greedy algorithm, known as the Discrete Empirical Interpolation Method (DEIM), that minimizes the nonlinear residual. The effectiveness of this technique relies then on the ability to cheaply evaluate nonlinear terms at the interpolation points. The reader interested in error estimates is referred to 

Chaturantabut and Sorensen (2012); Wirtz et al. (2014) for details. Over the last decade, numerous extensions of the DEIM have been proposed to either tailor this method to specific applications (see, for instance, the matrix Wirtz et al. (2014); Negri et al. (2015); Bonomi et al. (2017) or the unassembled Tiso and Rixen (2013) variants) or to improve the quality of the reduced-order representation (see, for instance, the adaptive Peherstorfer and Willcox (2015); Feng et al. (2017), localized Peherstorfer et al. (2014), trajectory-based Tan et al. (2019), non-negative Amsallem and Nordström (2016), QR-factorization based (Drmač and Gugercin, 2016) and weighted Drmač and Saibaba (2018) variants). In the context of fluid flows, this technique has been applied successfully to incompressible Xiao et al. (2014) and compressible Fosas de Pando et al. (2016) cases, as well as reacting flows Huang et al. (2018).

Despite notable success over the past years, the field of nonlinear model reduction for fluid systems is marked by empiricisms and heuristics for the choice of basis, the treatment of nonlinearities, the enforcement of physical constraints, or the selected model order. It is the objective of this paper to compare the performance of POD-DEIM models with the more traditional Galerkin methods. More specifically, we will study the convergence of the models as the number of POD modes increases; we are, however, not interested in the calibration of models of very small size or in data-driven regression techniques

(Peherstorfer and Willcox, 2016), based on sparsity-promoting techniques (Loiseau and Brunton, 2018) or on linear models (Brunton et al., 2017; Alomar et al., 2020). We also propose a new reduction technique for the nonlinear terms that takes advantage of a supplementary POD basis for the representation of the nonlinear terms: however, instead of using interpolation to determine the coefficients (like in the DEIM technique), we proceed straightforwardly by projecting the nonlinear terms onto this additional basis.

The outline of the article is as follows. Section 2 is devoted to the presentation of two test-cases: (i) cylinder flow at , and (ii) a model problem, the non-parallel Kuramoto-Sivashinsky equations. For each case, we will introduce different trajectories, a transient initialized by the fixed-point solution, a limit-cycle solution and a transient initialized by a mean-flow solution. These trajectories may be considered both for the building of the projection bases and for evalution of the models. Section 3 introduces the three model reduction techniques: (i) Galerkin projection with a single POD basis, (ii) Galerkin projection with two POD bases, and (iii) the POD-DEIM technique. In the same section, we will analyze the various properties of the reduced-order models, in particular in view of the energy-preservation of the nonlinear convective terms. Sections 4, 5 and 6 then apply the introduced techniques to the different test-cases, with various choices for the trajectories used for building the projection bases. Section 7 offers a summary of our main results and concluding remarks.

2 Test cases, trajectories, numerical discretization

The objective of this article is to generate low-order models that accurately, stably and robustly reproduce the full nonlinear dynamics of the underlying governing equations on given trajectories. We choose two systems that undergo a Hopf bifurcation, settling into a finite-amplitude limit-cycle behavior. First, as a fluid system of this type, we consider the incompressible flow past a cylinder at a supercritical Reynolds number of Second, we consider a simpler, one-dimensional model problem: the non-parallel version of the Kuramoto-Sivashinsky equation. This model equation mimics the behavior of cylinder flow, but allows a simpler and more straightforward analysis of various model-reduction techniques. We describe in the next two subsections (§ 2.1, and § 2.2) the respective governing equations, various considered trajectories and the adopted (spatial and temporal) numerical discretization.

2.1 Flow past a cylinder

We consider flow past a cylinder at a Reynolds number of which is governed by the incompressible 2D Navier-Stokes equations, made non-dimensional with the upstream velocity and the cylinder diameter. Below, we first (§2.1.1) describe the different equivalent formulations in perturbative form around the fixed point or the time-averaged flow around the limit-cycle and the spatial and numerical discretization details.

2.1.1 Governing equations

Considering finite elements, the semi-discretized form of the Navier-Stokes equations governing the composite velocity-pressure variable can be written as (Sipp and Lebedev, 2007):




and matrix designating the mass matrix linked to the finite element discretization of one velocity component. Throughout our study, we use the finite-element package FreeFEM++ (Hecht, 2012) to implement all computations. The unknown is discretized using finite elements (Arnold et al., 1984) on a triangular mesh. The mesh contains 33586 triangles and extends from to in the streamwise direction and to in the cross-stream direction, with the cylinder located at . Uniform Dirichlet boundary conditions are imposed at the inlet boundary, no-slip conditions at the cylinder surface, symmetric boundary conditions at the lateral boundaries and no-stress outflow conditions at the outlet boundary. A unknown holds 71499 degrees of freedom, the region in the vicinity of the cylinder exhibiting an isotropic mesh with triangles of size . This corresponds to a rather coarse mesh: yet, all relevant features of cylinder flow are sufficiently captured, the flow becoming unstable for

with a marginal eigenvalue appearing at

. These are classical values for such a configuration (Barkley, 2006).

In supercritical cylinder flow, there are two specific flow fields of interest, around which we may consider perturbations: w(t) &=& w_b+w’(t)
&=& ¯w+w”(t), where is a fixed point of the Navier-Stokes equations (base-flow) and the time-averaged flow (mean-flow): r(w_b)&=&0
¯w&=&lim_T →∞ 1T ∫_0^T w(t)  dt. The base-flow solution is computed in a classical manner (Sipp and Lebedev, 2007) using Newton’s method, based on a direct Lower-Upper (LU)-solver (in our case the MUMPS package Amestoy et al. (2001)).

These two changes of variables will provide two alternative governing equations in either the or variables, with a different split between linear and nonlinear operators. Depending on the case, and as shown below, the linear operator is either unstable when selecting the formulation involving , or nearly marginal if the formulation based on is chosen. As the modelling strategies may be different for the linear and the nonlinear operators, deciding on either of the two formulations will not be equivalent.

2.1.2 Perturbative form around the base-flow (BF formulation)

When considering a trajectory of the flow that starts in the vicinity of the base-flow, it may seem important to accurately capture the linear dynamics close to the base-flow and hence explicitly introduce the linearized operator around the base-flow in the equations that will be reduced. In doing so, we separate the dynamics into two distinct parts: a linearized part, which exhibits the well-known unstable global mode, and a nonlinear part which imposes a stabilizing effect. In so far, we aim at structurally reproducing the dynamics of a Stuart-Landau amplitude equationSipp and Lebedev (2007):


In such an equation, the first term on the right-hand side represents the exponential instability and induces perturbation growth; the second term exerts a restoring stabilizing force, and the equilibrium between these two effects results in a saturated limit-cycle of a given amplitude.

The semi-discretized form of the perturbed equations can then be written as (BF formulation in the following):


where matrix is the Jacobian of around ,


and the term refers to the quadratic convection term. The first term on the right-hand side corresponds to the (destabilizing) linear dynamics around the base-flow, the second one to the (stabilizing) nonlinear term that maintains the flow on the limit-cycle.

In the context of finite elements, weak forms should be favored wherever possible for the spatial discretization, as done in the definition of the residual (eq. (1)). We adopt a classical weak formulation for the linear dynamics . Yet, for computational efficiency when time-marching the equations, the convection term can also be discretized as


where and are the weak forms of the derivative matrices and denotes the element-wise (Hadamard) product. In view of the DEIM reduction, we note that this choice of discretization of the quadradic term is not pointwise, i.e., the evaluation of the nonlinearity at a given location does not correspond to the evaluation of a nonlinear function at that location

. Yet, it still provides a quick and efficient evaluation of the nonlinear term since only four sparse matrix-vector products and two element-wise products are required. For computational efficiency of the DNS solver, the various sparse matrices

, , and should be evaluated and stored in a preliminary stage of the DNS computations.

As mentioned in the introduction, a specific projection basis for the nonlinear term will be considered here (in addition to the projection basis for the state ). To this end, instead of considering the weak nonlinear term for the snapshot series, we prefer to choose a physically more meaningfull nonlinear field, , such that:


For a quick evaluation of the nonlinear term with a DEIM reduced order modelling technique (Chaturantabut and Sorensen, 2010), we finally choose the following (symmetrized, see advantages below for the implied structure of the nonlinear reduced order model) discretization scheme for :


The evaluation of this nonlinear term involves additional inverses of the mass matrix , which may be efficiently handled by a conjugate-gradient solver using a diagonal preconditioner. As mentioned earlier, this choice of discretization of the convection term stands as an approximation of the true weak-form discretization used in . We have verified that we obtain the same trajectories with both discretizations (the time-evolution of various signals in the flow field are indistinguishable), validating our choice. In what follows, we have considered, for all simulations, this implementation of the nonlinear quadratic terms.

For the time integration, we use a second-order semi-implicit scheme, with the linear operator in being inverted by the direct LU solver at each time-step (Cerqueira and Sipp, 2014). The time-step is set to , which ensures a maximum CFL number (based on the base-flow velocity) of 0.13.

2.1.3 Perturbative form around the mean-flow (MF formulation)

The evolution equation for the perturbation is given as (MF formulation):


where b &=& r(¯w)
A”&=&∂r∂w—_¯w, are, respectively, the residual of the discretized Navier-Stokes equations taken at the mean-flow and the Jacobian of the residual around this flow. We will use the same discretization choices for as those presented in the previous section.

These are the governing equations that naturally arise when the trajectory is taken along the limit-cycle (Bergmann et al., 2009; Fosas de Pando et al., 2016). It is also the natural choice if only the mean-flow is known (and not the base flow We remark that this equation is now inhomogeneous with a constant term

and that the linear operator now corresponds to the Navier-Stokes equations linearized around the mean-flow, which exhibits stability properties different from the ones involving the base-flow. It is expected that the flow snapshots on the limit-cycle are close to the features of the mean-flow marginal eigenvector. For a harmonic flow (which is a reasonable approximation for cylinder flow), the frequency of the flow field and its Fourier mode correspond to a marginal-eigenvalue and eigenvector of the linearized operator around the mean-flow

(Mezić, 2013; Turton et al., 2015). Therefore, eq. (9

) likely represents the best choice to reproduce the dynamics around the mean-flow, since most of the features of the limit-cycle are already captured by the linear operator, while the nonlinear operator only needs to correct small defects of the linear representation (for example, the eigenvalue might be slightly unstable). Also, this equation is at the heart of all successfull mean-flow-based resolvent studies, which justify the validity of the approach by different arguments (dominant singular value condition

(Beneddine et al., 2016)

, white noise approximation of the forcing

(Towne et al., 2018) or small-amplitude assumption of the fluctuation field (Leclercq et al., 2019)). Of course, in the present case, the dynamics close to the base-flow, that is the frequency shift of the dynamics and the different wavelengths of the perturbations, needs to be entirely modelled by the nonlinear terms. Hence, it is likely that the linear dynamics close to the base-flow will be more poorly represented, since nonlinear model reduction generally performs poorer than linear reduction. As in the previous section, we will consider snapshots for building the nonlinear projection basis.

2.1.4 Transient (TR), Limit-Cycle (LC) and Mean-Flow Transient (MFTR) trajectories

The base- and mean-flows at are represented in figures 1(a,b) with iso-values of streamwise velocity. The mean-flow solution shows a markedly reduced recirculation zone behind the cylinder, when compared to the base-flow solution.


x y u


x y u

Figure 1: (a): Base- and (b): mean-flow solutions for flow past a cylinder at a supercritical Reynolds number of visualized by contours of the streamwise velocity component.

Normal modes of the linear operators can be formulated as: and . These are solution of the eigen-problems A’^w^BF&=&λ^BF Q^w^BF
A”^w^MF&=&λ^MF Q^w^MF, which can be treated by Krylov methods linked to a shift-invert strategy based on the direct LU-solver (Sipp and Lebedev, 2007). As expected, we obtain a pair of unstable eigenvalues for and a pair of (nearly-) marginal eigenvalues for . These compare well with those in the literature (Barkley, 2006). Note also that the value of the frequency is very close to the frequency of the limit-cycle, as obtained from the DNS solution ( here). A close-up view of both spectra in the vicinity of the two previously mentioned eigenvalues is shown in figure 2(a), with black symbols for and red symbols for . Also, the horizontal black solid line depicts the frequency The real part of the streamwise velocity of the corresponding eigenvectors are shown with lines in corresponding colors in figure 2(b). The shapes of the two eigenvectors are very different: the base-flow eigenvector exhibits a gradual growth of the oscillation amplitude in the downstream direction, while the mean-flow eigenvector exhibits a peak around , before slowly decreasing. The last pattern is representative of oscillations of the flow on the limit-cycle. These observations are reminiscent of the property that an eigenvalue/eigenvector of a mean-flow solution exactly reproduces the frequency and spatial structure of the unsteady solution, if the unsteady solution exhibits a harmonic behavior (Barkley, 2006; Sipp and Lebedev, 2007; Mezić, 2013; Turton et al., 2015). This property is well satisfied in the case of cylindar flow (Turton et al., 2015).

(a) (b)
s w aBF bMF c x v aBF bMF
Figure 2: (a): Eigenspectrum of the linearized operator around the base-flow (black symbols) and mean-flow (red symbols). The horizontal solid line depicts the frequency of the limit-cycle obtained by DNS. (b): Real part of the streamwise component of the associated eigenvectors.

A direct numerical simulation has been carried out, starting with the small-amplitude unstable global mode as an initial condition, , with standing for a small amplitude. The evolution of the perturbation towards the limit-cycle is best visualized by the time-evolution of the kinetic energy given by the quadratic form Figure 3(a) clearly displays (see the black solid line labelled TR) an exponential instability over about time-units (at an amplification rate corresponding to twice the amplification rate of the unstable global mode), until saturation sets in, as the elimit-cycle behavior is reached (red solid line labelled LC). We also have shown with a solid magenta line the transient MFTR trajectory, which is initialized by the mean-flow solution about the limit-cycle: after a quick decrease of the kinetic energy, the perturbation grows again showing a near exponential growth rate, before saturation on the limit-cycle sets in. We will use this particular trajectory to assess the robustness of the reduced-order models, by evaluating their performance on a trajectory that was not considered for the building of the model. A more detailed view of the same evolutions is given by the time-trace of the streamwise velocity component at point in the wake of the cylinder (see figure 3(b)). Again, an exponential instability is clearly discernible for the black solid line, before convergence towards a limit-cycle behavior sets in (red solid line). The vortex shedding frequency is visible in the velocity trace.

(a) (b)
Figure 3: (a) For flow past a cylinder at , time evolution of the kinetic energy of the perturbation about the base-flow for the transient TR from the base-flow (black line), for the limit-cycle solution (red line) and the transient MFTR from the mean-flow (magenta line). (b) Time trace of the streamwise velocity extracted at for each trajectory TR, LC, MFTR. The TR and MFTR trajectories span the time-range and the LC trajectory the range .

In addition, we show representative snapshots of the perturbation field and on the limit cycle (see figures 4(a,b)), visualized by iso-contours of the streamwise velocity component. In the left snapshot, the dominant red colors in the central part of the wake represent the mean-flow deformation with respect to the base-flow (the recirculation bubble shortens, inducing positive mean-values of ), and the antisymmetric streamwise large-scale modulations depict the vortex-shedding mode. In the right snapshot, only the vortex shedding mode is visible, since the mean-deformation was subtracted by the change of variables. These differences will yield different projection bases depending on the considered formulation, since the snapshots are used in the BF formulation and the snapshots in the MF formulation.


x y u


x y u

Figure 4: (a): Representative snapshot of the perturbation field at from the limit cycle of flow past a cylinder at Reynolds number visualized by iso-contours of the streamwise velocity component. (b): Same snapshot but represented with the variable.

2.2 Non-parallel Kuramoto-Sivashinsky equation

This equation is commonly used as a model equation for complex fluid motion as it contains many features in a one-dimensional setting that have equivalents in higher dimensions. It is thus a valuable proxy for investigating analytical and computational techniques and for quantifying the influence of its ingredients on user-specified performance measures. We will use this equation for assessing the accuracy, stability and robustness of model-reduction methods in a simpler case. To this end, we consider the non-parallel version of the Kuramoto-Sivashinsky equation in the form




and as a real function defined on the interval To complete the problem, we provide the boundary conditions on according to &u’=∂_xx u’ = 0   &at  x=-x_∞,
&∂_x u’ =∂_xxx u’ = 0   &at  x=x_∞. The term models uniform convection at a prescribed speed of and constitutes the nonlinear, quadradic convection term which is also present in the full Navier-Stokes equations. The expression in the Kuramoto-Sivashinky equation models an instability with a strength of its origin located about and its spatial extent governed by the width parameter This term mimics the streamwise-localized instability mechanism acting in the recirculation bubble of cylinder flow. The final term provides a stabilizing hyper-diffusion which damps high-wavenumber (small-scale) structures.

In view of spatial discretization with classical finite-element methods, we consider the auxiliary variable to render the system second-order, so that:


with at and at . The second equation can be seen as a constraint on the state , reminiscent of the divergence-free constraint in the Navier-Stokes equations.

Semi-discretization of the above equation with finite elements yields an expression whose structure is similar to (4) for the state-variable . We obtain




as the weak-form of


and the symmetrized bilinear term:


where, again, the point-wise discretization scheme of the nonlinear term has been chosen.

We use second-order elements for the discretization of . For the mesh, we choose and the elements are of size , which yields 4000 elements. The time-step for the simulation is . The same semi-implicit time integration strategy as described for the Navier-Stokes equations is used here.

For our analysis, we choose the following constants


which leads to an unstable linearized operator (about the base flow ) with a single pair of unstable eigenvalues at The eigenvalue spectra for the linearized operators associated to the base-flow (BF) and mean-flow (MF) formulations are shown with red and black symbols in figure  5(a), and the (real part of the) global mode associated with the leading eigenvalue in figure 5(b) with the same color. Figures 5(c,d,e) depict the TR, MF and MFTR trajectories in a similar manner than in subsection § 2.1.4, i.e., the evolution of the perturbation kinetic energy (figure 5(c)) for the three trajectories TR, LC and MFTR, the time-trace of the variable at (figure 5(d)) and four representative snapshots from the limit cycle (figure 5(e)). We have also shown the mean-flow solution corresponding to the limit-cycle in figure 5(f). Coming back to the mean-flow eigenvalue , we observe that, similarly to the case of cylinder flow, the amplification rate of the linearized operator around the mean-flow (MF) is closer to marginality. Yet, the frequency of the mode has decreased well below the frequency of the limit-cycle, which is . The eigenvector of the mean-flow unstable eigenmode is now closer to the oscillations of the limit-cycle. The observed discrepancies are due to the fact that the limit-cycle for the KS equation exhibits many harmonics and is therefore far from harmonic. The ”harmonic property” is better satisfied in the case of cylinder flow. To conclude, despite these differences, we can state that the behavior of this model problem is qualitatively rather close to the case of flow past a cylinder: the solution undergoes a Hopf bifurcation and progresses towards a limit-cycle behavior.

(a) (b)
s w aBF bMF c x v aBF bMF
(c) (d)
(e) (f)
x w a b c d x w
Figure 5: Kuramoto-Sivashinsky equation. (a): Eigenvalue spectrum of linear operator in BF formulation ( variable) and MF formulation ( variable). The horizontal solid line depicts the frequency of the limit-cycle obtained by DNS. (b): Corresponding real parts of leading eigenvectors. (c): Temporal evolution of kinetic perturbation energy The initial condition corresponds to a small-amplitude global mode (see subfigure (b), in black). (d): Temporal evolution of evaluated at TR and LC refer, respectively, to the time intervals and . Trajectory MFTR (magenta line) has been initialized by the mean-flow solution about the limit-cycle. (e): Four snapshots within one period of the limit-cycle regime. (f): Mean-flow solution about the limit-cycle.

3 Model reduction methodology

In this section, we develop a mathematical framework for the projection-based model reduction of semi-discretized evolution problems of the form


where is a bilinear symmetric term. We express the perturbation as where the matrix contains basis vectors as its columns and denotes the vector of expansion coefficients, and project the resulting equation onto the subspace spanned by . The reduced-order model then reads


with W^H Q W &=& I,
N_ijk &=& N_ikj. The last relation stems from the symmetry of the nonlinear operator. The subscripts indicate the

-th component of the corresponding vectors, matrices and tensors. The values of the coefficients

, and will be given in the next sections. The operation count associated with this reduced-order model scales as , since for each of the degrees of freedom, we have a double-summation over terms. Yet, as shown below, in the case of the DEIM method, the scaling of the operation count linked to the evaluation of the nonlinear term can be drastically reduced in the case of a point-wise nonlinearity (as is the case here).

Based on the above approach, the building of the reduced-order model involves three steps: the computation of the POD basis for the representation of the state (section §3.1), the modelling of the constant term (section §3.2), of the linear term (section §3.3) and of the nonlinear term (section §3.4). Section §3.5 discusses the mathematical properties of the model while §3.6 presents all quality and error measures that will serve to assess the models.

3.1 The orthonormal basis

The orthogonal basis contained in the columns of the matrix are computed from snapshots gathered from numerical simulations of the full system. With denoting the snapshot matrix, consisting of columns that constitute the flow fields at equispaced instants in time, we form the temporal correlation matrix and determine its eigenvalues and eigenvectors according to


The matrix contains the non-negative eigenvalues along its diagonal, ranked in decreasing order, and the columns of provide the corresponding eigenvectors. The POD-modes, and thus the basis matrix is then formed following


By construction, the columns of are orthonormal; the corresponding matrix satisfies the orthonormality condition The POD-modes associated to the largest eigenvalues form an optimal basis to represent the snapshots contained in .

3.2 The constant term

The constant term in the reduced-order model can easily be determined from the constant term of the full-scale equation (19). It is simply obtained by a left-multiplication of the full-scale constant term by which is equivalent to a Galerkin projection of onto the basis spanned by the columns of i.e. the POD-modes. Mathematically, we have


3.3 The linear term

The projection of the linear term onto the reduced-order basis yields a reduced system matrix that may be computed following:


The number of POD modes in (i.e., the number of columns of ) may be chosen such that the spectrum of the reduced matrix displays unstable eigenvalues close to the unstable eigenvalues of the large-scale generalized eigenproblem defined by .

3.4 The nonlinear term

While there is little choice for the reduced expression of the constant and linear terms, the treatment of the nonlinear terms allows more choice and flexibility. In this article, three different techniques are investigated: the traditional Galerkin projection of the nonlinear terms onto the reduced basis (method 1, see section § 3.4.1), the Galerkin projection of the nonlinear terms onto a new dedicated basis representing only the nonlinearities (method 2, see section § 3.4.2) and the Discrete Empirical Interpolation Method (DEIM) applied to the dedicated basis (method 3, see section § 3.4.3).

3.4.1 Method 1: traditional Galerkin projection

This method follows the traditional derivation of a nonlinear dynamical system for the coefficients In it, the arguments of the bilinear function are expressed in terms of their Galerkin expansion, after which the nonlinear expression is left-multiplied by the matrix which is equivalent to a projection onto our orthonormal basis. We obtain the nonlinear coefficients as


where we have used the common notation indicating the -th column of Since is symmetric in its arguments, we have for all indices

3.4.2 Method 2: Galerkin projection with an additional nonlinear basis

Besides the common basis extracted from the snapshot sequence the second method introduces a second basis intended to represent the nonlinear terms and thus nonlinear effects. For this, we evaluate the nonlinearities for all snapshots in to form a second snapshot sequence which we refer to as Analogous to the first method, the correlation matrix based on this second sequence is then decomposed into its eigenvalues and eigenvectors according to


and the corresponding POD modes are determined as


As before, the structures contained in the columns of are orthonormal by construction, which is expressed mathematically as

The nonlinear terms may then be projected onto the basis


with as the coefficient vector. Using this new basis to express the nonlinearities in the governing equations, we arrive at the nonlinear terms of our reduced-order model as


Again, the symmetric nature of leads to

3.4.3 Method 3: Discrete Empirical Interpolation Method with an additional nonlinear basis

The basis introduced in section § 3.4.2 for the representation of the nonlinear term is considered again, but the coefficients are obtained differently. If contains columns representing structures onto which we project the nonlinear term, we enforce equation not in a least-squares sense as above, but instead by enforcing equality at selected interpolation points. Mathematically, we premultiply the above equation by a row-selector matrix which yields


The matrix contains columns, each displaying a single unit value at some row with the remaining entries as zero. The choice of these columns and interpolation points follows a greedy algorithm and is given below. The premultiplication by ensures that is invertible and, as a consequence, the above equation can be solved for the coefficient vector We thus have and invoking the bilinearity of the nonlinear operator:


we obtain the following representation of the nonlinear term


Note that this expression does not require the nonlinear term to be point-wise. It is the bilinearity of the nonlinear operator that ensures an operation count of the order .

In the case of a point-wise nonlinearity, the evaluation of the nonlinear term in the reduced order model can be achieved at a very small cost, that is the evaluation of the nonlinearity at points , where is a nonlinear operator taking a vector of size and giving back a vector of size . In such a case, the full reduced-order model is:


where and are matrices of size and respectively. Hence, the operation count scales as . This implementation is consistent with equation (32). If the nonlinearity is not pointwise, the sparsity argument of (Chaturantabut and Sorensen, 2010) or the introduction of auxiliary variables (Fosas de Pando et al., 2016) may also lead to a reduction of the evaluation cost of the nonlinear term.

Data: nonlinear basis
Result: row-selector matrix
= zeros();
= argmax();
= ;
for  =  do
       = ;
       = ;
       = argmax();
       = ;
end for
return ;
Algorithm 1 DEIM algorithm (adapted from Chaturantabut and Sorensen (2010)).

3.5 Mathematical properties of the reduced-order models

When comparing the stability, accuracy and robustness of various model reduction techniques, it is imperative to introduce quality measures and other mathematical properties to quantitatively assess their absolute and relative performance. The kinetic perturbation energy will serve as the quantity that will be monitored and compared for the full and the reduced-order model.

The kinetic energy of the reduced-order model is governed by


If we define the symmetric and antisymmetric parts of and as follows (with superscript denoting the symmetric and superscript denoting the anti-symmetric part) L_ij^S = Lij+ Lji2,  LijA = Lij- Lji2,
NijkS = Nijk+ Nikj+ Njik+ Njki+ Nkij+ Nkji6, NijkA = 5 Nijk- Nikj- Njik- Njki- Nkij- Nkji6 we can recast the evolution equation for the energy as


The energy evolution equation leads to the following tight bounds (i.e., there exist that achieve the bound):


where stands for the Frobenius norm of a matrix or a tensor, defined as


The first term on the right-hand-side of equation (36) may lead to an algebraic increase of of the form


the second term to an exponential growth following


and the third term to a finite-time blow-up according to


The nonlinearities that we consider in this article are of purely convective type and are therefore energy-preserving (under standard boundary conditions). It is straightforward to show that


which implies for the reduced-order model that


If this condition holds exactly, the finite-time singularity implied by the third term in equation (36) can be avoided; if, on the other hand, the reduced-order model does not satisfy one should expect a finite-time blow-up in energy for some initial condition. The condition thus constitutes an important and effective test for the robustness and long-term stability of the reduced-order model.

3.6 Quality measures of reduced-order models

In this subsection, refers to either if the BF formulation is chosen or if the MF formulation is selected. For a quantitative assessment of the quality of model reduction, we introduce various error measures.

We first look for a fixed point of the reduced-order model using a Newton method to solve:


We evaluate the relative error between the predicted base-flow and the actual base-flow ( with the BF formulation and with the MF formulation):


We compute the leading eigenvalues/eigenvectors of the operator obtained by linearizing the equations governing the reduced-order model (20) around . The number of unstable eigenvalues is denoted and should be equal to for any formulation. We evaluate both the relative error between the leading eigenvalue and and the relative alignment error between the predicted leading eigenvector and the actual one following:


We then compute the mean-flow by time-averaging the results of a temporal simulation of the reduced-order model (20):


In a similar way then before, we define the relative errors pertaining to the mean-flow and the leading eigenvalue / eigenvector of equations (20) linearized around the mean-flow :


The number of unstable eigenvalues should again be equal to two. The recovery of the least-damped eigenvalues and is an important quality measure for the model reduction procedure since, in the transient regime (TR), the solution behaves like , while in the limit-cycle case (LC), displays near marginal stability properties with a frequency close to the frequency of the limit-cycle Barkley (2006).

The remaining two errors, proposed for the assessment of model reduction quality, are errors linked to the reconstruction of the TR, LC and MFTR trajectories. In each case, the initial conditions for the reduced-order models are set to those of the large-scale simulation,


where denotes the projection of the snapshots of the unsteady simulation onto the POD modes representing the state:


The mean truncation error , which is the ratio between the integral energy (over time) of the retained modes and the integral total energy, is given by


and the relative mean model error by:


and respectively, assess the error due to truncation of the projection basis to POD modes and the error due to prediction of the time-evolution of the retained states of the model. In Alomar et al. (2020), it was shown that the total error can be deduced from and following .

Finally, we introduce the ratio between the Frobenius norm of the symmetric part of the nonlinear term and the Frobenius norm of


With the nonlinearities representing the convective term of the Navier-Stokes equations, this ratio should be small (or zero). In contrast, a large value of this quantity may indicate robustness problems, since the presence of a non-zero -term in the model may induce a finite-time singularity for some initial conditions. We have therefore also evaluated the performance of the models when considering in the reduced-order-model instead of , that is when removing the symmetric part of the quadratic term. In this case, we denote the model error .

4 Kuramoto-Sivashinsky equation

In this section, we assess the quality of the reduced-order models for the case of the Kuramoto-Sivashinsky equation presented in section § 2.2. After describing the POD bases for the state and the nonlinearities (section § 4.1), we analyze the different reduction techniques presented in section § 3 (section § 4.2).

4.1 POD bases

In the following subsections, we present the state and nonlinear POD modes associated with two different time-spans of the learning trajectory.

4.1.1 Bases and determined from snapshots on the transient and the limit-cycle ()

(a) (b)
n s n s
(c) (d)
a b c x w a b c x w
Figure 6: Model reduction of Kuramoto-Sivashinsky equation with projection bases obtained from snapshots within for the BF formulation. (a,b): Eigenvalues of the correlation matrices for the snapshots representing (a): the state and (b): the nonlinearity. (c,d): first three POD modes of (c): state and (d): nonlinearity. Colored circles show the corresponding DEIM points (only the position is relevant).

We take snapshots from the simulation shown in figures 5(c–e) over the time-span to build representative bases and of the transient and the limit-cycle. The sampling time is time units. Based on these snapshots, the leading eigenvalues of the correlation matrices for the state and the nonlinearities are shown, for the BF formulation, in figures 6(a,b) . We observe that the eigenvalues of the correlation matrix for the nonlinearity snapshots drop off markedly faster, indicating that the effective dimensionality of the space spanned by nonlinearity snapshots is lower than the analogous dimensionality formed by the state snapshots.

In figures 6(c,d), we present the first three POD modes of the state and of the nonlinearity . We observe that, while the state POD modes are non-zero over the entire region the support of the nonlinearity POD modes is far more compact (within ). This result is consistent with the pronounced drop-off in the eigenvalues of the respective correlation matrices shown in figure 6(a,b). It is also seen that the first state POD mode accounts for the mean-flow deformation shown in fig. 5(f). The other state POD modes allow us to represent, in an optimal way, the unstable growing mode shown in fig. 5(b) with a black solid line as well as the time-oscillations depicted in 5(e). In figure 6(d), we have displayed also the associated DEIM points with symbols. As is evident from the plot, the DEIM points approximately correspond to the location of maximum amplitude of the underlying POD modes.

4.1.2 Bases and determined from snapshots on the limit-cycle ()

(a) (b)
n s n s
(c) (d)
a b c x w a b c x w
Figure 7: Model reduction of Kuramoto-Sivashinsky equation with projection bases obtained from snapshots within for the MF formulation. Same caption as in fig. 6.

We take again snapshots from the simulation shown in figures 5(c–e) but over the time-span (not represented in the figure), which is only representative of the limit-cycling behavior. The sampling time remains unchanged, and the results for the MF formulation (which is consistent with the sole knowledge of the limit-cycle solution) are shown in fig. 7 in a similar way than in fig. 6. The dimensionality of the dynamics is drastically reduced in both cases, since only ten POD modes are required to achieve a decrease of four orders of magnitude of the eigenvalues (compared to approximately 30 POD modes in the previous section). The eigenvalues for the state variable approximately come in pairs and therefore represent spatio-temporal structures that are smoothly convected downstream. The state-POD mode representing the mean-flow deformation has disappeared since this feature is not visible in the variable over the time-span of the learning trajectory. Only POD modes representing limit-cycle oscillations are seen in fig. 7(c).

4.2 Model reduction

Based on the snapshot bases for the state and the nonlinearities, we are now in a position to explore the various reduced-order models and compare them with the corresponding large-scale unsteady simulations. We will consider for that the TR, LC and MFTR trajectories presented in fig. 5(c,d). For the first reduction method (Galerkin projection with a single basis), we choose a number of POD modes to represent the state ( or following the chosen formulation). For the second method (Galerkin projection with two bases) and third method (Galerkin projection for the linear term and DEIM for the nonlinear term), we choose the number of POD modes for either basis, labelled and .

1B-60 12 2 0.1 0.00 2 2 1 0.1 0.01 0.5 0.5 0.00 0.1 0.3 1
1B-50 13 2 0.5 0.00 2 2 1 0.1 0.03 3 3 0.02 0.2 0.3 3
1B-40 14 2 2 0.03 2 2 1 0.2 0.2 13 13 0.1 3 1 38
1B-30 12 2 2 10 0.3 2 2 2 1 85 85 1 3 4 49
1B-20 4 8 45 15 4 2 10 3 4 103 103 4 9 9 27
1B-10 0.3 4 45 15 29 2 28 10 19 100 100 7 75 75 81
1M-60 12 2 0.1 0.00 2 2 1 0.1 0.01 2 8 0.01 0.1 0.1 1
1M-40 14 2 3 0.04 2 2 1 0.2 0.3 119 117 0.2 4 4 69
1M-20 5 6 74 17 1 2 9 3 5 140 140 7 23 23 49
2B-60-60 14 2 0.1 0.00 2 2 1 0.1 0.01 1 12 0.00 0.1 13 1
2B-60-40 20 2 0.1 0.00 2 2 1 0.2 0.01 3 15 0.00 1 35 4
2B-60-20 24 2 0.1 0.00 2 2 1 1 0.01 12 18 0.00 4 26 10
2B-40-40 16 2 2 0.03 2 2 1 0.2 0.2 13 14 0.1 2 13 41
2B-40-20 18 2 2 0.03 2 2 1 1 0.2 17 18 0.1 5 14 45
2B-20-20 8 8 45 15 4 2 7 3 4 104 103 4 12 11 25
2B-10-10 19 4 45 15 35 0 25 7 19 100 100 7 88 106 71
2M-60-60 14 2 0.1 0.00 2 2 1 0.1 0.01 3 124 0.01 1 4 4
3B-60-60 62 2 0.1 0.00 1 2 1 0.1 0.01 1 10 0.00 1 26 5
3B-60-40 65 2 0.1 0.00 2 2 1 0.2 0.01 3 7 0.00 1 28 6
3B-60-20 67 2 0.1 0.00 2 2 1 1 0.01 16 61 0.00 3 47 12
3B-40-40 58 2 2 0.03 2 2 1 0.2 0.2 14 13 0.1 2 8 42
3B-40-20 60 2 2 0.03 2 2 1 1 0.2 21 29 0.1 4 14 43
3B-20-20 30 8 45 15 4 2 7 3 4 103 103 4 13 31 23
3B-10-10 23 4 45 15 39 0 97 75 19 100 100 7 79 96 80
3M-60-60 64 2 0.1 0.00 2 2 1 0.1 0.01 4 114 0.01 1 5 3
Table 1: Error analysis of various reduced-order modelling techniques for the Kuramoto-Sivashinsky equation with projection bases and based on snapshots taken from the transient and the limit-cycle. M, refers to methods 1, 2 and 3, F to base-flow (B) and mean-flow (M) formulations, and to the number of POD modes representing the state and the nonlinearity (when applicable), to the ratio between and . The four next columns concern the recovery of the base-flow (BF) properties: is the relative error with respect to the true base-flow solution, the number of unstable eigenvalues of the operator linearized around the base-flow, the relative error with respect to the most unstable eigenvalue , the relative alignment error with respect to the most unstable eigenvector . The next four columns give analogous results for the mean-flow properties (mean-flow solution , number of unstable eigenvalues, eigenvalue , eigenvector ). The next three columns deal with the recovery of the TR trajectory: , and are the mean truncation error, model error and model error when setting to zero the symmetric part of the reduced-order model. The next three columns provide the same information for the MF trajectory, while the last column deals with the recovery of the MFTR trajectory (transient initialized with the mean-flow solution). We have shown in green reduced-order simulations that achieve a relative error less than , in orange those achieving an error between and , and in red the remaining ones. All quantities are given in percentage.

All results from the simulations are summarized in table 1 for bases and determined from the transient and limit-cycle behavior (fig. 6) and in table 2 for bases solely obtained from the limit-cycle trajectory (fig. 7). These tables provide all parameters of the tested models and the associated errors. We display in green errors that are less than , in orange those between and and in red those who exceed or simulations that diverge. We have also indicated in green or red whether the model recovers or not the right number of unstable eigenvalues of the equations linearized around the fixed-point and time-averaged solutions. This is an important feature for subsequent use of the models in a flow control strategy, especially for the base-flow eigenvalue.

From table 1, which deals with the case of POD bases built with snapshots in the transient and on the limit-cycle, we can draw the following conclusions :

  • Comparing 1B-60, 2B-60-60 and 3B-60-60, all methods manage to achieve excellent and nearly equivalent results, if the number of POD modes, and are sufficiently high and close. In this case, we have . There thus seems to be no apparent gain in using the more elaborate two-bases strategies (methods 2 and 3).

  • Comparing 1B-60 to 1M-60, 2B-60-60 to 2M-60-60 and 3B-60-60 to 3M-60-60, it seems slightly better to use BF formulations than MF formulations, whichever method is used.

  • Comparing 1B-60 to 1B-50, the number of POD modes required for high precision for the TR trajectory, i.e., is rather large (). In particular, the recovery of the unstable eigenvalue/eigenvector ( and ) requires large values of , typically .

  • On both TR and LC trajectories, the truncation error is always small compared to the model error .

  • Considering the column , the most robust models are those provided by the first method and the second method with . The models deduced by the DEIM method are less energy-preserving with always high values of . However, a high value of does not necessarily prevent the model from being accurate (), see line 3B-60-60. With the third technique, decreasing the number of POD modes for the representation of the nonlinearities always alters the energy-preserving property. This is not the case, however, for the second method.

  • Errors are consistently smaller on the LC trajectory than on the TR trajectory.

  • Considering the recovery of the mean-flow properties (four sub-columns of MF column), it is seen that both the mean-flow solution and the linear properties around it are easily reproduced with low values of