1 Introduction
Reducedorder models play an important role in many fluid applications. Rapid evaluation of fluid systems in multiquery situations (such as MonteCarlo techniques or parameter sweeps) or the lowdimensional representation of inputoutput 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 fullsize 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 linearalgebra techniques, the reduction of nonlinear models is far less developed and understood.
Most commonly, model reduction for highdimensional fluid systems relies on a projection of the governing equations onto a proper, lowdimensional 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 PODmodes 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 shiftmode Noack et al. (2003) notably improved the stability and robustness of the modelreduction 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 PODGalerkin models greatly depend on the underlying structure of the governing equations of the fullorder 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 fullorder model.
To circumvent this limitation, an alternative technique known as PODDEIM 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 PODGalerkin models, an additional POD basis is introduced to represent the nonlinear terms in the reducedorder 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 reducedorder representation (see, for instance, the adaptive Peherstorfer and Willcox (2015); Feng et al. (2017), localized Peherstorfer et al. (2014), trajectorybased Tan et al. (2019), nonnegative Amsallem and Nordström (2016), QRfactorization 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 PODDEIM 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 datadriven regression techniques
(Peherstorfer and Willcox, 2016), based on sparsitypromoting 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 testcases: (i) cylinder flow at , and (ii) a model problem, the nonparallel KuramotoSivashinsky equations. For each case, we will introduce different trajectories, a transient initialized by the fixedpoint solution, a limitcycle solution and a transient initialized by a meanflow 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 PODDEIM technique. In the same section, we will analyze the various properties of the reducedorder models, in particular in view of the energypreservation of the nonlinear convective terms. Sections 4, 5 and 6 then apply the introduced techniques to the different testcases, 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 loworder 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 finiteamplitude limitcycle 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, onedimensional model problem: the nonparallel version of the KuramotoSivashinsky equation. This model equation mimics the behavior of cylinder flow, but allows a simpler and more straightforward analysis of various modelreduction 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 NavierStokes equations, made nondimensional 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 timeaveraged flow around the limitcycle and the spatial and numerical discretization details.
2.1.1 Governing equations
Considering finite elements, the semidiscretized form of the NavierStokes equations governing the composite velocitypressure variable can be written as (Sipp and Lebedev, 2007):
(1) 
with
(2) 
and matrix designating the mass matrix linked to the finite element discretization of one velocity component. Throughout our study, we use the finiteelement 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 crossstream direction, with the cylinder located at . Uniform Dirichlet boundary conditions are imposed at the inlet boundary, noslip conditions at the cylinder surface, symmetric boundary conditions at the lateral boundaries and nostress 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 NavierStokes equations (baseflow) and the timeaveraged flow (meanflow):
r(w_b)&=&0
¯w&=&lim_T →∞ 1T ∫_0^T w(t) dt.
The baseflow solution is computed in a classical manner (Sipp and Lebedev, 2007) using Newton’s method, based on a direct LowerUpper (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 baseflow (BF formulation)
When considering a trajectory of the flow that starts in the vicinity of the baseflow, it may seem important to accurately capture the linear dynamics close to the baseflow and hence explicitly introduce the linearized operator around the baseflow in the equations that will be reduced. In doing so, we separate the dynamics into two distinct parts: a linearized part, which exhibits the wellknown unstable global mode, and a nonlinear part which imposes a stabilizing effect. In so far, we aim at structurally reproducing the dynamics of a StuartLandau amplitude equationSipp and Lebedev (2007):
(3) 
In such an equation, the first term on the righthand 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 limitcycle of a given amplitude.
The semidiscretized form of the perturbed equations can then be written as (BF formulation in the following):
(4) 
where matrix is the Jacobian of around ,
(5) 
and the term refers to the quadratic convection term. The first term on the righthand side corresponds to the (destabilizing) linear dynamics around the baseflow, the second one to the (stabilizing) nonlinear term that maintains the flow on the limitcycle.
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 timemarching the equations, the convection term can also be discretized as
(6) 
where and are the weak forms of the derivative matrices and denotes the elementwise (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 matrixvector products and two elementwise 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:
(7) 
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 :
(8) 
The evaluation of this nonlinear term involves additional inverses of the mass matrix , which may be efficiently handled by a conjugategradient solver using a diagonal preconditioner. As mentioned earlier, this choice of discretization of the convection term stands as an approximation of the true weakform discretization used in . We have verified that we obtain the same trajectories with both discretizations (the timeevolution 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 secondorder semiimplicit scheme, with the linear operator in being inverted by the direct LU solver at each timestep (Cerqueira and Sipp, 2014). The timestep is set to , which ensures a maximum CFL number (based on the baseflow velocity) of 0.13.
2.1.3 Perturbative form around the meanflow (MF formulation)
The evolution equation for the perturbation is given as (MF formulation):
(9) 
where
b &=& r(¯w)
A”&=&∂r∂w—_¯w,
are, respectively, the residual of the discretized NavierStokes equations taken at the meanflow 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 limitcycle (Bergmann et al., 2009; Fosas de Pando et al., 2016). It is also the natural choice if only the meanflow 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 NavierStokes equations linearized around the meanflow, which exhibits stability properties different from the ones involving the baseflow. It is expected that the flow snapshots on the limitcycle are close to the features of the meanflow 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 marginaleigenvalue and eigenvector of the linearized operator around the meanflow
(Mezić, 2013; Turton et al., 2015). Therefore, eq. (9) likely represents the best choice to reproduce the dynamics around the meanflow, since most of the features of the limitcycle 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 meanflowbased 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 smallamplitude assumption of the fluctuation field (Leclercq et al., 2019)). Of course, in the present case, the dynamics close to the baseflow, 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 baseflow 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), LimitCycle (LC) and MeanFlow Transient (MFTR) trajectories
The base and meanflows at are represented in figures 1(a,b) with isovalues of streamwise velocity. The meanflow solution shows a markedly reduced recirculation zone behind the cylinder, when compared to the baseflow solution.
Normal modes of the linear operators can be formulated as:
and .
These are solution of the eigenproblems
A’^w^BF&=&λ^BF Q^w^BF
A”^w^MF&=&λ^MF Q^w^MF,
which can be treated by Krylov methods linked to a shiftinvert strategy based on the direct LUsolver (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 limitcycle, as obtained from the DNS solution ( here).
A closeup 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 baseflow eigenvector exhibits a gradual growth of the oscillation amplitude in the downstream direction, while the meanflow eigenvector exhibits a peak around , before slowly decreasing. The last pattern is representative of oscillations of the flow on the limitcycle. These observations are reminiscent of the property that an eigenvalue/eigenvector of a meanflow 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 
A direct numerical simulation has been carried out, starting with the smallamplitude unstable global mode as an initial condition, , with standing for a small amplitude. The evolution of the perturbation towards the limitcycle is best visualized by the timeevolution 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 timeunits (at an amplification rate corresponding to twice the amplification rate of the unstable global mode), until saturation sets in, as the elimitcycle 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 meanflow solution about the limitcycle: after a quick decrease of the kinetic energy, the perturbation grows again showing a near exponential growth rate, before saturation on the limitcycle sets in. We will use this particular trajectory to assess the robustness of the reducedorder 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 timetrace 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 limitcycle behavior sets in (red solid line). The vortex shedding frequency is visible in the velocity trace.
(a)  (b) 
T K aTR bLC cMFTR  ATR BLC aTR bLC cMFTR T U u v 
In addition, we show representative snapshots of the perturbation field and on the limit cycle (see figures 4(a,b)), visualized by isocontours of the streamwise velocity component. In the left snapshot, the dominant red colors in the central part of the wake represent the meanflow deformation with respect to the baseflow (the recirculation bubble shortens, inducing positive meanvalues of ), and the antisymmetric streamwise largescale modulations depict the vortexshedding mode. In the right snapshot, only the vortex shedding mode is visible, since the meandeformation 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.
2.2 Nonparallel KuramotoSivashinsky equation
This equation is commonly used as a model equation for complex fluid motion as it contains many features in a onedimensional 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 userspecified performance measures. We will use this equation for assessing the accuracy, stability and robustness of modelreduction methods in a simpler case. To this end, we consider the nonparallel version of the KuramotoSivashinsky equation in the form
(10) 
with
(11) 
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 NavierStokes equations. The expression
in the KuramotoSivashinky 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 streamwiselocalized instability mechanism acting in the recirculation bubble of cylinder flow.
The final term provides a
stabilizing hyperdiffusion which damps highwavenumber (smallscale)
structures.
In view of spatial discretization with classical finiteelement methods, we consider the auxiliary variable to render the system secondorder, so that:
(12)  
(13) 
with at and at . The second equation can be seen as a constraint on the state , reminiscent of the divergencefree constraint in the NavierStokes equations.
Semidiscretization of the above equation with finite elements yields an expression whose structure is similar to (4) for the statevariable . We obtain
(14) 
with
(15) 
as the weakform of
(16) 
and the symmetrized bilinear term:
(17) 
where, again, the pointwise discretization scheme of the nonlinear term has been chosen.
We use secondorder elements for the discretization of . For the mesh, we choose and the elements are of size , which yields 4000 elements. The timestep for the simulation is . The same semiimplicit time integration strategy as described for the NavierStokes equations is used here.
For our analysis, we choose the following constants
(18) 
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 baseflow (BF) and meanflow (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 timetrace of the variable at (figure 5(d)) and four representative snapshots from the limit cycle (figure 5(e)). We have also shown the meanflow solution corresponding to the limitcycle in figure 5(f). Coming back to the meanflow eigenvalue , we observe that, similarly to the case of cylinder flow, the amplification rate of the linearized operator around the meanflow (MF) is closer to marginality. Yet, the frequency of the mode has decreased well below the frequency of the limitcycle, which is . The eigenvector of the meanflow unstable eigenmode is now closer to the oscillations of the limitcycle. The observed discrepancies are due to the fact that the limitcycle 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 limitcycle behavior.
(a)  (b) 

s w aBF bMF c  x v aBF bMF 
(c)  (d) 
T K aTR bLC cMFTR  T W ATR BLC aTR bLC cMFTR 
(e)  (f) 
x w a b c d  x w 
3 Model reduction methodology
In this section, we develop a mathematical framework for the projectionbased model reduction of semidiscretized evolution problems of the form
(19) 
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 reducedorder model then reads
(20) 
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 reducedorder model scales as , since for each of the degrees of freedom, we have a doublesummation 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 pointwise nonlinearity (as is the case here).Based on the above approach, the building of the reducedorder 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
(21) 
The matrix contains the nonnegative eigenvalues along its diagonal, ranked in decreasing order, and the columns of provide the corresponding eigenvectors. The PODmodes, and thus the basis matrix is then formed following
(22) 
By construction, the columns of are orthonormal; the corresponding matrix satisfies the orthonormality condition The PODmodes 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 reducedorder model can easily be determined from the constant term of the fullscale equation (19). It is simply obtained by a leftmultiplication of the fullscale constant term by which is equivalent to a Galerkin projection of onto the basis spanned by the columns of i.e. the PODmodes. Mathematically, we have
(23) 
3.3 The linear term
The projection of the linear term onto the reducedorder basis yields a reduced system matrix that may be computed following:
(24) 
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 largescale 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 leftmultiplied by the matrix which is equivalent to a projection onto our orthonormal basis. We obtain the nonlinear coefficients as
(25) 
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
(26) 
and the corresponding POD modes are determined as
(27) 
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
(28) 
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 reducedorder model as
(29) 
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 leastsquares sense as above, but instead by enforcing equality at selected interpolation points. Mathematically, we premultiply the above equation by a rowselector matrix which yields
(30) 
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:
(31) 
we obtain the following representation of the nonlinear term
(32) 
Note that this expression does not require the nonlinear term to be pointwise. It is the bilinearity of the nonlinear operator that ensures an operation count of the order .
In the case of a pointwise 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 reducedorder model is:
(33) 
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.
3.5 Mathematical properties of the reducedorder 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 reducedorder model.
The kinetic energy of the reducedorder model is governed by
(34) 
If we define the symmetric and antisymmetric parts of
and as follows (with superscript denoting the
symmetric and superscript denoting the antisymmetric
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
(35) 
The energy evolution equation leads to the following tight bounds (i.e., there exist that achieve the bound):
(36) 
where stands for the Frobenius norm of a matrix or a tensor, defined as
(37) 
The first term on the righthandside of equation (36) may lead to an algebraic increase of of the form
(38) 
the second term to an exponential growth following
(39) 
and the third term to a finitetime blowup according to
(40) 
The nonlinearities that we consider in this article are of purely convective type and are therefore energypreserving (under standard boundary conditions). It is straightforward to show that
(41) 
which implies for the reducedorder model that
(42) 
If this condition holds exactly, the finitetime singularity implied by the third term in equation (36) can be avoided; if, on the other hand, the reducedorder model does not satisfy one should expect a finitetime blowup in energy for some initial condition. The condition thus constitutes an important and effective test for the robustness and longterm stability of the reducedorder model.
3.6 Quality measures of reducedorder 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 reducedorder model using a Newton method to solve:
(43) 
We evaluate the relative error between the predicted baseflow and the actual baseflow ( with the BF formulation and with the MF formulation):
(44) 
We compute the leading eigenvalues/eigenvectors of the operator obtained by linearizing the equations governing the reducedorder 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:
(45)  
(46) 
We then compute the meanflow by timeaveraging the results of a temporal simulation of the reducedorder model (20):
(47) 
In a similar way then before, we define the relative errors pertaining to the meanflow and the leading eigenvalue / eigenvector of equations (20) linearized around the meanflow :
(48)  
(49)  
(50) 
The number of unstable eigenvalues should again be equal to two. The recovery of the leastdamped 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 limitcycle case (LC), displays near marginal stability properties with a frequency close to the frequency of the limitcycle 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 reducedorder models are set to those of the largescale simulation,
(51) 
where denotes the projection of the snapshots of the unsteady simulation onto the POD modes representing the state:
(52) 
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
(53) 
and the relative mean model error by:
(54) 
and respectively, assess the error due to truncation of the projection basis to POD modes and the error due to prediction of the timeevolution 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
(55) 
With the nonlinearities representing the convective term of the NavierStokes 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 nonzero term in the model may induce a finitetime singularity for some initial conditions. We have therefore also evaluated the performance of the models when considering in the reducedordermodel instead of , that is when removing the symmetric part of the quadratic term. In this case, we denote the model error .
4 KuramotoSivashinsky equation
In this section, we assess the quality of the reducedorder models for the case of the KuramotoSivashinsky 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 timespans of the learning trajectory.
4.1.1 Bases and determined from snapshots on the transient and the limitcycle ()
(a)  (b) 
n s  n s 
(c)  (d) 
a b c x w  a b c x w 
We take snapshots from the simulation shown in figures 5(c–e) over the timespan to build representative bases and of the transient and the limitcycle. 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 nonzero over the entire region the support of the nonlinearity POD modes is far more compact (within ). This result is consistent with the pronounced dropoff 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 meanflow 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 timeoscillations 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 limitcycle ()
(a)  (b) 
n s  n s 
(c)  (d) 
a b c x w  a b c x w 
We take again snapshots from the simulation shown in figures 5(c–e) but over the timespan (not represented in the figure), which is only representative of the limitcycling behavior. The sampling time remains unchanged, and the results for the MF formulation (which is consistent with the sole knowledge of the limitcycle 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 spatiotemporal structures that are smoothly convected downstream. The statePOD mode representing the meanflow deformation has disappeared since this feature is not visible in the variable over the timespan of the learning trajectory. Only POD modes representing limitcycle 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 reducedorder models and compare them with the corresponding largescale 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 .
BF  MF  TR  LC  MFTR  
MF  
1B60  12  2  0.1  0.00  2  2  1  0.1  0.01  0.5  0.5  0.00  0.1  0.3  1  
1B50  13  2  0.5  0.00  2  2  1  0.1  0.03  3  3  0.02  0.2  0.3  3  
1B40  14  2  2  0.03  2  2  1  0.2  0.2  13  13  0.1  3  1  38  
1B30  12  2  2  10  0.3  2  2  2  1  85  85  1  3  4  49  
1B20  4  8  45  15  4  2  10  3  4  103  103  4  9  9  27  
1B10  0.3  4  45  15  29  2  28  10  19  100  100  7  75  75  81  
1M60  12  2  0.1  0.00  2  2  1  0.1  0.01  2  8  0.01  0.1  0.1  1  
1M40  14  2  3  0.04  2  2  1  0.2  0.3  119  117  0.2  4  4  69  
1M20  5  6  74  17  1  2  9  3  5  140  140  7  23  23  49  
2B6060  14  2  0.1  0.00  2  2  1  0.1  0.01  1  12  0.00  0.1  13  1  
2B6040  20  2  0.1  0.00  2  2  1  0.2  0.01  3  15  0.00  1  35  4  
2B6020  24  2  0.1  0.00  2  2  1  1  0.01  12  18  0.00  4  26  10  
2B4040  16  2  2  0.03  2  2  1  0.2  0.2  13  14  0.1  2  13  41  
2B4020  18  2  2  0.03  2  2  1  1  0.2  17  18  0.1  5  14  45  
2B2020  8  8  45  15  4  2  7  3  4  104  103  4  12  11  25  
2B1010  19  4  45  15  35  0  25  7  19  100  100  7  88  106  71  
2M6060  14  2  0.1  0.00  2  2  1  0.1  0.01  3  124  0.01  1  4  4  
3B6060  62  2  0.1  0.00  1  2  1  0.1  0.01  1  10  0.00  1  26  5  
3B6040  65  2  0.1  0.00  2  2  1  0.2  0.01  3  7  0.00  1  28  6  
3B6020  67  2  0.1  0.00  2  2  1  1  0.01  16  61  0.00  3  47  12  
3B4040  58  2  2  0.03  2  2  1  0.2  0.2  14  13  0.1  2  8  42  
3B4020  60  2  2  0.03  2  2  1  1  0.2  21  29  0.1  4  14  43  
3B2020  30  8  45  15  4  2  7  3  4  103  103  4  13  31  23  
3B1010  23  4  45  15  39  0  97  75  19  100  100  7  79  96  80  
3M6060  64  2  0.1  0.00  2  2  1  0.1  0.01  4  114  0.01  1  5  3 
All results from the simulations are summarized in table 1 for bases and determined from the transient and limitcycle behavior (fig. 6) and in table 2 for bases solely obtained from the limitcycle 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 fixedpoint and timeaveraged solutions. This is an important feature for subsequent use of the models in a flow control strategy, especially for the baseflow eigenvalue.
From table 1, which deals with the case of POD bases built with snapshots in the transient and on the limitcycle, we can draw the following conclusions :

Comparing 1B60, 2B6060 and 3B6060, 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 twobases strategies (methods 2 and 3).

Comparing 1B60 to 1M60, 2B6060 to 2M6060 and 3B6060 to 3M6060, it seems slightly better to use BF formulations than MF formulations, whichever method is used.

Comparing 1B60 to 1B50, 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 energypreserving with always high values of . However, a high value of does not necessarily prevent the model from being accurate (), see line 3B6060. With the third technique, decreasing the number of POD modes for the representation of the nonlinearities always alters the energypreserving 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 meanflow properties (four subcolumns of MF column), it is seen that both the meanflow solution and the linear properties around it are easily reproduced with low values of
Comments
There are no comments yet.