{'fraction_non_alphanumeric': 0.033786189450014055, 'fraction_numerical': 0.01229270120865736, 'mean_word_length': 4.460554589174255, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Entropy stable reduced order modeling of nonlinear conservation laws

Reduced order models of nonlinear conservation laws in fluid dynamics do not typically inherit stability properties of the full order model. We introduce projection-based hyper-reduced models of nonlinear conservation laws which are globally conservative and inherit a semi-discrete entropy inequality independently of the choice of basis and choice of parameters.

directly construct entropy stable ROMs for nonlinear conservation laws by combining hyper-reduction with a modified Galerkin projection of an appropriate full order model. The approach taken in this paper combines techniques from entropy stable finite volume schemes [34,35] and entropy stable summation-by-parts (SBP) discretizations [36,37,38,39,40,41] to produce entropy stable reduced order models.
The paper is organized as follows. Section 2 introduces a full order model on 1D periodic domains based on entropy stable finite volume schemes. Section 3 describes how to construct an entropy stable reduced basis approximation, while Section 4 discusses entropy stable hyper-reduction techniques to reduce costs associated with the evaluation of nonlinear terms. Section 5 describes how to extend the aforementioned approaches to non-periodic boundary conditions, and Section 6 describes the extension to higher dimensions. We conclude in Section 7 with numerical experiments which verifying the presented theoretical results.

The full order model: entropy stable finite volume schemes
We briefly summarize entropy inequalities associated with systems of nonlinear conservation laws. Let Ω denote some domain with boundary ∂Ω. Nonlinear conservation laws are expressed as a system of partial differential equations (PDEs) where u ∈ R n are the conservative variables, f i are nonlinear fluxes, and v(u) are the entropy variables. By multiplying (1) by the entropy variables, viscosity solutions [42,43] of popular fluid systems (e.g., shallow water, compressible Euler and Navier-Stokes [44,38]) can be shown to satisfy where n i denotes the ith component of the outward normal vector. The entropy inequality (2) is the analogue of energy stability for nonlinear conservation laws [45,46]. However, due to the use of inexact treatment of nonlinear terms (e.g., from collocation approximations or inexact quadrature), most numerical methods for nonlinear conservation laws do not satisfy a discrete analogue of this stability condition. The reduced order models in this work are constructed using full order models based on entropy stable finite volume methods [34], which reproduce a discrete version of the continuous entropy inequality (2). For simplicity, we illustrate the construction of full and reduced order models on a periodic 1D domain. We note that, while the reduced order model will be constructed from an entropy stable scheme, the solution snapshots do not need to be generated byentropy stable full order model. Throughout this work, we will assume that the both full and reduced models yield "physically relevant" solutions, such that the entropy is a convex function. For example, for the compressible Euler and Navier-Stokes equations, we will assume that the thermodynamic variables (density and energy/temperature/pressure) are positive. Guaranteeing both positivity and accuracy of numerical solutions for general discretizations remains challenging [47,48], and will be explored in future work.
A key ingredient of the full and reduced order models considered in this work are entropy conservative finite volume numerical fluxes. Let u L , u R denote left and right states. Then, a two-point numerical flux f S (u L , u R ) is entropy conservative if it satisfies the following three conditions These fluxes are used to construct entropy conservative and entropy stable finite volume schemes [34, 66? ]. Let the domain be decomposed into K cells of size ∆x, and let (u h ) i denote the vector containing mean values of the vector of conservation variables over the ith cell. An entropy conservative finite volume method results from discretizing the integral form of the conservation law as follows where periodicity is imposed through the equations for (u h ) 1 , (u h ) K . We first rewrite the system (4) in a matrix form which is more amenable to model reduction. We define the skew-symmetric differentiation matrix Q and flux matrix F such that Note that Q1 = 0, that the matrix F is symmetric (due to symmetry of f S ), and that the diagonal of F is equal to f (u) due to flux consistency. For a scalar nonlinear conservation law, the matrix-based formulation of (4) is then equivalent to where • denotes the Hadamard product, 1 denotes the vector of all ones, and each entry of u h corresponds to a point x i in the domain. In this setting, Q • F extracts and takes linear combinations of nonlinear flux interactions between different nodal values of (u h ) i and (u h ) j . The discretization of the nonlinear flux term using (Q • F ) 1 is commonly referred to as flux differencing [37,49,38,39,40]. 1 For a system of n nonlinear conservation laws, all matrices are treated in a Kronecker product fashion. Let the vector u h now correspond to the concatenated vector of components of the solution u h = [u 1 h , . . . , u n h ] T . Then, the matrix formulation of (4) is given by Here, F is a block diagonal matrix where the ith block corresponds to the evaluation of the ith component of the numerical flux f i,S (u L , u R ). From this point on, we drop explicit references to components and the Kronecker product to simplify notation, so that (6) applies to both scalar equations and systems. We next show that the full order model satisfies a discrete conservation of entropy.
Theorem 1. Let u h (t) be a solution of (2) for which the entropy S(u h ) is convex at each time t. Then, u h satisfies the semi-discrete conservation of entropy Proof. The proof can be found in the literature [34,35,50]. Inspired by entropy stable summation by parts (SBP) schemes [37,49,38,39,40], we present an alternative proof which relies only on matrix properties of Q. We test (6) with the vector of entropy variables Assuming continuity in time and using the definition of the entropy variables in (1), Using skew-symmetry of Q, Rearranging indices and using the symmetry of f S (u L , u R ) = f S (u R , u L ) exposes the entropy conservation condition in the sum where we have used that Q = −Q T and Q1 = 0 in the final step.

Viscosity and entropy dissipation
Theorem 1 shows that the formulation (2) preserves a semi-discrete conservation of entropy. However, in the presence of shock discontinuities, entropy should be dissipated instead of conserved. To mimic this at the semi-discrete level, we add entropy dissipation to (6) through appropriate viscosity terms. For example, it was shown in [51] that a centered approximation of the Navier-Stokes viscosity is entropy dissipative. In this work, we apply a simple Laplacian artificial viscosity to each component of the solution [52] ∆x where is the visosity coefficient. This choice is intended to simplify the presentation of entropy stable treatments of diffusion terms. Future work will analyze both physically relevant viscosities and more nuanced artificial dissipation mechanisms.
If v T h Ku h ≥ 0, then the solution satisfies a discrete dissipation of entropy The proof is similar to the one given in [51]. We rewrite the viscous term as Since v h and u h are related through a differentiable invertible mapping, the mean value theorem implies that where ∂u ∂v i,i−1 is the Jacobian matrix evaluated at some state between (u h ) i and (u h ) j . Since the Jacobian is related to the Hessian of the entropy through ∂ 2 S(u) ∂u 2 = ∂v ∂u = ∂u ∂v −1 , it is positive definite so long as S(u) is convex. Substituting these expressions into the viscous terms and using positive definiteness of ∂u ∂v yields 3. An entropy conservative reduced order model We can now formulate a reduced basis approximation on periodic domains using the entropy stable full order models described in the previous section. Boundary conditions will be addressed later in Section 5. We also note that this reduced order model is not practical, since the cost of evaluating nonlinear terms scales with the size of the full order model. Later sections will discuss how to reduce the cost of nonlinear evaluations using appropriate hyper-reduction techniques.
Let {φ j (x)} N j=1 denote a reduced basis for each component of the solution, which may be generated using (for example) principal orthogonal decomposition (POD) or the reduced basis procedure. Let V denote the generalized Vandermonde matrix whose columns contain evaluations of φ j at grid points x i We approximate grid values of the solution by u h = V u N , where u N denote "modal" coefficients of the solution in the reduced basis. Ignoring viscosity terms for now, plugging this expression into (6) and enforcing that the residual is orthogonal to all columns of V (Galerkin projection) yields a reduced system We will show that, while this formulation is not entropy conservative, a slight modification recovers semidiscrete entropy conservation. We motivate this modification by first considering the time derivative term. Because of the Galerkin projection, we can no longer directly test with the vector of entropy variables, which may not lie in the span of the reduced basis functions. However, we can test with an appropriate projection of the entropy variables. Let v N denote coefficients of the projection of the entropy variables where V † is the pseudoinverse of V . Then, we can recover the time derivative of the discrete entropy by testing with the coefficients

The remainder of the proof requires showing that
where we have introduced the grid values of the projected entropy variables as v. We can repeat the steps of the proof of Theorem 1 up to the point when we invoke the entropy conservation condition of (3) This is not possible due to the fact that the projected entropy variables v are no longer mappings of the grid values of the conservative variables u h = V u N . To remedy this, we follow approaches taken in [53,40] and replace the grid values of u h used to evaluate the flux matrix F with values of the entropy-projected conservative variables u The entropy projected conservative variables u are thus mappings of the projected entropy variables. We note that V V † can also be interpreted as a discrete approximation of the L 2 projection operator onto the span of the reduced basis φ j . One can then construct a semi-discretely entropy conservative reduced model.
Then, the solution satisfies the following semi-discrete conservation of entropy Additionally, if 1 lies within the range of the reduced basis matrix V (e.g., 1 is in the span of the reduced basis functions φ 1 , . . . , φ N ), solutions of (11) conserve global averages of the conservative variables.
Proof. The conservation of entropy follows from testing with v N and applying (10). The remaining steps are identical to those of Theorem 1. The global conservation results from testing with 1. If 1 is in the range of V , then the exist coefficients in the reduced basis matrix e such that V e = 1. Then, Since Q is skew-symmetric and F is symmetric, Q • F is skew-symmetric and which is a discrete statement of conservation.
Remark. In order for the formulation (11) to remain accurate, the entropy projected conservative variables must accurately approximate the conservative variables. This requires that the reduced basis accurately approximates the entropy variables. This can be taken into account, for example, by computing the POD basis vectors from snapshots of both the conservative and entropy variables.

Entropy conservative hyper-reduction
While the Galerkin projected formulation (11) is entropy stable, the computational cost involved in solving the semi-discrete system scales with the size of the full order model rather than the dimension of the reduced basis. For example, explicit time-stepping methods require the evaluation of the nonlinear term V T (Q • F ) 1. Since the matrices Q and F are the original full order model matrices, the cost of the reduced system is not lower than the cost of the full order model.
To reduce the cost of evaluating nonlinear terms, we introduce a second hyper-reduction step [54]. Here, terms involving vectors of nonlinear function evaluations are approximated by terms involving nonlinear function evaluations at a subset of points [55,56,57,13,58,33,59]. In this work, we utilize a sampling and weighting strategy [13,33,59] which enables proofs of discrete entropy stability. These hyper-reduction techniques can be interpreted as reduced quadratures, and produce approximations of the form Here, g(u) denotes a nonlinear function, I denotes a subset of N s row indices corresponding to sampled points, V (I, :) denotes the sub-matrix consisting of the N s sampled rows of V , and W = diag(w) is a N s × N s diagonal matrix of positive weights. We will describe algorithms for computing hyper-reduced points and weights in Section 4.2.
We will hyper-reduce the nonlinear term (Q • F ) by constructing a hyper-reduced sampled matrix Q s . Recall that the proof of entropy conservation in Theorem 1 required that matrix Q is skew symmetric and has zero row sums. If the hyper-reduced matrix Q s also satisfies the properties then the proof of entropy conservation in Theorem 2 will carry over to the hyper-reduced model. Unfortunately, common hyper-reduction techniques preserve either skew-symmetry or zero row sums, but not both. For example, Q can be decomposed into the sum of local skew-symmetric matrices. The hyper-reduction techniques of [13,60] approximate Q using a sparse linear combination of those local matrices. This hyperreduction technique preserves skew-symmetry, but does not necessarily preserve the zero row sum property.
Similarly, techniques such as gappy POD or empirical interpolation [61,55,57] may preserve the zero row sum property but not skew-symmetry.

Two-step hyper-reduction: compress and project
We apply a two-step hyper-reduction procedure which preserves both the zero row sums and skewsymmetry of the hyper-reduced matrix Q s . Rather than directly hyper-reducing the full order matrix Q, we first construct a compressed "modal" matrix, then combine this with a projection operator based on the hyper-reduced points. In the first step, we construct a compressed intermediate operator by combining Galerkin projection with the "expanded basis" approach of [33]. Let V t denote a test basis through which to extract the action of Q. We define the intermediate reduced operator Let c be some vector in the span of the test basis such that c = V t c.
Here, d are the coefficients of the orthogonal projection of Qc onto the test basis. Note that Q t is skewsymmetric by construction. Moreover, if the vector of all ones lies in the span of the test basis (e.g., 1 = V t e for some coefficient vector e), then Q t e = 0, since Q t exactly recovers the action of Q1 = 0 and 0 lies within the span of any basis.
Since Q t acts on coefficients of the test basis (a "modal" operator), it cannot be directly applied to evaluate nonlinear functions of the solution. However, we can construct "nodal" operator by combining Q t with appropriate mappings from hyper-reduced points to coefficients of the test basis. We make the following assumptions on the hyper-reduced points to enable the construction of such mappings: Assumption 1. Let I denote the index set of hyper-reduced points. We assume that the hyper-reduced test mass matrix V t (I, :) Note that if Assumption 1 holds, then the mass matrix V (I, :) T W V (I, :) is also non-singular, since the column space of V is contained in the column space of V t .
If Assumption 1 holds, we can define the projection matrix P t onto the test basis as follows: The projection matrix P t maps values at hyper-reduced points to coefficients in the test basis, and is welldefined under Assumption 1. The key property of P t which we will utilize is as follows: suppose that f = V t c for some coefficients c. Then, In other words, if a vector lies in the range of V t , the projection matrix P t recovers the coefficients exactly.
We also note that if the hyper-reduced set I is taken to include all the points and W is a multiple of the identity (for example, W = ∆xI), then P t reduces to the pseudo-inverse We can now construct a "nodal" differentiation operator which satisfies both skew-symmetry and zero row sum properties: Then, Q t is skew-symmetric and satisfies Q t 1 = 0.
Proof. Since Q is skew-symmetric, the matrix Q t is skew-symmetric by construction as well. By (14), if 1 is exactly representable in the test basis (such that 1 = V t (I, :) e for some coefficient vector e ), then P t 1 = e. Then, Q t 1 = P T t Q t e = 0, since Q t exactly differentiates elements of the test basis.
Remark. The projection matrix P t can also be replaced, for example, by the gappy POD projection (V t (I, :)) † [62,63]. We do not observe significant differences in numerical results when doing so. In principle, any matrix operator which maps from nodal values to modal coefficients and reproduces the test basis exactly can be used for P t . In this work, we restrict ourselves to the projection matrix defined by (13), and will investigate other definitions of P t in future work.

Choice of test basis
We have not yet specified the choice of test basis V t . The span of the test basis should include the reduced basis V , and Lemma 1 implies that the span of the test basis should also contain 1. However, these two conditions are insufficient to ensure accuracy for solutions of nonlinear conservation laws. Consider, for example, the 1D periodic Burgers' equation with initial condition The solution u(x, t) is a decaying stationary shock which is anti-symmetric across the origin at all times t. The POD modes are thus also anti-symmetric across the origin; however, their derivatives QV are nearly symmetric across the origin, as shown in Figure 1. Since the entries of V T QV are inner products of nearly symmetric and anti-symmetric functions, V T QV ≈ 0, and the resulting hyper-reduced model is highly inaccurate. This example implies that the reduced basis matrix V can do a poor job of sampling the range of QV , and that the test basis should contain additional vectors. To remedy this issue, we borrow techniques from least squares Petrov-Galerkin ROMs [64,22] and enrich the test basis V t with vectors spanning the range of QV . The test basisnow spans a space which contains 1, V , and QV , and thus has a dimension of at most 2N + 1. The dimension of V t may also be smaller if the intersection of the range of the reduced basis matrix V and the range of QV and 1 is non-empty.
We comment on why this choice of test basis restores accuracy. In the absence of hyper-reduction (where , the matrix Q t should be an accurate approximation to Q. The following lemma quantifies the relationship between accuracy in the 2-norm and choice of test space.
Thus, e T 2 e 2 = e 2 2 gives In other words, if QV is contained within the range of V t , the error in approximating Qf using the

Hyper-reduction algorithm and target space
In this work, we utilize the greedy algorithm for computing empirical cubature points and weights from [32,33], which is described in Algorithm 1. Let V target be some matrix whose columns span some target space of functions to be integrated, and let w target denote some reference target weights. The algorithm approximates V T target w target by V target (I, :) T w. The index set is constructed in a greedy fashion by selecting the row index of V target which is most positively parallel to the residual, then computing the weight vector which minimizes the residual through a (non-negative) least squares solve. The iteration terminates when the norm of the residual is smaller than some user-defined tolerance tol. In this work, we utilize an N -mode POD basis, and tol is determined by the singular values of the associated snapshot matrix. Suppose there are M total singular values of the snapshot matrix; then, we define tol as Add the new index to the set of hyper-reduced points I = I ∪ {i}.

7:
Compute w using linear least squares if min w ≤ 0 then 9: Recompute w using non-negative least squares arg min Recompute r = b − V target (I, :) T w.

12: end while
Algorithm 1 requires both V target and w target as inputs. For all numerical experiments, we set w target = ∆x1 and set V target to be the span of products of POD functions V , such that This choice of target space ensures that the mass matrix ∆xV T V is accurately approximated. In practice, since many of the pointwise vector products V (:, i) • V (:, j) are linearly dependent, we replace V target by a dimensionally reduced matrix [33]. In this work, we compute this second reduction using another application of the SVD. Let µ 1 , . . . , µ Nt denote the singular values of V target , and define the energy residual of the first i modes Note that E 1 ≥ E 2 ≥ . . . ≥ E Nt . We replace the matrix V target with its k leading left singular vectors, where k is the smallest index such that E k ≤ tol.

Conditioning of the test mass matrix
We note that the choice of target space is designed to accurately approximate the mass matrix ∆xV T V , and does not take the test basis matrix V t into account. This choice is motivated by the empirical observation that incorporating the test basis by approximating either ∆xV T V t or ∆xV T t V t produces a significantly larger number of hyper-reduced points without any significant improvement in accuracy. However, because the hyper-reduction strategy does not account for V t , it is possible to construct a hyper-reduced quadrature which accurately integrates the mass matrix but violates Assumption 1 by producing a singular test mass matrix V t (I, :) T W V t (I, :). The challenge is then to construct a hyper-reduced volume quadrature which accurately approximates the reduced basis mass matrix ∆xV T V while ensuring that the hyper-reduced test mass matrix is non-singular. We address this challenge by computing the spectra of the test mass matrix and adding additional "stabilizing" points if the condition number is larger than a specified tolerance.
Let z j denote the N z eigenvectors corresponding to small eigenvalues of the hyper-reduced test mass matrix. The entries of z j correspond to vectors of coefficients in the test basis V t . Let Z denote the matrix whose columns are point values of each eigenvector To ensure that the hyper-reduced test mass matrix is non-singular, we will employ a greedy heuristic based on adding additional points to approximate ∆xZ T Z. Suppose that we have already computed hyper-reduced points and weights to approximate the reduced mass matrix ∆xV T V . Let I V denote the index set of these stabilizing points. We compute a second group of hyper-reduced points with index set I Z using the greedy hyper-reduction algorithm applied to the target space Z target , which we define analogously to V target in (16). We note that, unlike the hyper-reduced points computed previously, we do not perform a dimensionality reduction on Z target when computing this set of stabilizing points.
Given I Z , we update the hyper-reduced index set as the union I = I V ∪ I Z , and compute a new set of hyper-reduced weights through the non-negative least squares solve where α Z > 0 is some scaling parameter which controls whether the weights w more accurately approximate the integration of the target space V target for the mass matrix, or Z target for the null space of the test mass matrix. We set α Z = 10 −2 in all experiments, which is empirically observed to control the condition number of the test mass matrix without significantly impacting the accuracy to which the reduced mass matrix ∆xV T V is approximated. For the numerical experiments in this paper, it was sufficient to add one single set of stabilizing hyperreduced points. However, this procedure can be repeated multiple times to further improve the condition number of the hyper-reduced test mass matrix.

An entropy conservative hyper-reduced model on periodic domains
To summarize, an entropy conservative reduced order model on periodic domainsrequires the following offline steps: 1. Compute a reduced basis matrix V from snapshots of both conservative variables and entropy variables. 2. Use the "test" basis matrix V t to construct a "modal" differentiation matrix Q t = V T t QV t . 3. Compute a hyper-reduced quadrature using Algorithm 1, adding stabilizing points as necessary to ensure that the test mass matrix V t (I, :) T W V t (I, :) is non-singular.
4. Construct the hyper-reduced nodal differentiation matrix Q t = P T t Q t P t using the projection matrix P t = V t (I, :) T W V t (I, :) Then, we have the following theorem: Theorem 3. Suppose that the hyper-reduced quadrature satisfies Assumption 1. Then, the semi-discrete formulation semi-discretely conserves the sampled and weighted average entropy Additionally, if the reduced basis V exactly represents 1, solutions of (18) conserve global averages of the conservative variables.
where we have used that W is diagonal and the temporal chain rule in the final step. Since the hyperreduction preserves the skew-symmetry and zero row sums of Q t , the remainder of the proof consists of showing that v T N V (I, :) T (Q t • F ) 1 = 0, and is identical in structure to the proof of Theorem 2.
Remark. It is not strictly necessary to use the hyper-reduced mass matrix to construct an entropy conservative ROM. For example, we can use the exact mass matrix M N = ∆xV T V rather than the hyper-reduced mass matrix. If the projection matrix P is defined accordingly, the resulting ROM remains entropy conservative. However, we observe significantly larger errors when using the exact mass matrix instead of the hyper-reduced mass matrix.

Online cost estimates for the hyper-reduced system
In practice, the flux matrix F is not computed explicitly, but is instead formed on the fly by evaluating the entropy stable flux f S . Due to the complexity of the formulas for f S , the computation of F constitutes a significant portion of the computational cost. For the full order model, this cost is typically ameliorated by exploiting the sparsity of Q. However, since the hyper-reduced matrix Q t is fully dense, the nonlinear convective term Q•F in (18) requires O(N 2 s ) nonlinear flux evaluations. In contrast, typical hyper-reduction approaches target a nonlinear vector, which requires only O(N s ) nonlinear function evaluations. We emphasize that this O(N 2 s ) online cost is specific only to the convective nonlinearity, and not to the hyper-reduced treatment of other nonlinear terms such as viscosity.
While the O(N 2 s ) nonlinear flux evaluations are trivially parallelizable, they constitute a significant increase in online costs compared to traditional hyper-reduction procedures. However, we note that the cost can be comparable to techniques for exactly treating higher order polynomial nonlinearities for a large number of modes [30]. For example, suppose one wishes to compute V T g (V u N ), where g(x) = x m for m ≥ 1 and m integer. We can rewrite this as Suppose first that m = 2, such that the nonlinearity is quadratic. Define V j = diag (V (:, j)) as the diagonal matrix with the jth column of V on the diagonal. Then, diag (V u N ) = N j=1 V j (u N ) j , and The matrices V T V j V can be precomputed, such that the cost of evaluating the above expression is N matrix-vector multiplications, or O(N 3 ) operations. For m = 3, repeating the same steps yields Assuming flux evaluations dominate costs, entropy stable schemes require α 2 βN 2 operations. Thus, rather than increasing the asymptotic computational complexity, the entropy stable hyper-reduction presented here increases the associated constant.

Hyper-reduction and entropy dissipation
Given an entropy conservative ROM with hyper-reduction, we can construct an entropy stable ROM by adding appropriate entropy-dissipative viscosity terms. Recall that the artificial diffusion matrix K in (7) does not depend nonlinearly on u. Thus, typical model reduction techniques compress such operators using Galerkin projection, e.g., V T KV . However, since one cannot in general show that this operator dicretely dissipates entropy, we discuss two treatments of the viscous terms which are provably entropy dissipative.
This section introduces hyper-reduction procedures for K which provide a provable discrete dissipation of entropy. We seek to mimic the entropy balance of the full order model in (8). Suppose the reduced order model is given by where d(u) is some hyper-reduced approximation to the viscous terms. Theorem 3 gives that the entropy balance of the reduced order model is We will design hyper-reduced treatments of the viscous terms such that v T N d(u N ) ≥ 0, which implies that the reduced model is entropy stable.
Let D denote the (K − 1) × K difference matrix The artificial diffusion matrix in (7) can be decomposed as K = D T D. To construct hyper-reduced diffusion terms, we first compute (using the SVD) an auxiliary basis V D for DV , where V is the reduced basis matrix. We then perform hyper-reduction to compute weights and sampled row indices such that We can now construct entropy-dissipative diffusion terms using this second hyper-reduction of D. The first treatment is relatively straightforward, and simply samples rows of D corresponding to indices I D . The viscous terms Ku in the full order model are approximated using the sparse hyper-reduced matrix D (I D , :) Here, the entropy-projected conservative variables u are mappings of grid values of the projected entropy variables V v N . Repeating the steps used to derive (9) yields that the entropy dissipation of (19) is The second treatment of the diffusion terms is motivated by approaches taken in [37,38,68,52]. Recall that, using the mean value theorem, the ith entry of Du is where ∂u ∂v i,i+1 is ∂u ∂v evaluated at some intermediate state between u i and u i+1 . This motivates a hyperreduction based on sampling and weighting rows of D applied to the vector of entropy variables v(u). The diffusion terms can be approximated by where v N are coefficients of the projected entropy variables in (18) This implies that the dissipation of entropy produced by the hyper-reduced approximation of the viscous terms (21) also mimics the dissipation of entropy (9) derived for the full order model. The difference between the hyper-reduction treatments of viscosity (19) and (21) are that (19) relies explicitly on the fact that D is a two-point finite volume difference matrix, which enables the use of the mean value theorem in deriving a statement of entropy dissipation. The first approach (19) results in a simpler online phase, though it involves larger matrices compared to the second approach. The second approach (21) is not limited to finite volume methods, and is entropy dissipative for any choice of differentiation matrix D. Both approaches require computing an additional set of hyper-reduced points during the offline phase.
Remark. The hyper-reduction treatments of viscosity described in this section can also be applied to nonlinear interface dissipation models, such as local Lax-Friedrichs penalization. We do not focus on such approximations for this paper, as the amount of dissipation applied is typically proportional to ∆x and thus changes based on the dimension of the full order model.
We note that, in practice, a naive approximation of diffusive terms Ku by is often effective. This approach is simpler to implement than (19), (21), and does not require a second set of hyper-reduced points. We are unable to show that this treatment is provably entropy dissipative. However, in all experiments, we observe that v T N V T KV P u ≥ 0, which implies that entropy is discretely dissipated.

Weak imposition of boundary conditions
Until now, we have assumed periodic domains. In this section, we discuss how to extend the construction of entropy stable ROMs to the non-periodic case. Boundary conditions are weakly imposed through a numerical flux and a "hybridized" summation by parts operator [40].
In the non-periodic case, the matrix Q is nearly skew-symmetric, such that it satisfies a summation by parts property A non-periodic full order model can be expressed as where f * = [f * 0 , 0, . . . , 0, f * K ] T is a vector containing values of the nonlinear flux at the boundary points. Following [38], a boundary numerical flux f * is entropy stable if The flux is referred to as entropy conservative if the inequality in (25) is an equality. If an entropy stable or entropy conservative flux f S (u L , u R ) is used to evaluate both the volume and surface flux, then the resulting full order model given by (24) is also entropy stable or entropy conservative [34,69,38]. Extending entropy conservative treatments of boundary conditions to the reduced order model is less straightforward due to fact that reduced matrices do not satisfy the same SBP property (23). If the full order matrix Q satisfies the SBP property, then the hyper-reduced matrix Q t satisfies a generalized SBP property [70] Q where the matrices B and E encode scaling by normals and evaluation at boundary points. In 1D, these matrices are given by where V bt interpolates test basis functions to boundary points, and V t (:, 1), V t (:, K) denote the first and last columns (corresponding to points on the boundary) of the test basis matrix V t . For linear problems, it is possible to impose energy stable boundary conditions by combining generalized SBP operators with appropriate simultaneous approximation terms (SATs) [71]. However, due to the presence of E within the boundary term E T BE, the accurate and entropy stable imposition of nonlinear boundary conditions using generalized SBP-SATs remains an open problem [39,72,73]. We address this issue by constructing a hybridized SBP operator [40,41] (also referred to as a decoupled SBP operator in [40]) Note that Q h satisfies a block version of the standard summation by parts property by construction where B h is diagonal and does not involve the boundary evaluation matrix E.
It is less straightforward to show that Q h 1 = 0, which is also necessary for the proof of entropy stability [41]. Note that the full order matrix Q in (23) satisfies Q1 = 0. Then, since 1 is exactly representable in the test basis V t , V t P t 1 = 1, and This implies that Q h 1 simplifies to where we have used that Q t 1 = 0 and that (14) implies E1 = 1. Here, we abuse notation by letting 1 denote the vector of ones with appropriate dimension. Furthermore, since Q t satisfies a generalized SBP property, By replacing Q t in (18) with Q h , we can impose boundary conditions weakly through a numerical flux f * . We define V b as the matrix which evaluates the reduced basis at boundary points {−1, 1}, and define V h as the matrix which evaluates the reduced basis at both boundarypoints and hyper-reduced volume points We then have the following theorem: If the boundary flux f * is entropy stable in the sense of (25), then the hyper-reduced ROM is also entropy stable such that with equality holding for an entropy conservative flux.
Proof. The proof is identical in structure to the proof of Theorem 1 in [40]. We reproduce it here for completeness. Testing with v N and using the summation by parts property of Q h then yields Here, we have used that, by the consistency of the entropy conservative flux (3), the diagonal entries of F are F ii = f ( u i ). Thus, since B h is also diagonal, Expanding out the term v T Q h − Q T h • F 1 and using skew- where we have used that Q h 1 = 0 and the SBP property in the last two lines. Using that B is diagonal, straightforward manipulations imply that If the flux is entropy stable, then 6. Extension to higher dimensions

Periodic domains
For periodic domains, the extension to higher spatial dimensions is relatively straightforward. We provide a concrete construction of matrices and formulations in two dimensions. For finite volume methods on a K × K structured quadrilateral elements of size ∆x × ∆x, differentiation matrices along each coordinate direction can be constructed as Kronecker products of one-dimensional matrices. Let Q 1D denote the onedimensional periodic differentiation matrix defined in (5). The differentiation matrix along the ith coordinate direction Q i can be constructed as follows Given a basis matrix V and indices of hyper-reduced points I, we can introduce test basis matrices V i t such that the range of V i t is equal to the direct sum of the ranges of V and Note that, unlike the 1D case, the test bases differ along each coordinate direction. However, the dimension of the range of each test basis V i t is still at most 2N + 1. We can now define hyper-reduced projection matrices for the ith coordinate direction These matrices can be used to construct reduced nodal differentiation matrices Q i which are used to construct an entropy conservative reduced order model in d dimensions Here, f i S (u L , u R ) is the entropy conservative flux in the ith coordinate direction and M N is the mass matrix as defined in (18). We can show that (27) is entropy conservative by repeating steps in the proof of Theorem 3 for each coordinate direction.
The diffusive matrices K can be similarly extended to higher dimensions via Kronecker products An hyper-reduction of K analogous to the one described in Section 4.6 can be used to guarantee a discrete dissipation of entropy.

Non-periodic domains in higher dimensions
Non-periodic domains and boundary conditions can be treated by combining hybridized SBP operators and a weak imposition of boundary conditions as was done for the 1D hyper-reduced formulation (26). Without loss of generality, we assume a square domain with 4 boundary faces, each of which is discretized by K 1D intervals of size ∆x. Denote values of the outward normal vector on the domain boundary by n = [n 1 , . . . , n d ] T . We enforce non-periodic boundary conditions through numerical fluxes at boundary points. However, the number of boundary points scales with K 1D . To ensure that the cost of the reduced order model scales with the number of modes N rather than the dimension of the full order model, we approximate boundary terms by a hyper-reduced weighted combination of sampled boundary points.
Let V b denote the matrix whose columns correspond to values of the reduced basis at boundary points. Let I b denote a sub-sampled set of N b boundary points, and let g denote a nonlinear function. We then approximate the boundary inner product is a diagonal matrix whose entries consist of positive weights. We also define N b × N b operators B 1 , B 2 as the diagonal matrices whose entries consist of weighted component values of outward normals on the boundary where n i is a vector containing the values of the component n i at boundary points.
We can now construct hybridized SBP operators using hyper-reduced points on both the volume and boundary of the domain. Let V bt be the matrix which evaluates the ith test basis at all points on the boundary. Then, we can define E i as the matrix which extrapolates from hyper-reduced volume points to hyper-reduced boundary points through projection onto the ith test basis The hybridized SBP operator for differentiation along the ith coordinate is then These operators can be used to construct a hyper-reduced formulation in higher dimensions where M N and F i are as defined in (27), and f i (u), f * i are the ith components of the flux function and boundary numerical flux, respectively.
In contrast to the 1D case, additional steps are necessary to ensure that entropy stability is preserved under this boundary hyper-reduction. Recall that the proof of entropy stability for 1D non-periodic domains in Theorem 4 requires that Q h 1 = 0. The analogous proof of entropy stability in multiple dimensions requires that Q i h 1 = 0 for i = 1, . . . , d. However, this is not automatically satisfied under an arbitrary hyper-reduction of the boundary points. Expanding out Q i h 1 yields In general, − Q i t T 1 + E T B i 1 = 0, so Q i h 1 = 0 and the proof of entropy stability does not hold. To enforce that Q i h 1 = 0, we impose constraints on the boundary weights w b . Note that where we have used that V i t P i t 1 = 1. Thus, to ensure that Q i h 1 = 0, it is sufficient to guarantee that The dN constraints encoded in (29) can be interpreted as enforcing a discrete version of the fundamental theorem of calculus relating approximate integrals of reduced basis derivatives to boundary integrals of reduced basis values. We enforce these constraints on the boundary weights w b (29) directly into a hyperreduction approach based on the solution of a linear programming problem using the dual simplex method, and refer the reader to [59] for details. We note that it may also be possible to enforce these constraints, for example, by augmenting the non-negative least squares solve in Algorithm 1 with equality constraints, and will explore these directions in future work.

Numerical experiments
In this section, we study the behavior of entropy stable reduced order models for the compressible Euler equations. In d dimensions, these are given by: Here, ρ is density, u i is the ith component of velocity, and E is the total energy. The pressure p and specific internal energy ρe are given by There is a unique entropy which symmetrizes the viscous heat conduction term in the compressible Navier-Stokes equations [44]. This entropy S(u) is given by where s = log p ρ γ is the physical specific entropy, and the dimension d = 1, 2, 3. The entropy variables in d dimensions are given by while the conservation variables in terms of the entropy variables are given by where the quantities ρe and s in terms of the entropy variables are Let f i , f j denote two arbitrary values. We define the average and logarithmic average .
To ensure numerical stability when the denominator is close to zero, we evaluate the logarithmic average using the algorithm of [65]. Explicit expressions for entropy conservative fluxes are given by Chandrashekar [66]. In 1D, these fluxes are where the auxiliary quantities are In 1D, u 2 avg = 2 {{u}} 2 − {u 2 } . In two dimensions, the x and y fluxes are given by The viscous terms are treated as described in Section 2.1, and utilize evaluations of the Jacobian ∂u ∂v . Explicit expressions for Jacobians in the compressible Euler and Navier-Stokes equations are given in [44]. The solution is evolved using an explicit 5-stage 4th order Runge-Kutta scheme.

1D Euler equations
We begin by examining the 1D Euler equations with reflective wall boundary conditions. We utilize a full order finite volume model with 2500 cells on [−1, 1] with a CFL of .75, and set the artificial viscosity parameter to = 2 × 10 −4 . It was shown in [77,38] that wall boundary conditions can be imposed using a "mirror state". We use a boundary numerical flux f * = f S (u + , u) augmented with local Lax-Friedrichs penalization, where the exterior state u + is defined by The initial condition is set to be The solution exhibits a viscous shock some time after T = .25.
We first examine the impact of enriching solution snapshots with snapshots of the entropy variables. We run the full order model until final time T = .7 and store the solution at 2801 time-steps. We compute a reduced basis using SVD by subsampling every 10 snapshots. Figure 2a shows the decay of the singular values with and without entropy variable enrichment. In both cases, the decay is similar, with entropy variable enrichment resulting in a slightly slower decay for singular values smaller than 10 −8 . The difference becomes more pronounced for a coarser subsampling of the snapshots. We note that the enrichment does affect the form of the resulting singular vectors (Figure 2b). We also compute projection errors for each solution snapshot using 25 modes in Figure 2c, along with differences between the projected solutions with and without entropy variable enrichment. The projection errors are nearly indistinguishable.
We next examine solutions produced by the reduced order models. It is known that reduced models can sometimes utilize a larger CFL compared the full order model [74,75]; however, to ensure that temporal errors are small, all solutions are computed using the same CFL used to generate the solution snapshots. Figure 3 shows density computed using 25, 75, and 125 modes and the hyper-reduced treatment of viscosity in (21). For all resolutions, the shock is under-resolved and the solution possesses Gibbs-type oscillations. However, despite this under-resolution, the ROM remains stable and does not blow up. Moreover, the reduced order solution converges uniformly to the full order solution as the number of modes increases. Table 1 shows the number of hyper-reduced points computed for different numbers of modes. We report on three different sets of points: empirical cubature points (used to approximate the mass matrix ∆xV T V ), stabilizing points (added to control the condition number of the test mass matrix; see Section 4.3), and viscous points (used to approximate viscous terms; see Section 4.6). All point sets are computed using Algorithm 1 with dimensionality reduction as described in Section 4.2. In all cases, we observe that the number of empirical cubature and viscous points grows roughly as 2N (similar observations were made in [33]), while the number of stabilizing points hovers around 30 for N > 25.
Next, we verify that (in the absence of viscosity), the ROM is semi-discretely entropy conservative. The proof of this property uses that the entropy contribution of the convective term v T    (26) is discretely entropy conservative in the absence of viscosity and boundary flux penalization terms. We also run the ROM with the viscosity parameter = 0. For 125 modes, the solution is significantly more oscillatory than the case of = 2 × 10 −4 ; however, the reduced solution does not blow up.
Next, we examine the difference in the hyper-reduced treatments of entropy dissipation (19), (21), and (22) described in Section 4.6 (recall that (22) is not provably entropy dissipative). If the discrete entropy dissipation v T N d(u N ) is positive, then the simulation is entropy stable. Figure 5 shows the computed entropy dissipation over the time interval [0, .75]. All hyper-reduced treatments produce similar results, with the entropy dissipation produced by the naive approximation (22) differing most significantly. The naive approximation differs from the provably entropy stable approximations (19) and (21)  We now compare the evolution of both the L 2 error and discrete entropy, and examine the effect of hyper-reduction on the discrete solution in Figure 6. Both error and entropy are computed on the original grid of the full order model (e.g., without hyper-reduction), and the error is computed as the discrete relative L 2 error over all the conservative variables. We observe that the evolution of the average discrete entropy for the reduced model rapidly approaches the average discrete entropy for the full order model as the number of modes increases from 25 to 75. For 125 modes, the average entropies for the ROM and FOM  are indistinguishable and are not shown. The errors behave similarly, decreasing as the number of modes increases. For reference, we also include errors computed without hyper-reduction using (11). We observe that the hyper-reduced errors are virtually identical to errors without hyper-reduction in all cases.

2D Euler equations
We now study the Euler equations in 2D. We begin by considering periodic boundary conditions on the domain [−1, 1] 2 and a smoothed version of the Kelvin-Helmholtz instability. The initial condition is adapted from [76,30]    In the following experiments, we set α = .1 and σ = .1.
We utilize a full order model with 200×200 grid (40000 cells), which is run until time T = 3. We generate 301 solution snapshots with which to compute POD modes. The decay of the singular values of the solution snapshots with and without enrichment by the entropy variables is shown in Figure 7. We observe that, compared to the 1D case, the singular values decay more slowly when enriched with snapshots of the entropy variables. However, as in 1D, the differences between the singular values with and withoutentropy variable enrichment are significantly larger for lower energy modes as compared to higher energy modes. The reduced order model utilizes 75 POD modes to achieves a 1.02% relative L 2 error. We note that the difference between the 75th singular value of the snapshot matrix with entropy variable enrichment is only about 6% larger than the 75th singular value without entropy variable enrichment. This suggests that entropy variable enrichment does not significantly impact the accuracy of the 75-mode POD approximation.
Since little difference is observed between the hyper-reducted treatments of viscosity, the naive treatment (22) is used. No additional stabilizing points were necessary for this setting, as the condition numbers of the hyper-reduced x and y test mass matrices were both O(1). Figure 8 compares results for the full and reduced order models, and shows the hyper-reduced points selected by the greedy empirical cubature algorithm.
We next consider a 2D domain with wall boundary conditions. Entropy stable reflective wall boundary  conditions are again imposed using "mirror states". Let n denote the outward normal at a wall, and let u n = un x + vn y and u n ⊥ = un y − vn x denote the normal and tangential components of the velocity at the wall. We use a boundary numerical flux f * = f S (u + , u) augmented with local Lax-Friedrichs penalization, where the exterior state u + is defined by ρ + = ρ, u + n = −u n , u + n ⊥ = u n ⊥ , p + = p.
We set a Gaussian pulse initial condition ρ = 1 + e −50(x 2 +(y+1/2) 2 ) , The full order model is computed on a 150 × 150 grid with viscosity coefficient = 1e − 3. The POD basis is computed from 75 solution snapshots with a CFL of .5. Figure 9 shows the decay of the singular values with and without entropy variable enrichment. As before, we observe that the singular values decay more slowly under entropy variable enrichment, but that the difference in the singular values with and without enrichment is more evident for low energy modes. The 25th singular value with entropy variable enrichment is less than 5% larger than the 25th singular value without entropy variable enrichment, suggesting that entropy variable enrichment does not significantly affect the accuracy of the 25-mode POD approximation. to approximately 5.71. The solution does not form shock discontinuities, and is chosen instead to test the imposition of boundary conditions. We note that, due to the nature of the linear programming software used to compute the boundary weights, the constraints (29) on the boundary weights w b are imposed to a tolerance of 5e − 8. Despite this inexact enforcement of constraints, the computed entropy RHS remained small, oscillating around 10 −11 in the absence of entropy dissipative terms.

Conclusion
We have presented a methodology for constructing projection-based reduced order models for nonlinear conservation laws by combining a reduced basis, a modified Galerkin projection, and tailored hyper-reduction techniques. The main novelty of these new reduced models is the approximation of the nonlinear convective term, which combines a "flux differencing" approach with an appropriate hyper-reduced approximation of the differentiation matrix.
Future work will aim to address computational costs associated with entropy stable reduced models. Unlike standard hyper-reduction techniques, the number of nonlinear evaluations necessary scales with O(N 2 s ) rather than O(N s ), where N s is the number of hyper-reduced sampling points. This discrepancy is due to the fact that the hyper-reduction presented here approximates a nonlinear matrix, rather than a nonlinear vector. These costs can be reduced by combining domain decomposition [78] and a discontinuous Galerkintype discretization [40] to produce multi-domain reduced order models. Suppose there are k = 1, . . . , K subdomains. Then, a multi-domain reduced order model is expected to decrease costs if the sum of the squares of the number of subdomain sampling points k (N k s ) 2 is significantly smaller than the global number of sampling points N 2 s . Future work will also investigate implicit time-stepping for entropy stable reduced models, as additional costs associated with the spatial formulation may be offset by the reduction in cost for solving smaller matrix systems during the linearization process.

Acknowledgments
Jesse Chan gratefully acknowledges support from the National Science Foundation under awards DMS-1719818 and DMS-1712639, and thanks Matthias Heinkenschloss, Masayuki Yano, Matthew Zahr, and Irina Tezuar for informative discussions.