Energy-conserving hyper-reduction and temporal localization for reduced order models of the incompressible Navier-Stokes equations

A novel hyper-reduction method is proposed that conserves kinetic energy and momentum for reduced order models of the incompressible Navier-Stokes equations. The main advantage of conservation of kinetic energy is that it endows the hyper-reduced order model (hROM) with a nonlinear stability property. The new method poses the discrete empirical interpolation method (DEIM) as a minimization problem and subsequently imposes constraints to conserve kinetic energy. Two methods are proposed to improve the robustness of the new method against error accumulation: oversampling and Mahalanobis regularization. Mahalanobis regularization has the beneﬁt of not requiring additional measurement points. Furthermore, a novel method is proposed to perform energy-and momentum-conserving temporal localization with the principle interval decomposition: new interface conditions are derived such that energy and momentum are conserved for a full time-integration instead of only during separate intervals. The performance of the new energy-and momentum-conserving hyper-reduction methods and the energy-and momentum-conserving temporal localization method is analysed using three convection-dominated test cases; a shear-layer roll-up, two-dimensional homogeneous isotropic turbulence and a time-periodic inviscid ﬂow consisting of a vortex in a uniform background ﬂow. Our main ﬁnding is that energy conservation in combination with oversampling or regularization leads to a robust method with excellent long time stability properties. When any of these two ingredients is missing, accuracy and/or stability is signiﬁcantly impaired.


Introduction
Computational fluid dynamics (CFD) has become an integral part of many modern engineering applications.The increase in computational power in recent decades has allowed engineers to model increasingly larger fluid dynamical systems.However, many modern applications are of a multi-query or real-time nature e.g.design optimization [82] and uncertainty quantification [29] or realtime control [81,71] and digital twin technology [43].These applications still pose prohibitively large computational costs.Reduced order models (ROMs) have been proposed as a solution to this problem.A ROM is a type of surrogate model that approximates the high-dimensional full scale model, or full order model (FOM), in a low-dimensional way by finding approximate formulations of involved quantities or operators.The low-dimensionality of the model consequently makes the ROM significantly cheaper to evaluate than the FOM.
Traditionally in the CFD community, these ROMs have been constructed by projecting the fluid dynamics equations of interest onto low-dimensional linear spaces obtained from the proper orthogonal decomposition (POD) algorithm [47,93], using either Galerkin [87,81,95,50,64] or Petrov-Galerkin [18,38,17] projection.More recently, alternatives have also been explored.ROMs have been constructed without availability of a FOM by inferring the ROM from data using operator inference methods [76,92,12,63].Machine learning methods like convolutional autoencoders have been used in a projective sense [57,80] and also inference methods [69,62] have been applied to obtain nonlinear low-dimensional approximations.Other nonlinear dimensionality reduction methods like quadratic manifolds [10,34] and diffusion maps [94] have also been leveraged.However, the traditional POD-Galerkin methods remain powerful and this article will primarily deal with these methods and their natural extensions like the principal interval decomposition (PID) [14,49].
Nevertheless, the traditional methods possess several limitations.For turbulent (and convection-dominated) flows that are of engineering interest it is well-known that linear, projection-based ROMs suffer from stability and accuracy issues [87,33,8,5,64].Efforts have been made to solve this issue, an overview is provided by [33,87].A promising solution is structure-preservation, this entails constructing ROMs such that the underlying physics of fluid flows are respected.Especially conservation of kinetic energy is an important physical principle to uphold in a ROM with regards to stability as it bounds the norm of the solution [87,2,21,96].
Current attempts at constructing energy-conserving ROMs of fluid flow can be classified into four categories.The first category is constrained optimization projection [19,91,39,13,48]; here model reduction is cast into an optimization problem and constrained to preserve structure.A second category is formed by symmetry-preservation [87,2,21,96] where symmetries of the continuous analogues of ROM operators have been preserved resulting in the conservation of energy or entropy and associated stability.In [36,23,78,46,1,77,45] Hamiltonian physics are preserved in a low-dimensional setting, forming a third category.Finally, in [67,99,59,68,9,4] physics-informed data-driven approaches using inference and machine learning methods have been adopted establishing a fourth category.
The idea of structure preservation is elegant but does not resolve directly the well-known issue that reducing nonlinear models using POD requires intermediate lifting of the reduced representation to the high-dimensional (FOM) spaces.Thus, although the ROM is low-dimensional, the computational effort to evaluate it is still high-dimensional, defeating its purpose.Methods to overcome this problem are referred to as hyper-reduction methods [16].For sufficiently simple equations exact hyper-reduction methods exist that eliminate computational dependence on the FOM dimensions [87,3].This can be done when the underlying nonlinearity is of polynomial nature.Yet, these methods can become prohibitively expensive for larger ROMs as noted in [87].Hence, approximate hyper-reduction methods may be considered since these generally have better scaling properties.Examples of such approximate hyper-reduction methods are the empirical interpolation method (EIM) [11,37] and its discrete counterpart the discrete empirical interpolation method (DEIM) [24,26,25,11], Gauss-Newton with approximated tensors (GNAT) [20,18], energy-conserving sampling and weighting (ECSW) [38,32,31] and missing point estimation (MPE) [7].
Many existing hyper-reduction methods do not preserve the structure of the operators to which they are applied and thus stability can be lost.The field of model reduction of Hamiltonian systems has offered some solutions for this; [23] proposes a Hamiltonianconserving DEIM variant for systems with Hamiltonian functionals with non-quadratic terms; [97] improves this DEIM variant and [66] proposes a Hamiltonian-conserving DEIM variant for nonlinear Hamiltonian operators that preserves skew-symmetry.An approach similar to [97] was proposed in [73].However, the incompressible Navier-Stokes equations considered in this work have a quadratic Hamiltonian functional in the inviscid case [72,6] (kinetic energy) making the first two methods ( [23,97]) inapplicable.Furthermore, the approximate method proposed in [66] scales computationally the same as the prohibitively expensive exact method used in [87].
In this article we propose a novel energy-and momentum-conserving hyper-reduction method.The method is similar to the DEIM and gappy-POD [30] in which FOM operators are projected on low-dimensional spaces based on a small (carefully chosen) subset of spatial operator evaluations referred to as measurements.The new method is capable of robustly, accurately and efficiently dealing with convection-dominated flows described by the incompressible Navier-Stokes equations.The main idea will be to relax the DEIM interpolation requirement and use the resulting available degree of freedom to enforce energy conservation.Furthermore, we will propose two methods to further enhance the robustness of the energy-and momentum-conserving hyper-reduction method.Namely, oversampling [75,83], resembling the gappy-POD approach [30], and the new approach of regularization using the Mahalanobis distance (which, to our knowledge, has not yet been applied to the DEIM or gappy-POD).Finally, we will extend the feasibility of the hyper-reduced order model (hROM) with a newly proposed energy-and momentum-conserving localization method based on the PID.Though the methods proposed in this work pertain to the way the operators are projected based on some measurementswhich is not specific to either DEIM or gappy-POD -we will mostly refer to them as variants of the DEIM.This could be considered inaccurate as the DEIM is often characterized by the way the measurement points are chosen [24].In theory our methods can also be applied in the more general context of gappy-POD.However, our work has only been carried out and tested with DEIM-based measurement point selection procedures, hence we prefer to limit our naming to the DEIM.This article will be organized as follows.First, the governing equations will be introduced in section 2, together with an energyand momentum-conserving FOM and ROM.Following this, the energy-and momentum-conserving hyper-reduction method will be proposed in section 3.In section 4 the new energy-and momentum-conserving interval decomposition approach for convectiondominated flows will be discussed.Finally, in section 5 the performance of the proposed energy-and momentum-conserving hyperreduction and temporal localization methods will be analysed using three convection-dominated test cases.

Incompressible Navier-Stokes equations and kinetic energy conservation
In this work we consider the incompressible Navier-Stokes equations: ∇ ⋅  = 0, on a domain Ω ⊂ ℝ  , where  ∈ {2, 3}, with periodic boundary conditions.We denote with ℝ + ⊂ ℝ the set of positive real numbers including 0. In equations ( 1)-( 2)  ∶ Ω × ℝ + → ℝ  is the velocity field,  ∶ Ω × ℝ + → ℝ the modified pressure field (pressure scaled by density),  ∈ Ω,  ∈ ℝ + denotes time and  ∈ ℝ + is the kinematic viscosity.Several quantities related to the flow will also be considered, namely kinetic energy  ∶ ℝ + → ℝ + and momentum  ∶ ℝ + → ℝ  .Kinetic energy is defined as: with  evaluated at  and where: are the  2 -inner product and induced  2 -norm, respectively.Momentum is defined as: Denoting (, ) ∶= ∇ ⋅ ( ⊗ ),  ∶= Δ,  ∶= ∇ and  ∶= ∇ ⋅ , an evolution equation of the total kinetic energy can be found by differentiating (3) and substituting (1): This evolution equation can be further simplified when the following is considered.Since the convection operator is skew-adjoint given condition (2), it holds that ⟨, (, )⟩  2 = − ⟨(, ), ⟩  2 = − ⟨, (, )⟩  2 , meaning that ⟨, (, )⟩  2 = 0. Additionally, since the gradient operator and the divergence operator are each other's negated adjoints it can be stated that ⟨, ⟩  2 = − ⟨, ⟩  2 = 0 due to equation (2).Finally, the negative-definiteness of the diffusion operator can be employed to simplify the evolution equation of the total kinetic energy further to: Equation (5) implies that the kinetic energy, or the norm of the velocity field , is a monotonically decreasing quantity for periodic or homogeneous Dirichlet boundary conditions and is conserved in the inviscid case ( = 0).It will become clear soon that mimicking this property discretely will be crucial to nonlinear stability of the FOM and ROMs.
Furthermore, it follows from integrating equation (1) over the domain Ω that: for periodic boundary conditions.It is now our goal to mimic these conservation properties at the discrete and especially the reduced level.

Energy-and momentum-conserving FOM
Equation ( 1) is discretized with a finite volume method (FVM) on a staggered grid, resulting in a system of coupled ordinary differential equations (ODEs) complemented by a set of linear constraints: Here  ℎ ∶ ℝ + → ℝ  are the numerical velocity values arranged in a vector, Ω ℎ ∈ ℝ × is a diagonal matrix containing finite volume sizes,  ℎ ∶ ℝ  → ℝ  is the spatial discretization of the nonlinear convection operator,  ℎ ∶ ℝ   → ℝ  is the spatial discretization of the gradient operator,  ℎ ∶ ℝ + → ℝ   are the numerical pressure values arranged in a vector,  ℎ ∶ ℝ  → ℝ  is the spatial discretization of the diffusion operator,  ℎ ∶ ℝ  → ℝ   is the spatial discretization of the divergence operator,  =   +   is the total number of velocity unknowns and   and   are the numbers of velocity unknowns in the  and  directions respectively (for the case of two-dimensional domains) and   is the number of pressure unknowns.Relations ( 7)-( 8) are in turn complemented by a vector of initial conditions  0 =  ℎ (0) and suitable boundary conditions.The choice of boundary conditions for the flows considered in this work are periodic.
We use the discretization as described in [85,87] such that the following three important properties of the continuous operators are inherited by their spatial discretizations.First, the convection operator, which can be written in quasi-linear form as  ℎ ( ℎ ) = Cℎ ( ℎ ) ℎ (see [87] for details), is skew-symmetric: The quasi-linear form is not essential for the development of our hyper-reduction methods, but makes the proof of energy conservation straightforward.Second, the duality between gradient and divergence operators holds: Third, the diffusive operator is symmetric negative-definite, allowing the general decomposition (e.g., a Cholesky decomposition): We remark that the conservation properties discussed in this article are not specific to the method described in [85,87], any other finite volume discretization of ( 1)-( 2) satisfying the above properties can also be used.
Using these properties and the telescoping property of the FVM, discretization ( 7)-( 8) conserves discrete analogues of mass, total momentum and total kinetic energy.To show this, we first define discrete total kinetic energy and momentum on the basis of an inner product.To this end, the Ω ℎ -inner product is introduced.The Ω ℎ -inner product and its induced norm are defined by: where ⟨⋅, ⋅⟩, without subscript, denotes the standard Euclidean inner product.After deriving an evolution equation for the discrete total kinetic energy it will be clear why this is a natural choice.Using the Ω ℎ -inner product, the discrete total kinetic energy  ℎ ∶ ℝ + → ℝ + is defined as: where   ∈ ℝ  is a vector of ones at the vector indices where  ℎ () stores velocity components in direction  ∈ {1, ..., } and zeros elsewhere.
An evolution equation for discrete kinetic energy is found by temporal differentiation of  ℎ (): where || ⋅ ||, without subscript, denotes the Euclidean norm induced by the Euclidean inner product.The simplification in the second line follows from the application of properties ( 9)- (11).Indeed, evolution equation ( 12) is a discrete analogue to equation ( 5) and shows that  ℎ (), and so the norm of  ℎ (), is a monotonically decreasing quantity.The solution of the semi-discrete system of ODEs (7) can thus be bounded from above by: It is also a well-known property of the FVM, referred to as the telescoping property [98], that: This constitutes a discrete analogue to (6) and shows that discrete total momentum is also conserved.

Energy-and momentum-conserving ROM
To construct a reduced order model as in [87] we make the assumption that for any  ∈ [0,  ] there are elements   () of a lowdimensional linear subspace  ⊂ ℝ  that accurately approximate  ℎ ().Here, we denote  ∶= dim() and the low dimensionality of  implies that  ≪ .Furthermore, we consider  ⊂ ker( ℎ ) such that all   () automatically satisfy condition (8).Let Φ ∈ ℝ × be a POD basis for  that is orthonormal in the Ω ℎ -inner product, i.e.
Here  ∶ ℝ + → ℝ  are the generalized coordinates in  associated to the POD basis Φ. Substituting this approximation into the FOM (8) we obtain an overdetermined system.To obtain a solvable system, we test the overdetermined system against the POD modes in Φ and impose the Galerkin condition.This results in a pressure-free 2 ROM: where C () ∶= Φ  Cℎ (Φ)Φ and   ∶= Φ   ℎ Φ.
In order to consider energy and momentum conservation we first define the following quantities.A reduced order representation of kinetic energy is defined as: For total momentum a reduced order representation is: Noting that symmetry properties ( 9) and ( 11) are still satisfied by the reduced operators in (15), the evolution equation for   () follows by differentiating (16): where   ∶=  ℎ Φ.This equation is the ROM analogue of equation ( 5) and shows that the norm of the generalized coordinates satisfies: with equality for  = 0.
To conserve total momentum at the reduced level we require   ∈  ∀ ∈ {1, ..., }.In short, this is achieved by applying the SVD to an adapted snapshot matrix, in which the projection of the snapshot data on the   's is subtracted from the snapshots, and the resulting POD modes are subsequently enriched with the   's to form the final basis.For more details, see [87].An evolution equation for   () is then found as: ( = 0.
In the second line the telescoping property of the FOM operators was evoked.
Finally, a POD basis Φ satisfying all required criteria (span(Φ) ⊂ ker( ℎ ), Ω ℎ -orthonormality and containing all   ∕||  || Ω ℎ in its columns) is constructed from snapshot data using an altered POD method as described in [87].The FOM and the ROM can both be integrated in time with energy-conserving Runge-Kutta methods to conserve energy fully discretely (see [87,86] for further details).

Introduction to DEIM
The ROM (15) as constructed in the previous section is -dimensional, but evaluating the convective terms in a naive fashion still requires a computational effort that scales with the FOM dimension  .One approach to alleviate this issue is the exact tensor decomposition.This method results in an exact low-dimensional representation of the convection operator but scales with ( 3 ), which becomes computationally expensive for realistic values of  [87].A more efficient alternative approach, which will be considered here, is the discrete empirical interpolation method (DEIM) [24].Using the DEIM we will construct a so-called hyper-reduced order model (hROM).
The DEIM approximates  ℎ () with elements from a low-dimensional linear subspace  ⊂ ℝ  .Here, we denote  ∶= dim() and the low dimensionality of  implies that  ≪ .Let  ∈ ℝ × be an orthonormal basis for  constructed from snapshot data of the operator using the POD algorithm, where the snapshot data is gathered in a matrix Ξ ∈ ℝ ×  given as: ] .
The basis  thus satisfies ⟨  , ,  , ⟩ =   .The DEIM approximation then takes the form: where  ∈ ℝ  are the generalized coordinates of the approximation in  associated to basis  and are referred to as DEIM coordinates.The linear nature of the approximation allows the Galerkin projection Φ   to be precomputed in an offline step.
Hence, with the values of the DEIM coordinates available, the ROM can be evaluated in a completely low-dimensional fashion.
Determining a function  ∶ ℝ  → ℝ  to obtain the DEIM coordinates from the POD coordinates so that (20) holds exactly is an overdetermined problem.Furthermore, the least-squares solution to system (20), () =    ℎ (Φ), requires an expensive projection operation.To solve this, the DEIM bases its approximation on only an -dimensional subset of evaluations of  ℎ (Φ) at points on the grid, referred to as the measurement points, and denoted by .More specifically, a measurement matrix  ∈ {0, 1} × is defined that consists of selected columns of the  ×  identity matrix corresponding to the measurement points.The DEIM coordinates are obtained from an interpolation condition of which the formal solution is () = (  ) −1    ℎ (Φ).A greedy measurement point selection procedure that ensures that (  ) −1 exists and that provides near-optimal approximations of  ℎ (Φ) was proposed in [24].Using this method, the DEIM coordinates can be determined using only  evaluations of  ℎ (Φ).
We will analyse the conservation properties of this 'conventional' DEIM approach.Because the telescoping property of operators discretized using the FVM is a linear property, it is conserved under application of the POD algorithm.This may be shown by analysing the eigenvectors of the correlation matrix ΞΞ  .Using this property it can be shown that momentum is conserved for a hROM using the DEIM: = 0.
We now arrive at the key issue addressed in this article.This issue is that the hROM does not conserve reduced total kinetic energy, i.e. it can occur that: ⟨ , Φ  () ⟩ does not equal zero.This is because the DEIM-interpolated convection operator: is no longer skew-symmetric.This can lead to instability as the norm of the generalized coordinates () is not bounded from above.
An approach to retain the skew-symmetry of a DEIM approximation and therefore the energy-conservation property of the ROM has been proposed in [66].However, this method scales equivalently to the exact tensor decomposition and is therefore considered infeasible for realistic fluid dynamics applications.We will present a new approach to make the hROM energy-conserving, while also retaining the momentum conservation property.

A novel energy-and momentum-conserving DEIM
The method we propose will be named the constrained least-squares discrete empirical interpolation method (CLSDEIM).It is based on the realization that skew-symmetry of the interpolated quasi-linear convection operator ( 23) is a sufficient but not a necessary condition for energy-conservation.Instead, a necessary condition to conserve reduced total kinetic energy is: Condition ( 24) can be satisfied even if the Galerkin-projected DEIM operator Φ   ∶ ℝ  → ℝ  is not skew-symmetric and, when satisfied, results in the correct reduced kinetic energy evolution equation (17).The new idea of the CLSDEIM is to enforce this condition by posing the DEIM as a constrained minimization problem.
The conventional DEIM finds the DEIM coordinates () by minimizing the Euclidean norm between the nonlinearity and the DEIM approximation in the measurement space.The Euclidean norm can be considered as minimized since the difference in  i.e.    ℎ (Φ) −   (), is zero.The idea of the proposed CLSDEIM is to employ this view of DEIM as a minimization problem and constrain it to take place over the set  () of DEIM approximations satisfying condition (24), defined using the DEIM coordinates as: The set  () is referred to as the feasible set.As DEIM approximations with () ∈  () satisfy condition (24), the CLSDEIM produces approximations that conserve reduced total kinetic energy.The constrained minimization problem to find the DEIM coordinates () will be posed as the following linearly constrained least-squares problem: This means the CLSDEIM relaxes the interpolation condition in  between the FOM's convection operator and the DEIM approximation imposed by the conventional DEIM.Rather, the CLSDEIM minimizes the difference between the approximation and the FOM operator in  while simultaneously constraining the approximation to be energy-conserving.
The constrained minimization problem (25) is solved using the method of Lagrange multipliers [15].The Lagrangian  ∶ ℝ  ×ℝ → ℝ is defined as: where  ∈ ℝ is a Lagrange multiplier.Taking partial derivatives of the Lagrangian and setting them to zero leads to the following system for the optima Note that this matrix is symmetric and has constant coefficients with the exception of the last row and column which depend on .Solving ( 26) can be done explicitly using the inverse of a symmetric 2 × 2 block matrix [60]: where  is assumed symmetric.Taking: where  is indeed symmetric and nonsingular due to the construction of  [24], the optimal DEIM coordinates   solving (25) define the function  for the CLSDEIM coordinates as: We use the notation  ∶= () to explicitly denote the dependence of  on  and the fact that it can be treated as a vector.
Since    −1  = ()   −1 () and    −1 2(  )     ℎ (Φ) = ()  (  ) −1    ℎ (Φ) are simply scalars, the final form of the function  is: CLSDEIM: The DEIM coordinates found from equation ( 28) result in a DEIM approximation () satisfying condition (24) and consequently also equation (17).All terms in (28) are either inner products or matrix-vector products between low-dimensional vectors or matrices and vectors respectively, hence the CLSDEIM can be solved in a low-dimensional fashion.Assuming that  ≥  and that (  ) −1 and  −1 are precomputed, the most expensive operations in (28) are calculating (  ) −1    ℎ (Φ) once    ℎ (Φ) is obtained and  −1 () once () is obtained.The computational complexity of both steps scales identically as ( 2 ).Evaluating the DEIM coordinates using the CLSDEIM therefore scales computationally as ( 2 ), the same scaling as the conventional DEIM.Note that the first term on the right hand side of equation ( 28) is the solution to the DEIM interpolation condition.Hence, equation ( 28) can be interpreted as the conventional DEIM with an extra term arising from the constraint.The extra term projects the conventional DEIM coordinates on ker(⟨(), ⋅⟩) such that condition ( 24) is satisfied and minimization problem ( 25) is solved.
Considering the objective function and the geometric interpretation of the DEIM [26,24], the CLSDEIM can be interpreted geometrically as an oblique projection of  ℎ (Φ) on the subspace   ⊂ ℝ  of all DEIM approximations c with c ∈  ().The subspace   is defined as: where ⟨Φ, ⋅⟩ ∶ ℝ  → ℝ is the functional taking the Euclidean inner product with Φ.Particularly,   are all vectors in the range of  that satisfy condition (24).Since   is the intersection between two linear subspaces of ℝ  it is also a linear subspace of ℝ  .Contrary to the conventional DEIM, the CLSDEIM also projects obliquely through the measurement space .Hence, the orthogonal projection of the DEIM residual () ∈ ℝ  on , as given by: is generally not zero.However, the CLSDEIM projector Π   ∶ ℝ  →   does project  ℎ (Φ) onto   such that the orthogonally projected residual   () is minimal (with respect to the elements of  ()) in the Euclidean norm.This statement follows trivially from the formulation of the CLSDEIM minimization problem (25).The geometrical interpretation of this statement is that, in the space ℝ  , the CLSDEIM calculates the coefficients with respect to the (non-orthogonal) basis    of the orthogonal projection of    ℎ (Φ) onto the subspace ker (⟨   Φ, (  ) −1 (⋅) ⟩) ⊂ ℝ  .Finally, a proof that Π   is in fact an oblique projection operator, existence and uniqueness results for CLSDEIM minimizers and a proof of some error bounds are provided in Appendix B, Appendix A and Appendix C, respectively.The error bounds show that the CLSDEIM is also nearly optimal in the same sense as described for the conventional DEIM in [24].
Furthermore, the CLSDEIM basis  and the measurement space  can simply be found following the procedures of the conventional DEIM algorithm.Indeed, as we are using the conventional DEIM basis the reduced total momentum will also remain a conserved quantity for the CLSDEIM as was shown in equation (21).
Remark 1.Both implicit and explicit Runge-Kutta (RK) methods will be considered to integrate the proposed hROMs.To achieve exact energy conservation (for inviscid flows) energy-conserving RK methods can be used.Many families of energy-conserving RK methods are implicit [86,87], however explicit energy-conserving RK methods exist too (see [51,79] and the references therein).Using such methods and energy-conserving hyper-reduction, the change in reduced total kinetic energy   is given by: which is the time-discrete equivalent of (17).A derivation of ( 29) has been added in Appendix D. Here,  ∈ ℕ and   ∈ ℝ are parameters of the given RK method and   ∈ ℝ  are the generalized coordinates at the  th RK stage.In the inviscid case it can be seen from ( 29) that the hROMs conserve reduced kinetic energy.In the viscous case the resulting hROMs will be nonlinearly stable when   ≥ 0 ∀  ∈ {1, ..., }, which holds for example for the Gauss-Legendre RK methods (see [87,86] for further details).

Preventing overfitting errors
Interpolation-based methods like the DEIM and CLSDEIM can be prone to overfitting when noise is present [44].It was shown in [75] that adding Gaussian noise  ∈ ℝ  to the measurements, i.e.    ℎ () + , results in error bounds on the DEIM approximation of  ℎ () that grow as ( √ ).Similar overfitting errors may occur in the presence of numerical errors in  as we will show using experiments later.As a result of such overfitting errors, numerical errors in the solution may accumulate, further harming the capability of the DEIM to reconstruct the correct  ℎ (Φ) of the following time steps.This is especially notable when the interpolation procedure is carried out many times, which is the case when solving time-dependent problems.We remark that this can happen even in the case of fully-discrete energy-conserving schemes since energy conservation is a restriction on the norm of the solution and not on its accuracy.In this case, overfitting errors typically manifest themselves in non-convergence of the nonlinear solver.To address overfitting we therefore propose two further stabilization approaches.The first one is oversampling, which has recently gained significant attention [75,83].The second approach is regularization, which is a well-known technique from optimization and has thus far, to the author's knowledge, not been applied to the DEIM.This latter approach has the benefit that additional measurement points are not necessary and therefore reduces the computational burden compared to oversampling.

Oversampling
When using oversampling we replace the interpolation of the DEIM and CLSDEIM (25) with a regression procedure.This is achieved by considering a larger measurement space than the dimension of the DEIM space , thus   ∶= dim () > .As a result, the measurement matrix becomes  ∈ {0, 1} ×  .The formulation of the underlying optimization problem remains the same as optimization problem (25), but now the larger measurement matrix is used.Due to the fact that we are trying to fit a linear combination of  DEIM modes to   measurements with   > , the resulting DEIM approximation () is less sensitive to noise in measurements of  ℎ (Φ) in .This was also analysed in [75], where it was shown that oversampling using randomized measurement points resulted in a contribution of noise to the error bound of the approximation of  ℎ () that scales with  (√ ∕  ) .
Hence, overfitting errors as a result of noise can be mitigated by choosing   sufficiently large compared to .We will refer to CLSDEIM being extended with oversampling as OCLSDEIM.The DEIM coordinates as calculated using the OCLSDEIM require solving optimization problem (25) with the new measurement matrix  .The solution procedure of this problem is completely analogous to that of the CLSDEIM (25) and the final expression for the optimal   defines the OCLSDEIM function: OCLSDEIM: where instead of (  ) −1 the Moore-Penrose pseudoinverse (  ) † appears.Assuming  ≥  and that (  ) † ∈ ℝ ×  and  −1 are precomputed, the most expensive step in evaluating (30) is computing (  ) †    ℎ (Φ) once    ℎ (Φ) is obtained.As a result, calculating the DEIM coordinates using the OCLSDEIM scales as (  ).
The measurement point selection procedure will be as follows.The first  measurement points are found using the classic DEIM procedure described in [24].The next   −  points are generated from a uniform random distribution where care is taken to not select repeating indices.This procedure strikes a good balance between efficiency and accuracy and is simple to implement.Alternatives are possible, see for example the following articles for deterministic approaches to measurement point selection: [75,83,100,61].

Generalized Tikhonov regularization
When overfitting occurs, DEIM modes interpolate measurements containing noise or numerical errors resulting from time integration, negatively effecting accuracy in subsequent time steps.Although the DEIM measurement points are chosen to minimize ||(  ) −1 || thereby mitigating the effect of overfitting [24], accumulation of errors can still take place as we will see during the numerical experiments.This accumulation can eventually result in severe errors and even numerical instability.To prevent this issue we propose regularizing the least-squares problem using a generalized Tikhonov regularization [44].Using such a regularization procedure we can penalize undesired properties of solutions to the DEIM optimization problem by adding extra terms to its associated cost function.In this article we will propose to use regularization terms of the form: Here ||c|| Γ ∶= √ ⟨c, Γc⟩ with Γ ∈ ℝ × symmetric positive definite (SPD),  ∈ ℝ  is some desired or reference state and  ∈ ℝ + is a hyperparameter.The minimization problem can again be solved using the method of Lagrange multipliers.The Lagrangian  ∶ ℝ  × ℝ → ℝ of this minimization is defined as: where  ∈ ℝ is a Lagrange multiplier.Taking partial derivatives of the Lagrangian and setting them to zero leads to the following system for the optima Again, this matrix is symmetric and has constant coefficients with the exception of the last row and column which depend on . Denoting: the DEIM coordinates as calculated by the regularized CLSDEIM are given as: Regularized CLSDEIM: −1 Γ ().(32) Note that this equation also applies to a combination of oversampling and regularization.Taking  = 0 in equation ( 32) we obtain the normal CLSDEIM since  −1 Γ (  )  = (  ) † = (  ) −1 in this case.A data-driven method to calculate Γ and  that uses the knowledge we have on the DEIM coordinates of our snapshot data Ξ will be proposed in the following section.The hyperparameter tuning of  was carried out manually for simplicity.More advanced approaches to find a priori well-performing hyperparameters for regularized minimization problems like (31) are generalized cross-validation [35] and L-curve based methods [42], which we suggest implementing in future work.Assuming  ≥ , that  −1 Γ and    are precomputed and that we are combining oversampling and Tikhonov regularization, the most expensive step in evaluating (32) is calculating (  )     ℎ (Φ) once    ℎ (Φ) is obtained.Due to oversampling the computational complexity of the previous step can generally be (  ), where   =  when we are strictly using Tikhonov regularization.

Mahalanobis distance
The method we propose to determine the regularization term is to use the Mahalanobis distance   ∶ ℝ  → ℝ + [44] of the solution () with respect to the data matrix   Ξ ∈ ℝ ×  .The Mahalanobis distance is computed as: and thus corresponds to taking Γ =  −1 and  =   in optimization problem (31).Here, the matrix  ∈ ℝ × and vector   ∈ ℝ  are the covariance matrix and the mean associated to the data in   Ξ, respectively.This method is based on the idea of promoting solutions () that are similar to those observed in the data matrix   Ξ (which has the DEIM coordinates of the snapshots in Ξ as its columns).Here, similar is in the sense of originating from the same distribution.In fact, regularization using the Mahalanobis distance   corresponds to assuming a prior distribution on the optimal () that can be accurately characterized by a mean   and a covariance matrix .The Mahalanobis distance is then the multidimensional generalization of the number of standard deviations an observation is removed from the distribution's mean.

R.B. Klein and B. Sanderse
In this paper we will concern ourselves with solution reproduction problems for simplicity, and not with parameterized problems.We will see that, as numerical errors accumulate during time integration, DEIM coordinates take increasingly different values from their reference values in   Ξ when using the DEIM or CLSDEIM.Hence, such values will contribute to a large   (()) and as a result be less optimal than a more well-behaved solution given the regularized CLSDEIM loss-function (31).We assume that this idea can be extended to the parametric case, i.e. that the DEIM coordinates of a parametric problem can also be treated as generated by a single prior distribution.We will refer to the CLSDEIM regularized by the Mahalanobis distance as the MCLSDEIM.
The Mahalanobis distance   (()) and the residual term in (31) can differ by orders of magnitude.In the current work we solve this issue by taking very small hyperparameter values .We will validate the choice for such a small hyperparameter for one of our experiments in the results section by means of plotting a number of L-curves [42] and showing that the chosen hyperparameter produces solutions in the heuristically optimal L-curve corner.
Remark 2. L-curves provide a heuristic to choose good hyperparameters for (generalized) Tikhonov regularization terms in leastsquares problems [42].They are curves parameterized by the hyperparameter  ∈ ℝ + tracing the quantities (||   ℎ (Φ) −     ||, ||  − || Γ ) as in (31) with   ∈ ℝ  being the solution to (31) for a specific hyperparameter value .L-curves are typically given in log-log plots.When  becomes large L-curves for problems like (31) tend to level off as the regularization term ||  − || Γ begins to dominate the solution and the residual term ||   ℎ (Φ) −     || increases due to a bad fit to the data.This results in a horizontal region in the L-curve.When  becomes small the residual term decreases, but the regularization term increases sharply due to bad properties of the fit, resulting in a vertical region.The vertical and horizontal parts of these L-curves are connected by a corner where the heuristic indicates that the residual term and regularization term are balanced.Intuitively, this makes the range of associated hyperparameters in the L-curve corner good choices.We suggest [42,41,56] for a more detailed discussion on L-curves and Tikhonov regularization.

An energy-and momentum-conserving temporal localization method
Having constructed new energy-conserving and robust DEIM methods, we now extend the feasibility of our approach with a new energy-and momentum-conserving temporal localization method.Namely, it is well known that as flows become more convection dominated ( → 0) the Kolmogorov N-width decay becomes increasingly slower.This makes it more difficult for conventional POD to generate a low-dimensional basis set that accurately captures a significant amount of energy contained in both Ξ and the solution snapshots  ∈ ℝ ×  .To deal with this problem researchers have proposed several solutions, one of which is the application of the principal interval decomposition (PID) for both the construction of DEIM and POD spaces,  and  respectively [22,74,49,14,3,38,27].However, as far as we know energy and momentum conservation has not been taken into account in existing PID approaches.In what follows, we will propose a method to preserve energy and momentum throughout a full time integration using a hROM with temporal localization.

Temporal localization using PID
The PID decomposes the snapshot sets over  intervals in time [  ,  +1 ] and applies the POD algorithm to the individual intervals.
The idea is that by calculating modes tailored to specific intervals, the local timescales within the respective intervals are captured significantly better than by a set of modes calculated from the full set of snapshots [14,22].Based on snapshot sets: the PID provides sets of POD modes: Φ  ,   ,  ∈ {1, ..., }, (33) applicable to use at times  ∈ [  ,  +1 ] within their respective intervals.The local POD modes Φ  follow from an SVD of the local snapshot matrix   ∈ ℝ ×   , where    ∈ ℕ is the number of snapshots in the  th interval.The temporally localized DEIM measurement space   for the  th interval is determined solely based on the DEIM modes in   using the algorithm as in [24] specifically for   .Furthermore, we have   ∶= span(Φ  ) ⊂ ℝ  and   ∶= span(  ) ⊂ ℝ  .
Setting up the hROM using the PID, the dynamical system takes the form: where    ∶= Φ    ℎ Φ  .Note that we also introduce subscripts for the generalized and DEIM coordinates   and   respectively, as they are only valid during the interval [  ,  +1 ].The DEIM coordinates are calculated using any of the previously discussed methods using the appropriate measurement spaces and DEIM and POD bases.It should be noted that all the conservation properties of the energy-and momentum-conserving DEIM methods hold within intervals.This is trivial as nothing changes in the way the DEIM coordinates are calculated.Instead, the art is to develop interface conditions such that during the transition from one interval to the next, these conservation properties are still satisfied.
Note that in the subsequent sections a hROM using the PID and a certain DEIM algorithm will be referred to as PID-CLSDEIM hROM (when for example the CLSDEIM is used without regularization).
We will define how to transition between intervals based on an interface condition.Using the conventional PID [3] the following interface condition is defined: which is stable in the sense that it does not lead to an energy increase, as we will show.This condition is satisfied by the coefficients with respect to Φ +1 of the Ω ℎ -orthogonal projection of  −  onto  +1 : Note that the interface matrix   ∶= Φ  +1 Ω ℎ Φ  ∈ ℝ  +1 ×  can be precomputed in an offline stage to obtain a cheap online evaluation.The kinetic energy due to the Bessel inequality [54], which states that for the (complete) inner product space defined as and orthonormal sequence with (Φ +1 ) , ∈  Ω ℎ ,  ∈ {1, ...,  +1 } we have the following: Bessel Inequality in  Ω ℎ : Hence, application of the conventional PID using interface condition (35) cannot lead to an increase in kinetic energy for any number of intervals , making it nonlinearly stable.Furthermore, the solution is also optimal in the sense that ( 36) is also the solution to the minimization problem: In summary, the conventional PID is nonlinearly stable.However, kinetic energy and momentum are in general not exactly conserved when transitioning between intervals using condition (35), which can lead to artificial (numerical) dissipation of the solution.To solve this issue of lack of conservation, we propose the following new energy-and momentum-conserving interface condition.

Energy-and momentum-conserving interface conditions
To consider energy and momentum conservation over interfaces we define and similarly . In order to conserve energy and momentum over interfaces we benefit from the optimization formulation (37) underlying the non-conserving interface condition (35).We propose to constrain the optimization problem (37) to conserve   () and   (), leading to: The constraints on this problem can be simplified by realizing that the momentum component (  ())  in direction  is fully determined by the  th component (())  of the generalized coordinates, namely: due to Ω ℎ -orthogonality.Furthermore, the first  POD modes for every interval are the same and given by the momentum-conserving modes ]  .
Since the momentum conservation constraint fully determines the first  components of the generalized coordinates  +1 () we can simplify the optimization problem (38) by formulating it in terms of the new (⋅) variables.Note that at any interface:

𝒖
due to the momentum-conserving interface conditions.The kinetic energy constraint can also be simplified to ||ã +1 ( +1 )|| As before, we solve this constrained minimization problem using the method of Lagrange multipliers [15].The Lagrangian  ∶ ℝ  +1 − × ℝ → ℝ of this minimization is as follows: where T ∶= Φ +1 Ω ℎ Φ and  ∈ ℝ is a Lagrange multiplier.Taking partial derivatives of the Lagrangian and setting them to zero leads to the following nonlinear system for the optima (  ,   ):      − ã ( +1 )  ã ( +1 ) = 0. Rewriting the equation in the first line, an expression for the generalized coordinates at the optimum   is obtained.This expression can be substituted into the second optimality condition to solve for 1 1+  .Finally, this results in: To interpret this solution, we consider the subset () ∶= { ∈  Ω ℎ | |||| Ω ℎ = } of all vectors in  Ω ℎ with the same norm of  ∈ ℝ + (and thus the same kinetic energy) and define the feasible set of minimization problem (39): The solution then corresponds to an orthogonal projection (in  Ω ℎ ) of Φ ã ( +1 ) on span( Φ+1 ) and a subsequent rescaling to obtain an element of  .Solution (40) then produces the generalized coordinates with respect to Φ+1 of this scaled projection.A visual interpretation for  = 3 and Ω ℎ =  is provided in Fig. 1.
Using equation (40) and the momentum conservation constraint the generalized coordinates of  +  as determined by the fully constrained minimization problem (38) are given by: We will refer to using the PID with this energy-and momentum-conserving interface condition as SP-PID.
It may happen that || T ã ( +1 )|| = 0, in this case we cannot scale the result to have the same norm, and thus kinetic energy, as Φ ã ( +1 ).This situation occurs when Φ ã ( +1 ) ⟂ Ω ℎ span( Φ+1 ) and it can be shown that, if dim(span( Φ+1 )),  > 1, minimization problem (39) has an infinite number of solutions.Namely, writing out the loss function of minimization problem (39), taking into  case we have an infinite number of solutions.A visual representation of this case for  = 3 and Ω ℎ =  is given in Fig. 2. In principle any element of   can be used to determine  +  in this situation.We note that in general || T ã ( +1 )|| = 0 is an unlikely scenario since the snapshots at the end of   tend to resemble those at the beginning of  +1 .As a result,   and  +1 are typically capable of resolving snapshots close to either side of their shared interface.However, to make sure that we prevent situations where || T ã ( +1 )|| = 0 (or close to zero) we suggest to use overlapping snapshots.This entails concatenating   ∈ ℕ of the rightmost snapshots in  −1 , Ξ −1 and   ∈ ℕ of the leftmost snapshots in  +1 , Ξ +1 with   and Ξ  respectively for the construction of Φ  and   , as illustrated in Fig. 3.This will promote non-empty and larger intersections between the reduced spaces and as a result discourage || T ã ( +1 )|| = 0, because both reduced spaces   and  +1 are constructed to resolve   on either side of the interface.
In addition, the overlapping of snapshot sets tends to improve the accuracy of the transition, as will be shown in the results section.

Results
In this section we will discuss the results of three test cases to be introduced in what follows.The source code for these experiments is a custom finite volume incompressible flow solver that has been written in the C++ programming language using the Armadillo linear algebra library [89,90] and the LIS library of iterative solvers for linear systems [70].All experiments in this paper have been carried out on a simple workstation with a single Intel(R) Core(TM) i7-6700HQ CPU.

Problem set-up
The first test case concerns the shear layer roll-up (SLR) [65,52], which is a flow on a double-periodic domain [0, 2] × [0, 2] with a central band of flow in a positive coordinate direction and neighbouring bands of flow in the opposite direction.The bands are joined with a thin region of strong velocity gradients and hence strong shear forces.The initial conditions are: , > , (, , 0) =  sin(), where  =  15 determines the initial thickness of the shear layers and the parameter  = 0.05 determines the initial amplitude of an unstable perturbation in the second coordinate direction to trigger the so-called roll-up.
Using the SLR the proposed hROMs will be tested for their energy-and momentum-conserving capabilities, accuracy and computational performance.Energy and momentum conservation will be tested by considering the temporal evolution of errors in conserved quantities (momentum and kinetic energy) for the first four seconds of the inviscid SLR flow.Accuracy will be tested by considering two different errors, namely the error with respect to the FOM, which is defined as and the best-approximation error which forms a lower bound to the accuracy that can be obtained using  POD modes.Comparing   and   will give an indication on how close to optimal the hROM is.We will also analyse whether overfitting takes place by considering the temporal evolution of the Mahalanobis distance   ().To validate our choice of hyperparameter for the MCLSDEIM we will inspect a set of L-curves.Finally, the computational performance will be tested by measuring execution times associated to the offline and online phases of the hROMs and comparing these to the FOM.
Based on a grid convergence study it was concluded that a 256 × 256 grid and time step Δ = 0.01 were sufficiently fine to have a reliable FOM reference result, and these settings will be used in the numerical experiments unless otherwise mentioned.The same time step is also used to build the snapshot matrix and to simulate the ROMs.The hyperparameter of the MCLSDEIM will be set to  = 1 ⋅ 10 −14 and the hyperparameter of the OCLSDEIM will be set to   = 2.

Conservation properties
The proposed hyper-reduction methods will be tested for their energy-and momentum-conserving properties.We expect   to be conserved for the inviscid case only when using the proposed hROMs in conjunction with energy-conserving Runge-Kutta (RK) methods (see Remark 1).Thus, the proposed methods and the conventional DEIM will be tested for the case  = 0 using the (implicit) energy-conserving Gauss-Legendre 4 (GL4) method [86], and compared also to the classical (explicit) Runge-Kutta 4 (RK4) method [88].Time integration will take place until  = 4, since for longer time intervals numerical oscillations develop in the inviscid FOM simulation.These numerical oscillations do not destabilize the FOM solution due to the nonlinear stability property, but they render the simulation results inaccurate and corrupt the snapshot data sets and hence the ROM.Until  = 4 the FOM solution has a quickly decaying Kolmogorov N-width, and 8 POD and 8 DEIM modes suffice to accurately capture most of the energy in the snapshots as can be seen in Fig. 4. Note that the POD and DEIM modes are obtained from snapshots taken at every individual time step of the FOM simulation.Since the initial conditions of the SLR have a net total momentum of zero in all coordinate direction we will plot the absolute value of its components |(  )  |.We expect this value to remain zero up to machine precision for all hROMs.For a momentum-conserving semi-discretization any RK method trivially conserves momentum, hence we only use RK4 (and not GL4) in order to not clutter the results.For the kinetic energy we will plot the value: which denotes the absolute deviation of   with respect to its initial value.As kinetic energy is conserved for  = 0 by the energyand momentum-conserving hROMs, this value is also expected to remain zero up to machine precision using energy-conserving RK methods.In Fig. 5 the error in the kinetic energy is displayed.Colors are used to distinguish the different hROM methods and solid lines and dashed lines are used to distinguish RK4 from GL4 time integration.The key observation is that when using the standard DEIM algorithm (with either RK4 or GL4),   deviates significantly from its initial value, while our proposed energy-conserving hROMs exactly conserve kinetic energy up to machine precision (using GL4 time integration).When using the RK4 time integrator the energyconserving methods deviate slightly from their initial value, due to the non-energy-conserving nature of RK4.These results are in line with [87] in which similar small energy errors were observed when using high-order explicit RK methods to integrate systems with energy-conserving semi-discretizations.
In Fig. 6 the error in the momentum components is displayed.In this figure colours are used to distinguish which hROM method is applied and solid and dashed lines are used to indicate the first and second component of   , respectively.It can be observed that the theoretical results on momentum conservation hold for all hROMs including the DEIM-hROM.

Convergence behaviour
To assess the behaviour of   as , the dimension of the POD space, and , the dimension of the DEIM space, are varied we perform a convergence study.This will be done for the case of viscous flow, with  = 0.001, allowing integration of the SLR flow until  = 8.Given the small differences between RK4 and GL4 in the study of conservation properties, we only consider RK4 here for efficiency.Our primary interest is the behaviour of   as a function of , hence we will sweep this parameter for some values of .
To determine what values of  and  are necessary such that the corresponding reduced spaces  and  can accurately resolve the snapshot data we study again the singular value decay of  and Ξ.The singular value decays of the snapshot data for a FOM simulation with  = 0.001 and until  = 8 have been plotted in Fig. 7 for the data sets  and Ξ.From this figure we determine a   sweep will be performed over all values  ∈ {1, 2, ..., 60} corresponding to a relative information content (RIC) [55] ranging from 89.86% to 100%.The reduced space dimension will be taken as  ∈ {30, 40} both corresponding to RIC ≈ 100%.
The results of the convergence study have been plotted in Fig. 8 and Fig. 9 for  = 30 and  = 40 respectively.In these figures we also plot   for the non-hyper-reduced order model (dashed lines).This value will be referred to as the 'ROM error'.Clearly, the ROM error is independent of the dimension  of .Although the ROM error is not strictly a lower bound for   of hROMs, the methods should converge to it if the convection operator is approximated accurately enough.All methods show erratic error behaviour as  increases for  < .This shows that the convection operator cannot be expected to be approximated well throughout a long time-integration procedure when  is not sufficiently large compared to .When  ≥  better accuracy is obtained with the OCLSDEIM and MCLSDEIM, and the ROM error is closely approximated.This means that the convection operator is approximated well enough throughout the time integration such that Φ   ℎ (Φ) and Φ  () are nearly indistinguishable.It can be seen that   as obtained using the DEIM and CLSDEIM remains behaving erratically when  ≥  and generally the ROM error is not obtained.This shows that exactly (DEIM) or near exactly (CLSDEIM) fitting the approximation () to  ℎ (Φ) in  when   =  is not sufficient to guarantee accuracy in a long time integration process.Indeed, using more measurements (OCLSDEIM) or smarter choices of DEIM coordinates (MCLSDEIM) is required to obtain the ROM error.

Temporal error evolution
To obtain a better understanding of the performance of the methods we will analyse the full temporal evolution of several errors.We will examine a particularly bad case for the DEIM and CLSDEIM such that the errors are further exaggerated, specifically we will look at  = 30,  = 50, for which Fig. 8 shows high errors for both DEIM and CLSDEIM.For illustration, the velocity and vorticity fields as obtained using  = 30,  = 50 are shown in Fig. 10.For this case the error due to the DEIM becomes unstable and due to the CLSDEIM becomes very large.Besides   , we will also analyse the temporal evolution of the operator error: where ũ ∶ ℝ + → ℝ  is the solution to the non-hyper-reduced ROM.These errors provide insight in how well a hyper-reduction method is approximating the convection operator.
In Fig. 11 the error evolution of both   and   have been displayed for  = 30,  = 50.The ideal projection error   has also been plotted as function of time.We once more use colours to distinguish different methods, solid lines are associated to the evolution of   using the RK4 integrator and dashed lines to   also using the RK4 integrator.We can see that the errors due to DEIM and CLSDEIM grow gradually as time progresses.Furthermore, as the solution becomes more dominated by errors it becomes more difficult for the DEIM and CLSDEIM to accurately approximate the correct convection operator output.The OCLSDEIM and MCLSDEIM are much more accurate throughout the full time integration process.It can be seen that their respective   values stay close to the lower bound   and that the operator output is approximated very accurately.
Plotting the Mahalanobis distance   (()) in Fig. 12, the poor accuracy of the DEIM coordinates as determined by the DEIM and CLSDEIM is accentuated.As time progresses the values of the DEIM coordinates becomes increasingly less similar to the reference values in   Ξ, whereas the DEIM coordinates as determined by the OCLSDEIM and MCLSDEIM remain close to the reference  distribution.Here, similar is in the sense of the Mahalanobis distance.Furthermore, Fig. 12 shows that the Mahalanobis distance, at least in the reproduction case, is a good indicator of a poor choice in DEIM coordinates as it becomes increasingly larger when   becomes larger.Combining these results and the previous results on error evolution, we see that by exactly (DEIM) or near exactly (CLSDEIM) fitting the approximation () to  ℎ (Φ) in  with   = , DEIM coordinates are determined that cause gradual accumulation of errors in the time integration process.Instead, more robustness is required to deal with demanding situation like time integration where accumulation of errors can cause problems.As shown in our experiments, robustness of the DEIM and CLSDEIM can be significantly improved by: 1. performing a regression through a larger set of measurements than free parameters (OCLSDEIM); 2. using knowledge on the prior distribution of DEIM coordinates (MCLSDEIM).

MCLSDEIM hyperparameter validation
We will briefly validate the choice of  = 1 ⋅ 10 −14 for the MCLSDEIM (31) by means of L-curves.To do so, we will consider a single reduced space configuration  = 30,  = 50 and calculate the residual norm, ||   ℎ (Φ) −     ||, and the Mahalanobis distance,   (  ) = ||  −   ||  −1 , for the MCLSDEIM coefficients   associated to hyperparameter values in the range  ∈ [ 1 ⋅ 10 −17 , 9 ⋅ 10 −10 ] at all times   ∈ {1, 2, ..., 8}.We will calculate the required value of  with the accurate OCLSDEIM method also with  = 30,  = 50 and   = 100.As a result, we will have an L-curve for every point in time   ∈ {1, ..., 8} as shown in Fig. 13.To validate the choice of  = 1 ⋅ 10 −14 we specifically mark the points on the L-curves in Fig. 13 with this hyperparameter value.By inspection, the dots are roughly in the corners of the L-curves (with exception of  = 8), indicating a good choice of hyperparameter.Note that in Fig. 13    −1 should be of the same order of magnitude.
We will not repeat these experiments for the coming test-cases, but assume that if  = 1 ⋅ 10 −14 gives acceptable values for the error with respect to the FOM,   , that the hyperparameter choice is good.

Performance
To end our discussion on the SLR experiment we will show that the hyper-reduction methods exhibit similar computational performance as the conventional DEIM.This will be done by measuring the on-and offline phase wall-clock times and the associated speed-ups with respect to the FOM.Here, the offline phase is constituted by the necessary SVDs to construct  and , the construction of  and the precomputation of the necessary matrices and LU-decompositions.The online phase is constituted by the numerical integration of the hROMs.We will evaluate this for  = 30,  = 40 as all methods provide accurate results for this configuration of reduced space dimensions.To see the effect of increasing   for the OCLSDEIM we will test for   ∈ {2, 4}, where we expect to see a decline in performance as   is increased.
The results are tabulated in Table 1 for the online and offline phase.It can be seen that the computational performance is similar for all methods, both in terms of online and offline performance.The DEIM, CLSDEIM and MCLSDEIM are especially close due to the fact that for all these cases   = .The effect of the larger number of measurement evaluations required can be seen for the OCLSDEIM(2m) and OCLSDEIM(4m), which are slower than the other methods.Furthermore, using   = 2 of the OCLSDEIM indeed provides better computational performance than   = 4, as expected.

Problem set-up
The second test case is two-dimensional homogeneous isotropic turbulence (2DT).2DT is a flow on a double-periodic domain [0, 2] × [0, 2].We will generate initial conditions following the procedure outlined in [84].In short, we construct an initial energy spectrum: , where is the wave-vector magnitude.The maximum value of the initial spectrum takes place at   = 12,   ∈ ℝ is a normalization factor given by   ∶= (2 + 1) +1 ∕(2  !) and  ∈ ℝ is a shape parameter taken as  = 3.From this spectrum we follow [84] and calculate vorticity Fourier-coefficients like: Consequently, we generate a vorticity spectrum ω ∶ ℝ 2 → ℂ by providing random phases to the Fourier-coefficients using a process for which we refer to [84].Finally, we obtain a divergence-free velocity field by solving the relation: numerically.Note that this procedure is strictly for two-dimensional flows.
2DT is a fluid flow that has been studied extensively [3,28,84].Furthermore, it is a highly convection-dominated flow with many convected spatial structures and as a result it has a slowly decaying Kolmogorov N-width.For this reason the purpose of this test case is two-fold.Namely, the complexity of the flow makes it an adequate test case to see how well our proposed hyper-reduction methods will perform for more turbulent flows of engineering interest.Additionally, the slow Kolmogorov N-width decay will allow us to compare the PID and our newly proposed energy-and momentum-conserving alternative, the SP-PID, in a meaningful setting.
The FOM simulations are performed with a 1024 × 1024 grid and a time step of Δ = 0.0002.These settings are based on the observation that these also led to converged results for similar spatially second-order energy-conserving schemes in [84].Saving snapshots every time step for this test case is not feasible given the simple workstation used to perform all calculations.Hence, the simulation will be sampled such that the Nyquist-Shannon criterion is met [3]; this criterion is satisfied by sampling every Δ  = 0.005.

Hyper-reduction methods
Due to the previous results for the SLR flow, in this section we will only test the two best-performing methods, MCLSDEIM and the OCLSDEIM.We simulate the FOM until  = 4 with  = 0.001 and compare how well the resulting energy spectra () of the FOM and hROMs align at  ∈ {1, 2, 3, 4}.We will also provide a visual comparison of the vorticity fields at these time instances.To deal with the slow Kolmogorov N-width decay we will apply the SP-PID with  = 8 evenly sized intervals of 100 snapshots and for all intervals we will take overlap parameters   =   =  = 20.The singular values of all intervals and the full snapshot sets have been displayed in Fig. 14 and Fig. 15 for the POD and DEIM, respectively.As discussed in [3], it can be seen that the application of temporal localization indeed leads to a faster decay of singular values within intervals.This effectively increases the Kolmogorov N-width decay.For the reduced space dimensions we will use  = 30 and  = 40 for all intervals.Like [84], we will calculate the energy spectra using the following equations: where ψ ∶ ℝ 2 × ℝ + → ℂ are the Fourier-coefficients of the numerically calculated stream-function associated to a discrete velocity field  ℎ ().
We have provided the energy spectra for the different hyper-reduction methods in Fig. 16 and Fig. 17.Excellent agreement at all points in time can be observed for both the hyper-reduction methods.Furthermore, the results are in line with those observed    in [84] for  = 1000 in terms of spectra, where an accurate high-order pseudospectral method was used.We can see clearly that the spectrum is shifting to the left as time progresses, implying the presence of an inverse energy cascade.Furthermore, the inertial range nearly scales as ( −3 ) which should be the case as  → 0 [28].We have also shown the vorticity fields as calculated using the FOM and the different hROMs in Fig. 18.We can see that, besides the energy spectra matching between the FOM and hROMs, the vorticity fields are also almost identical.From this we conclude that the hROMs using the MCLSDEIM and OCLSDEIM and the SP-PID: 1. can accurately reconstruct the FOM results; 2. adhere to the physics underlying the 2DT flow.

Performance
We will briefly summarize the computational performance of the hROMs compared to the FOM for the 2DT flow.The execution times for both the online as offline phases of the hROMs and the full execution time of the FOM have been tabulated in Table 2.As can be seen the FOM takes a very long time to complete.This can be attributed to the naive and serial implementation of the C++ code which is not optimized for relatively large two-dimensional problems.Furthermore, the wall-clock times of the offline phases for both hROMs were heavily influenced by the large amount of data associated to storing the snapshot matrices in memory while calculations were performed.Extrapolating computation times while having less memory occupied indicate that a wall-clock time of ≈ 85 for the offline phases of both MCLSDEIM as OCLSDEIM(2m) seems more accurate for a memory efficient implementation of the 2DT case.

Temporal localization methods
We will now study in more detail the conservation properties of SP-PID at interfaces and compare its transition accuracy against the orthogonal projection of the standard PID as in equation (36).Furthermore, we will briefly investigate the impact of the overlap parameter .We will do this by simulating the FOM until  = 1 with  = 0.001 and only considering  = 2 evenly sized intervals.We will then construct the reduced spaces associated to these intervals with overlap parameters  ∈ {0, 20} and analyse the differences this makes on the accuracy of the transition at  = 0.5.We will also consider the accuracy of the conserved quantities energy and momentum, which should be exactly conserved by the SP-PID and not by the PID.However, for the PID, we expect these quantities to be better conserved for larger .We also expect the accuracy of the transition to be improved by taking larger , for both SP-PID and PID.During this experiment we will use the MCLSDEIM with  = 1 ⋅ 10 −14 and  = 40.Additionally, we will use a reduced space of  = 30.Finally, we will determine the accuracy of the transition using the transition error: The results of the experiment have been displayed in Table 3.It can be seen that the transition accuracy indeed increases as the overlap parameter is increased, since the transition error   decreases significantly.Furthermore,   using the SP-PID and the PID is nearly identical for both overlap parameter values.This shows that even in the unconstrained setting the SP-PID interface condition ( 38) is nearly optimal.As expected, the conserved quantities are exactly conserved over interfaces using the SP-PID.Increasingly smaller errors are made using the standard PID as the overlap parameter increases.We have also shown the energy spectra of the velocity fields before ( −  ) and after ( +  ) the interface in Fig. 19 and Fig. 20 for  = 0 and  = 20, respectively.We see good agreement between the new and old spectra for both methods, however for larger overlap parameters the smaller scales in the flow are captured more accurately.This can be seen in the zoomed Fig. 21 and Fig. 22 for  = 0 and  = 20, respectively.From this we conclude that: 1. the SP-PID exactly conserves both energy and momentum over interfaces, whereas the PID does not; 2. the overall accuracy of the transition (conserved quantities and   ) improves with greater overlap parameters  for the SP-PID and PID.
We suspect there will be an optimal overlap parameter, this will be the subject of future research.

Problem set-up
The third test case is a vortex in a uniform background flow adapted from [52].We will refer to the flow as AV for advected vortex.The flow is inviscid, i.e.  = 0, and takes place on a double periodic domain [−4, 4] × [−4, 4], so that kinetic energy will be conserved.The initial velocity field is given as: ) , with (, ) = √  2 +  2 and  = 0.6944.The initial velocity and vorticity fields have been displayed in Fig. 23 and Fig. 24 respectively, for later reference.The inviscid background flow will advect the vortex to the right without deformation or decay.The centre of the vortex will return at its initial position   = [0, 0]  at  = 8 after which the process will repeat.
Its time-periodic nature makes AV an ideal flow to test the stability of the proposed hyper-reduction methods and to compare them to the standard DEIM when integrating over long time intervals.Namely, if the training data contains at least one period of the flow we in principle have captured all states the flow takes as  → ∞.So, ideally, capturing one period of training data is equivalent to capturing the flow indefinitely in time.In practice, the FOM will contain some numerical error and hence will not be perfectly time-periodic, but nonetheless we expect the obtained POD and DEIM bases Φ and  obtained from one period of data to be accurate in terms of ideal projection errors for all .
We will integrate the hROMs for several reduced space sizes and beyond the training horizon of  = 8 to show that energyconserving methods have superior stability properties when integrated for long times compared to non-energy-conserving methods.
Stability will be tested by monitoring the kinetic energy error   () which should remain constant to machine precision for fullydiscretely energy-conserving methods while it can increase arbitrarily upon instability.Although the time-periodic nature of the flow is ideal to test long time stability of the hROMs, we note that AV is actually not well-suited for linear model reduction methods like those proposed in this article.This is because AV is a pure convection problem and therefore it has very slow Kolmogorov N-width decay.We will deal with this problem by applying the SP-PID, taking into account the time-periodicity of AV.
The FOM simulations are performed on a 64 × 64 grid and with a time step of Δ = 0.01.These settings do not results in numerical wiggles and are sufficient to give accurate results.The same time step is used to build the snapshot matrices and to simulate the hROMs.

Long time stability
We will first study the long time stability of the proposed energy-conserving methods by comparing them with non-energyconserving methods.To deal with the slow Kolmogorov N-width decay we will apply the SP-PID with  = 8 evenly sized intervals of 100 snapshots and for all intervals we will take overlap parameters   =   =  = 20.Note that we will also take an overlap between the first and last intervals as the system is time-periodic.The singular values of all intervals and the full snapshot sets have been displayed in Fig. 25 and Fig. 26 for the POD and DEIM, respectively.The experiment will be performed for two configuration of reduced spaces, the first configuration is  = 15 and  = 20 for all intervals, the second is  = 20 and  = 30 for all intervals.Both configurations have RIC ≈ 100% for both POD and DEIM spaces and for every interval.We will use oversampling to improve the accuracy of both the non-energy-conserving as well as the energy-conserving hyper-reduction methods, i.e. oversampled DEIM (ODEIM) [75] and OCLSDEIM, respectively.Precisely, we will take   = 100 for both reduced space configurations and all intervals.
Randomly placing measurement points for this flow is not sufficient since most points outside the vortex hardly experience any temporal changes and therefore contribute no information to the time evolution of the system at other places.Therefore, we will determine all measurement points using the DEIM algorithm [24] and  = 100 and after constructing  for all intervals truncate the DEIM space back to either  = 20, 30.We will integrate the hROMs until  = 10 = 80, far beyond the time for which the hROMs are trained.We will perform the experiment with both RK4 and GL4.
The resulting energy conservation error   as a function of time has been displayed in Fig. 27.Failure due to instability can be observed for the non-energy-conserving ODEIM for all testing configurations.While in the training domain the ODEIM has   ≈ (10 −7 ), the error sharply increases when leaving the training domain and then continues to accumulate while integrating in time until it fails catastrophically.For RK4 the kinetic energy becomes infinite, while for GL4 the Newton iterations fail to converge.It can be seen that for both RK4 and GL4 the energy-conserving hyper-reduction methods remain stable, with small deviations from energy conservation for RK4 and exact energy conservation up to machine precision for GL4.The hROMs do not only remain stable  but also remain very accurate as can be seen in Fig. 28, suffering only of a slight phase error compared to Fig. 23 and Fig. 24.These results strongly highlight the importance of energy conservation, especially when it is desired to integrate the hROMs for longer periods in time, and they stress that oversampling alone is not sufficient to guarantee stability.

Performance
We will briefly summarize the computational performance of the hROMs compared to the FOM for one period of the AV flow.The execution times of the online and offline phases of the hROMs and the execution time of the FOM, all for one AV period, have been tabulated in Table 4. Significant speed ups can observed for the hROMs using the explicit RK4 time integrator even while the FOM is computed on a rather coarse 64 × 64 computational mesh.All offline phases are comparable in duration, with those of the ODEIM being slightly shorter because less operator precomputation is necessary.The hROMs using the GL4 time integrator take significantly longer to finish due to the fact that a computationally expensive Newton iteration procedure is carried out every timestep.The execution times of the hROMs using the RK4 time integrator scale weakly with the reduced space configurations, while the execution times of the hROMs using GL4 scale relatively strongly.This scaling is a result of the Jacobian evaluations and inversions which are necessary for the Newton iterations used during the GL4 time integration.These performance results suggest the use of Jacobian-free nonlinear solvers in the future [58,53], if implicit time integration is desired.

Conclusion
In this article we have proposed new energy-and momentum-conserving hyper-reduction methods referred to as the CLSDEIM, OCLSDEIM and MCLSDEIM.In these new hyper-reduction methods the conventional DEIM is posed as a minimization problem and subsequently constrained to conserve kinetic energy, in addition to conserving momentum.The OCLSDEIM and MCLSDEIM enhance the robustness of the CLSDEIM by oversampling and by Mahalanobis regularization, respectively.
The second novelty in this work is that we have proposed a method to conserve energy and momentum over interfaces in order to use the new hyper-reduction methods in combination with the PID temporal localization method.The method is based on solving a constrained minimization problem at every interface, and its solution can be interpreted as a projection followed by a scaling to conserve energy.This new SP-PID allows the construction of fully energy-and momentum-conserving and temporally localized hROMs of convection-dominated flows.
We have performed three numerical experiments to show the performance of the new energy-and momentum-conserving hyperreduction methods and the energy-and momentum-conserving temporal localization method.We used the shear layer roll-up (SLR) to test the energy-and momentum-conserving properties, accuracy and computational performance of the proposed hyper-reduction methods.We confirmed the theoretical results that the CLSDEIM, OCLSDEIM and MCLSDEIM conserve energy and momentum in a fully discrete setting when the Gauss-Legendre 4 (GL4) Runge-Kutta (RK) time integration method is used.It was also shown that the standard DEIM caused a significant error in kinetic energy in this case.Convergence experiments showed that both the DEIM and CLSDEIM exhibit erratic behaviour as a function of the DEIM space dimension , due to error accumulation during the time integration process.This erratic behaviour and error accumulation is alleviated by the robust OCLSDEIM and MCLSDEIM.Finally, the computational performance of the DEIM, CLSDEIM and MCLSDEIM was comparable, whereas the OCLSDEIM was slower than the other methods.This shows the benefit of using regularization instead of oversampling to increase robustness of the CLSDEIM.
The second numerical experiment was the two-dimensional homogeneous isotropic turbulence (2DT).The results of this test case further highlighted the accuracy and robustness of the new hyper-reduction methods.Furthermore, the 2DT experiment showed that using the SP-PID energy and momentum can be conserved over interfaces, in contrast to the classical PID.Finally, it was shown that overlapping snapshots in the construction of local reduced spaces can increase the transition accuracy.
Finally, the third experiment was an advected vortex in a uniform background flow (AV).Due to its time-periodic nature, AV was suitable to show the superior stability of energy-conserving methods compared to non-energy-conserving methods like the oversampled DEIM (ODEIM), especially when integrating for long periods of time.The ODEIM failed due to instability after a few periods, while the energy-conserving OCLSDEIM not only remained stable but also very accurate.Furthermore, the hROMs using the RK4 time integration method were significantly faster than those using the GL4 method.
There are several directions of future research that we believe to be valuable or interesting.First of all, we have mainly concerned ourselves with simple, two-dimensional test cases due to hardware limitations.We strongly believe that the proposed methods should also be tested on larger scale test cases like those in [10,39,38].We suggest this as future research.Second, is the generalization of the MCLSDEIM and SP-PID to parametric problems.For the MCLSDEIM this entails confirming that the DEIM coordinates of a parametric problem can indeed be thought of as originating from a single prior distribution.For the SP-PID this entails finding methods to efficiently and accurately distribute the full parametric snapshot set over multiple subsets.Besides these methods, the extension to arbitrary boundary conditions (instead of periodic conditions) will be very valuable.Third, the methods proposed in this article could be merged with ideas from e.g.[23,97] to obtain energy-and momentum-conserving hROMs for problems with general first integrals and general nonlinear skew-symmetric operators (like the convection operator in this article or those in [66]).We believe this could be done by suitably altering condition (24) to conserve a general hyper-reduced first integral.Finally, we also suggest generalizing the norms used in the minimization problems of the proposed methods.For instance, using (  ) − (  ) −1 or the inverse covariance matrix of    ℎ ( ℎ ()) may prove useful.Since   ⊂  is a linear subspace of  of dimension  − 1, we can write an orthogonal basis for   as   ∈ ℝ ×−1 for some   ∈ ℝ ×−1 , where every basis vector of   is a linear combination of the basis vectors in  .Since ⟨, Φ⟩ = 0 ∀ ∈   Finally, substituting the bound for || DEIM || found in [24] and using that ||  || = 1 we obtain the desired inequality.

Fig. 4 .
Fig. 4. The normalized POD and DEIM singular value decay of snapshot data from a FOM simulation with  = 0 and  = 4.

Fig. 5 .Fig. 6 .
Fig.5.The error   in kinetic energy with respect to its initial value for the inviscid SLR using the proposed hROMs and the DEIM-hROM, the solid lines are for RK4 and the dashed lines for GL4.(For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.) R.B. Klein and B. Sanderse

Fig. 7 .
Fig. 7.The normalized POD and DEIM singular value decay of snapshot data from a FOM simulation with  = 0.001 and  = 8.

Fig. 8 .
Fig. 8. Convergence study of   evaluated at  = 8 for different hROMs and different values of  for  = 30, including the non-hyper-reduced model results (dashed).

Fig. 9 .
Fig. 9. Convergence study of   evaluated at  = 8 for different hROMs and different values of  for  = 40, including the non-hyper-reduced model results (dashed).

Fig. 11 .
Fig.11.The evolution of   (solid lines) and   (dashed lines) for  = 30 and  = 50 using different hROMs and using RK4, compared to the projection errors   (purple solid line) and   (purple dashed line).

Fig. 14 .
Fig. 14.Singular value decay of solution snapshots of intervals versus the full snapshot set.

Fig. 15 .
Fig. 15.Singular value decay of operator snapshots of intervals versus the full snapshot set.

Fig. 25 .
Fig. 25.Singular value decay of solution snapshots of intervals versus the full snapshot set.

Fig. 26 .
Fig. 26.Singular value decay of operator snapshots of intervals versus the full snapshot set.

Fig. 27 .
Fig.27.The error   in kinetic energy with respect to its initial value for AV using the proposed hROMs and the ODEIM-hROM and the different reduced space configurations indicated as (, ).Note that results using RK4 are given with dashed lines and GL4 with solid lines.

Fig. 28 .
Fig. 28.The velocity and vorticity fields of AV at  = 80 calculated using the OCLSDEIM using different reduced space configurations indicated as (, ) and time integrators.

Table 2
Computational performance (offline and online) of the different hyper-reduction methods for  = 30,  = 40.

Table 3
Transition errors in velocity field and conserved quantities at the interface on  = 0.5.

Table 4
Computational performance (offline and online) of the different hyper-reduction methods and different reduced space configurations.