Machine-precision energy conservative reduced models for Lagrangian hydrodynamics by quadrature methods

TL;DR

This paper introduces an energy conservative reduced modeling approach for Lagrangian hydrodynamics that ensures near machine-precision energy conservation while maintaining high accuracy, using a quadrature-based model reduction framework.

Contribution

The paper develops a strongly energy conservative variant of empirical quadrature for reduced models of Euler equations, ensuring exact energy conservation during model reduction.

Findings

Conserves total energy to near machine precision.

Maintains accuracy comparable to basic EQP methods.

Validated on four benchmark problems.

Abstract

We present an energy conservative, quadrature based model reduction framework for the compressible Euler equations of Lagrangian hydrodynamics. Building on a finite element discretization of the governing equations, we develop reduced models using data based reduced basis functions and the empirical quadrature procedure (EQP). We introduce a strongly energy conservative variant of EQP that enforces exact energy conservation in the reduction process. Numerical experiments for four benchmark problems -- Sedov blast, Gresho vortex, triple point and Taylor-Green vortex -- demonstrate that the numerical implementation of our proposed method conserves total energy to near machine precision, while maintaining accuracy comparable to the basic EQP formulation.

Tables

Table 1. Table 1 : Spatial discretization parameters for each simulated problem. Parameter m m controls the mesh size h = 2 − m ​ h 0 h=2^{-m}h_{0} , where h 0 h_{0} denotes the coarsest mesh size; k k is the polynomial order of the basis used for the kinematic finite element space 𝒱 \mathcal{V} , with order k − 1 k-1 for the thermodynamic space ℰ \mathcal{E} ; N v N_{v} and N e N_{e} denote the global number of basis elements for the finite element spaces 𝒱 \mathcal{V} and ℰ \mathcal{E} respectively.

Run	R1	R2	R3	R4
Problem	Sedov	Gresho	Triple point	Taylor-Green
$m$	2	4	2	2
$k$	2	2	2	2
$N_v$	14739	18818	38475	14739
$N_e$	4096	9216	10752	4096

Table 2. Table 2 : Parameters and results of the reduced simulations. In the table t f t_{f} denotes the final time of the simulation; e σ e_{\sigma} the POD energy ratio ( 3.2 ); N s N_{s} the number of samples per window; N w N_{w} the total number of time windows. The symbols ϵ i \epsilon_{i} denote the absolute value of the relative error of the reduced solution for position ( x x ), velocity ( v v ) and energy ( e e ); Δ ​ E \Delta E the absolute value of the relative energy difference of each simulation (start to finish); S S the speedup of the reduced simulation relative to the full one. The ratio J ^ / J \hat{J}/J denotes the rounded ratio of quadrature points of the reduced rule over that of the full rule, averaged over all windows; for basic EQP there are two separate ratios for velocity ( v v ) and energy ( e e ); for conservative EQP there is only one ratio for the combined quadrature rule. For all results subscript B B refers to basic EQP; C C to energy conservative EQP.

Run	R1	R2	R3	R4
Problem	Sedov	Gresho	Triple point	Taylor-Green
$t_f$	0.3	0.4	0.8	0.4
$e_{\sigma }$	0.9999	0.9999	0.9999	0.9999
$N_s$	10	10	10	10
$N_w$	91	132	39	418
${ϵ}_{x,B}$	2.30e-03	7.03e-07	2.19e-07	9.85e-09
${ϵ}_{v,B}$	1.37e-02	1.19e-05	1.46e-05	1.41e-05
${ϵ}_{e,B}$	2.04e-04	2.15e-06	4.74e-06	2.27e-07
$\Delta E_B$	4.23e-04	1.72e-08	1.15e-06	3.03e-11
${\hat{J}}_v/J_v$	$103/32768$	$47/16384$	$46/86016$	$33/32768$
${\hat{J}}_e/J_e$	$61/32768$	$67/16384$	$45/86016$	$48/32768$
$S_B$	2.21	2.74	1.45	2.24
${ϵ}_{x,C}$	4.98e-05	3.04e-03	4.32e-06	2.38e-05
${ϵ}_{v,C}$	5.93e-03	5.61e-02	9.81e-05	1.13e-03
${ϵ}_{e,C}$	2.34e-03	8.87e-03	5.43e-05	1.36e-05
$\Delta E_C$	2.28e-13	9.28e-15	4.43e-14	1.38e-14
$\hat{J}/J$	$179/32768$	$133/16384$	$112/86016$	$104/32768$
$S_C$	1.68	2.60	1.38	1.69

Equations

$$\dot{\hat{u}}={\Phi }^{\top }f(u)={\Phi }^{\top }f(\Phi \hat{u}).$$

$$\text{mass conservation}:$$

$$\text{momentum conservation}:$$

$$\text{energy conservation}:$$

$$\text{equation of motion}:$$

$$p=(\gamma -1)\rho e.$$

$$\rho (t,x_0)=\frac{\rho (0,x_0)}{|J(t,x_0)|}$$

$$v(t,x_0)\approx v_f(t,x_0)$$

$$e(t,x_0)\approx e_f(t,x_0)$$

$$x(t,x_0)\approx x_f(t,x_0)$$

$$\text{momentum conservation}:$$

$$\text{energy conservation}:$$

$$\text{equation of motion}:$$

$$M_{v,ij}={\int }_{\Omega (t)}\rho {\theta }_i^v{\theta }_j^{v}dx{M}_{e,ij}={\int }_{\Omega (t)}\rho {\theta }_i^e{\theta }_j^{e}dx$$

$$F_{v,i}(\mathbf{w})={\int }_{\Omega (t)}\sigma (w_f):\nabla {\theta }_i^{v}dx{F}_{e,i}(\mathbf{w},{\mathbf{v}}^{\prime })={\int }_{\Omega (t)}(\sigma (w_f):\nabla v^{\prime }){\theta }_i^{e}dx.$$

$$F_{ij}(\mathbf{w})={\int }_{\Omega (t)}(\sigma (w_f):\nabla {\theta }_i^v){\theta }_j^{e}dx$$

$${\mathbf{F}}_v(\mathbf{w})=\mathbf{F}(\mathbf{w})1_{\mathcal{E}}{\mathbf{F}}_e(\mathbf{w},{\mathbf{v}}^{\prime })=\mathbf{F}(\mathbf{w}{)}^{\top }{\mathbf{v}}^{\prime }$$

$${\int }_{\Omega (t)}a(t,x(t,x_0))dx(t,x_0)={\int }_{{\Omega }_0}a(t,x_0)|J(t,x_0)|d{x}_0$$

$${\mathbf{v}}_{k+1/2}$$

$${\mathbf{e}}_{k+1/2}$$

$${\mathbf{x}}_{k+1/2}$$

$$v(t,x_0)\approx \tilde{v}(t,x_0)$$

$$e(t,x_0)\approx \tilde{e}(t,x_0)$$

$$x(t,x_0)\approx \tilde{x}(t,x_0)$$

$$\frac{{\sum }_{i=0}^{n_k-1}{\sigma }_i}{{\sum }_{i=0}^{N_{snap}-1}{\sigma }_i}\ge e_{\sigma }$$

$$\mathbf{v}(t)\approx \tilde{\mathbf{v}}(t)$$

$$\mathbf{e}(t)\approx \tilde{\mathbf{e}}(t)$$

$$\mathbf{x}(t)\approx \tilde{\mathbf{x}}(t)$$

$${\mathbf{M}}_v{\mathbf{\Phi }}_v\frac{d\hat{\mathbf{v}}}{dt}$$

$${\mathbf{M}}_e{\mathbf{\Phi }}_e\frac{d\hat{\mathbf{e}}}{dt}$$

$${\mathbf{\Phi }}_x\frac{d\hat{\mathbf{x}}}{dt}$$

$$\frac{d\hat{\mathbf{v}}}{dt}$$

$$\frac{d\hat{\mathbf{e}}}{dt}$$

$$\frac{d\hat{\mathbf{x}}}{dt}$$

$$F_{v,i}^{\psi }(\tilde{\mathbf{w}})={\int }_{\Omega (t)}\sigma (\tilde{w}):\nabla {\psi }_i^{v}dx{F}_{e,i}^{\psi }(\tilde{\mathbf{w}},{\tilde{\mathbf{v}}}^{\prime })={\int }_{\Omega (t)}(\sigma (\tilde{w}):\nabla {\tilde{v}}^{\prime }){\psi }_i^{e}dx$$

$${\hat{\mathbf{M}}}_v\frac{d\hat{\mathbf{v}}}{dt}$$

$${\hat{\mathbf{M}}}_e\frac{d\hat{\mathbf{e}}}{dt}$$

$$\frac{d\hat{\mathbf{x}}}{dt}$$

$${\hat{\mathbf{M}}}_v={\mathbf{\Phi }}_v^{\top }{\mathbf{M}}_v{\mathbf{\Phi }}_v\in {\mathbb{R}}^{n_v\times n_v}{\hat{\mathbf{M}}}_e={\mathbf{\Phi }}_e^{\top }{\mathbf{M}}_e{\mathbf{\Phi }}_e\in {\mathbb{R}}^{n_e\times n_e}.$$

$${\hat{M}}_{v,ij}={\int }_{\Omega (t)}\rho {\varphi }_i^v{\varphi }_j^{v}dx{\hat{M}}_{e,ij}={\int }_{\Omega (t)}\rho {\varphi }_i^e{\varphi }_j^{e}dx$$

Full text

\affiliation

[1]organization=Department of Mathematics, Dartmouth College, city=Hanover, state=NH 03755, country=USA \affiliation[2]organization=Center for Applied Scientific Computing, Lawrence Livermore National Laboratory, city=Livermore, state=CA 94550, country=USA

1 Introduction

High fidelity simulations of multiscale, multiphysics systems often involve a large number of degrees of freedom evolving over widely separated temporal and spatial scales, making them prohibitively expensive for tasks such as design, control or uncertainty quantification [1, 2]. Model reduction methods—and the related closure or parametrization methods—aim to construct low dimensional surrogate models that retain the essential dynamics of the original system, while reducing the computational cost associated with its numerical simulation [3, 4, 5]. Reduced models can be derived using both physical modeling and data based approaches, depending on the complexity of the underlying dynamics and the availability of high quality data, among other determining factors. In either case, reduction strategies employ a wide variety of mathematical tools, including both deterministic and stochastic methods.

In this work we derive reduced models that govern the evolution of the linear projection of the original state variables onto a subspace spanned by a set of reduced basis functions. Focusing on Lagrangian hydrodynamics, our approach combines physical knowledge about the underlying system in the form of its governing equations, with basis functions computed in a data based fashion.

In the remainder of Section 1 we discuss projection based reduction methods and modeling approaches for Lagrangian hydrodynamics, followed by a summary of the main results of this work. In Section 2 we present the governing equations and the discretization used to obtain the full model. The reduced model is derived in Section 3, followed by a presentation of the EQP hyper-reduction method and its implementation. The energy conservative EQP method is developed in Section 4. Section 5 presents our numerical results for four different problem cases using both EQP methods, followed by a conclusion in Section 6. Finally, A includes additional details on the numerical implementation of the EQP methods.

1.1 Projection based reduction

Projection based reduction methods derive reduced models by replacing the original state space of the dynamics with a subspace defined as the linear span of a set of reduced basis functions [1, 6]. The reduced state variables are then defined by the linear projection of the original state variables to the constructed subspace. The reduced basis functions that define the state space of the reduced dynamics can be computed by various data based methods, which usually employ the singular value decomposition (SVD) or eigenvalue decomposition of an appropriate data matrix.

The most common method is arguably the proper orthogonal decomposition (POD), which uses the SVD of a snapshot data matrix to construct an orthonormal basis for its range space [7, 8, 9]. Given a collection of samples of the state variables, POD yields the best linear fixed-rank approximation to the data in $\mathbb{R}^{N}$ with the euclidean inner product, where $N\in\mathbb{N}$ denotes the state space dimension. An alternative is the dynamic mode decomposition (DMD) [10], which can be used to approximate the eigenspaces of the system’s Koopman operator restricted to the finite dimensional subspace spanned by the sampled data [11, 12, 13].

The POD and DMD methods can be viewed as specific instances of kernel methods, where one employs a symmetric kernel function $k$ to compute the pairwise distances matrix $K_{ij}=k(u_{i},u_{j})$ from state samples $u_{i}$, $u_{j}\in\mathbb{R}^{N}$, and uses the eigenvectors of the kernel matrix $K$ as the reduced basis functions. Compared to POD, general kernel methods replace the correlation kernel based on the euclidean inner product with the kernel function $k$, thereby embedding the state samples into a feature space before computing their pairwise distances [14]. Kernel methods can also be used to approximate the eigenspaces of the restricted Koopman operator by using kernels that act on delay embedded state samples [15].

The effectiveness of linear subspace projection methods is generally limited to problems where the solution manifold can be accurately approximated by a low dimensional subspace. This assumption is not satisfied in problems with a slow decay in the Kolmogorov $n$-width, such as advection dominated problems that feature sharp gradients, moving shock fronts and turbulence [16, 17]. One approach toward addressing this issue has been to construct local reduced models that approximate the original dynamics for a specific subset of the considered parameter-time domain [18, 19]. A reduction approach based on time windows was introduced in [20, 21] to build temporally local reduced models that are small but accurate over a short time period for advection dominated problems. This approach has since been further extended and applied to problems in various domains [22, 23, 24].

The methods outlined above form a subspace of the original state space and define the reduced state variables by linear projection. An alternative approach is to construct a nonlinear projection of the original state variables onto an appropriate subspace, often called a latent space. The projection map is often represented by a neural network, with the latent space defined by the architecture of the employed neural network and its training procedure [25, 26, 27, 28].

In addition to accuracy and reduced complexity, it is often desirable for reduced models to preserve aspects of the qualitative structure of the original dynamics, such as conserved quantities or symplectic structure. Several reduction schemes have been proposed that incorporate structure preservation constraints directly into the formulation of the reduced model, ensuring that the relevant structure is not destroyed in the reduction process [29, 30, 31, 32, 33, 34, 35, 6, 36, 37, 38, 39, 40].

1.2 Hyper-reduction methods

We consider a dynamical system with state space $\mathbb{R}^{N}$, $N\in\mathbb{N}$, governed by the differential equation $\dot{u}=f(u)$ with state vector $u\in\mathbb{R}^{N}$ and nonlinear function $f\colon\mathbb{R}^{N}\to\mathbb{R}^{N}$. We assume access to an orthonormal reduced basis matrix $\Phi\in\mathbb{R}^{N\times n}$ with $n\ll N$, where each column represents a reduced basis function. Using the representation $u\approx\Phi\hat{u}$, $\hat{u}\in\mathbb{R}^{n}$, we derive the governing equation for the reduced state vector $\hat{u}$ by projection

[TABLE]

Although this is an equation for the $n\ll N$ reduced state variables, the evaluation of the nonlinear function $f$ at every timestep depends on $\Phi\hat{u}\in\mathbb{R}^{N}$ and is expected to dominate the model’s computational cost. For this reason, various hyper-reduction methods have been developed to approximate the evaluation of $f$ in different ways.

Hyper-reduction methods can be divided into interpolation and quadrature based methods. Interpolation methods employ a set of collocation points to approximate the evaluation of nonlinear function $f$, an idea that can be traced back to [41] in the context of image reconstruction. Examples of this class of methods include the empirical interpolation method [42] and its discrete variant [43], and employ various sampling strategies to select the collocation points based on different optimality criteria [44, 45, 46, 47]. Quadrature methods offer an alternative strategy for efficiently approximating integrals of nonlinear functions arising in finite element models. They derive reduced quadrature rules that approximate the spatial integration of functions over the computational domain by requiring their evaluation at a limited number of quadrature points. The origin of this approach can be traced back to [48] in the context of computer graphics applications. Examples in this class of methods include the energy conserving sampling and weighting method [29, 31], the empirical cubature method [49], the EQP method [50, 51, 52] and more [53, 54].

1.3 Lagrangian hydrodynamics

We consider the compressible Euler equations of fluid dynamics in a three dimensional spatial domain, which can be used to model inviscid, high speed flows and shock wave propagation [55, 56, 57]. The numerical solution of the equations can generally be achieved by two classes of methods: Eulerian and Lagrangian methods. In Eulerian methods the unknown variables are treated as fields over a fixed computational mesh, whereas in Lagrangian methods the mesh moves according to the local fluid velocity [58].

Lagrangian methods have traditionally been formulated using approaches such as the staggered grid hydrodynamics and the cell centered hydrodynamics methods, which usually employ finite difference or finite volume spatial discretizations. The accuracy of these methods relies heavily on the quality of the underlying mesh, which moves and deforms over time and can lead to phenomena such as mesh tangling. One way to address this issue has been the use of the arbitrary Lagrangian-Eulerian method, which involves a Lagrangian phase followed by a remesh and remap phase. In the Lagrangian phase, the computational mesh evolves according to the fluid velocity until its quality deteriorates in a prescribed sense. Once that happens, the mesh is adjusted accordingly and the solution fields are remapped onto the adjusted mesh, allowing the Lagrangian simulation to resume. We refer the reader to [57] for details on the methods outlined above and their historical development.

In this work we employ the Lagrangian method developed in [57], which uses a high order finite element discretization of the Euler equations in two and three dimensional cartesian grids. Among other features, the employed method offers strong conservation of mass and discrete total energy. We refer the reader to [57, 59] for the detailed development of the method and its application to various benchmark problems.

1.4 Main results

We resume the line of research initiated in [22, 23], where the authors derived projection based reduced models for Lagrangian hydrodynamics problems using interpolation hyper-reduction methods. In the present work we use quadrature hyper-reduction methods, which are naturally well suited to the finite element discretization used in the full model. The main features of this work can be summarized as follows.

Projection based reduction. Starting from a finite element (FEM) discretization of the compressible Euler equations, we develop projection based reduction models using reduced basis functions computed from simulation data. We derive two reduced problem formulations, demonstrating that the reduction process corresponds to replacing the original FEM trial and test functions with the reduced basis functions. We use the EQP method to reduce the cost of the spatial integration of nonlinear force functions that arises in the FEM discretization process. 2. 2.

Conservative quadrature. We develop an energy conservative variant of the basic EQP method that enforces exact conservation of the discrete total energy in the strong sense in the hyper-reduced model. This is achieved by deriving a reduced quadrature rule that preserves the energy conservation property of the full model. 3. 3.

Applications. We implement the hyper-reduced models derived using the basic and energy conservative EQP methods in the Laghos111Code repository: https://github.com/CEED/Laghos Lagrangian hydrodynamics simulation code using the libROM222Code repository: https://github.com/LLNL/libROM library for model reduction. We apply the methods to four standard benchmark problems, reporting their accuracy and obtained computational speedup. Our numerical results show that the numerical implementation of the energy conservative EQP method leads to conservation of total energy to near machine precision for all simulated problem cases.

2 Full model

We consider the system of Euler equations of fluid dynamics in a Lagrangian reference frame with no external body forces

[TABLE]

where $d/dt$ denotes the material derivative, $\rho$ the density of the fluid, $x$ and $v$ the position and velocity of the particles in a deformable medium $\Omega(t)$ in the Eulerian coordinates, $\sigma$ the deformation stress tensor, and $e$ the internal energy per unit mass. These physical quantities are treated as functions of time $t\geq 0$ and the initial configuration of particles $x_{0}\in\Omega_{0}=\Omega(0)$ [57, 56, 55]. The evolution equation for the position field $x(t,x_{0})$ governs the Lagrangian advection of the fluid particles; in the spatially discretized system, this leads to a computational mesh that deforms according to the local fluid velocity.

We use the isotropic stress tensor $\sigma=-pI+\sigma_{a}$, where $p$ denotes the thermodynamic pressure and $\sigma_{a}$ the artificial viscosity stress. The thermodynamic pressure is governed by a constitutive equation and can be expressed as a function of the density and internal energy. In the present work we focus on the case of polytropic ideal gases with an adiabatic index $\gamma>1$, which yields the constitutive equation

[TABLE]

The system is prescribed with an initial condition and the wall boundary condition $v\cdot n=0$, where $n$ denotes the outward normal unit vector on the domain boundary.

2.1 Spatial discretization

We follow [57] to derive a semidiscrete Lagrangian variational formulation based on a computational mesh that deforms according to the fluid velocity. By the Reynolds transport theorem, the mass conservation in (2.1) is related to the Jacobian $\nabla_{x_{0}}x(t,x_{0})$ of the Lagrangian transformation $x(t,x_{0})$ by

[TABLE]

where $J(t,x_{0})=\nabla_{x_{0}}x(t,x_{0})$ is the deformation gradient and $|J|$ its determinant [57]. The above equation can be used to determine the density field for all $t\geq 0$. As a result, from now on we focus on the differential equations for the state variable $w=(v,e,x)$.

For the spatial discretization of (2.1) we employ a finite element method (FEM) using a kinematic space $\mathcal{V}\subset[H^{1}(\Omega_{0})]^{d}$ with spatial dimension $d\in\mathbb{N}$ and basis $\{\theta_{v,i}\}_{i=0}^{N_{v}-1}$ for approximating the position and velocity fields, and a thermodynamic space $\mathcal{E}\subset L_{2}(\Omega_{0})$ with basis $\{\theta_{e,i}\}_{i=0}^{N_{e}-1}$ for approximating the energy field, where $N_{v}$, $N_{e}\in\mathbb{N}$ denote the number of the global degrees of freedom in the corresponding spaces. In what follows, the FEM basis functions are used as both trial and test functions.

We approximate the full solution fields by

[TABLE]

We use $\bm{v}$, $\bm{x}\in\mathbb{R}^{N_{v}}$ and $\bm{e}\in\mathbb{R}^{N_{e}}$ to denote the FEM coefficient vectors of the corresponding solution fields and $\bm{w}=(\bm{v},\bm{e},\bm{x})^{\top}\in\mathbb{R}^{N}$, $N=2N_{v}+N_{e}$ to denote the full state vector. The semidiscrete Lagrangian conservation laws can be derived by testing equations (2.1) against the Lagrangian extensions of the FEM basis functions on $\Omega(t)$

[TABLE]

with mass matrices $\bm{M}_{v}\in\mathbb{R}^{N_{v}\times N_{v}}$ and $\bm{M}_{e}\in\mathbb{R}^{N_{e}\times N_{e}}$

[TABLE]

and vector valued nonlinear force functions $\bm{F}_{v}:\mathbb{R}^{N}\to\mathbb{R}^{N_{v}}$ and $\bm{F}_{e}:\mathbb{R}^{N}\times\mathbb{R}^{N_{v}}\to\mathbb{R}^{N_{e}}$

[TABLE]

Introducing the matrix valued force function $\bm{F}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N_{v}\times N_{e}}$ with entries

[TABLE]

the nonlinear force functions can be compactly represented as

[TABLE]

with $\bm{1}_{\mathcal{E}}\in\mathbb{R}^{N_{e}}$ the coefficient vector of the unity function in $\mathcal{E}$.

In the above spatial integrations we write the integrals as taking place over the control volume $\Omega(t)$ instead of $\Omega_{0}=\Omega(0)$. For any scalar valued field $a(t,x_{0})$, $x_{0}\in\Omega_{0}$, this is equivalent to using its Lagrangian extension $a(t,x(t,x_{0}))$, $x(t,x_{0})\in\Omega(t)$, to perform the spatial integration; without misunderstanding, we use the same symbol to denote the two field descriptions. The two approaches are connected by performing the change of variables $x(t,x_{0})\mapsto x_{0}$

[TABLE]

where $|J(t,x_{0})|$ is the determinant of the deformation gradient defined in the previous section. To simplify the notation we use integrations over $\Omega(t)$ throughout this work, implicitly understanding that the above change of variables is carried out to perform the integrations over $\Omega_{0}$.

2.2 Temporal discretization

The time domain is discretized as $\{t_{k}\}_{k=0}^{K-1}$ with $t_{0}=0$ and $t_{K-1}=t_{f}$. All quantities at time $t_{k}$ are denoted by subscript $k$. We use the two-stage average Runge-Kutta (RK2A) scheme derived in [57], which has been shown to conserve the discrete total energy of the system. Applying the RK2A scheme leads to the discrete system of equations

[TABLE]

with $\bar{\bm{v}}_{k+1/2}=(\bm{v}_{k}+\bm{v}_{k+1})/2$. To ensure stability of the explicit timestepping scheme we also employ the timestep control algorithm described in [57].

3 Reduced model

In this section we derive the projection based reduced model for the semidiscrete Lagrangian conservation laws (2.3). We introduce the reduced bases for velocity $\{\phi_{v,i}\}_{i=0}^{n_{v}-1}\subset\mathcal{V}$, energy $\{\phi_{e,i}\}_{i=0}^{n_{e}-1}\subset\mathcal{E}$ and position $\{\phi_{x,i}\}_{i=0}^{n_{x}-1}\subset\mathcal{V}$ with reduced sizes $n_{v}$, $n_{x}\ll N_{v}$ and $n_{e}\ll N_{e}$. To approximate each solution field we restrict our solution space to the linear span of each basis

[TABLE]

with offset fields $v_{os}$, $x_{os}\in\mathcal{V}$ and $e_{os}\in\mathcal{E}$.

The FEM coefficient vectors of the reduced basis functions are represented by the orthonormal basis matrices $\bm{\Phi}_{v}\in\mathbb{R}^{N_{v}\times n_{v}}$, $\bm{\Phi}_{e}\in\mathbb{R}^{N_{e}\times n_{e}}$ and $\bm{\Phi}_{x}\in\mathbb{R}^{N_{v}\times n_{x}}$, where each column holds the FEM coefficients of one basis function. We compute the reduced basis matrices from simulation data of the full model using the POD method [9]. To decide how many POD basis functions to use to form each basis matrix we employ the POD energy ratio

[TABLE]

where $\sigma_{i}\geq 0$ denotes the $i$th largest singular value of the employed data matrix consisting of $N_{snap}$ solution snapshots, $e_{\sigma}\geq 0$ a chosen threshold and $n_{k}$ the size of the corresponding reduced state vector (velocity $v$, energy $e$ or position $x$). Using the introduced reduced basis matrices we rewrite (3.1) in terms of FEM coefficient vectors

[TABLE]

where $\hat{\bm{v}}\in\mathbb{R}^{n_{v}}$, $\hat{\bm{e}}\in\mathbb{R}^{n_{e}}$ and $\hat{\bm{x}}\in\mathbb{R}^{n_{x}}$ denote the reduced solution state vectors. We also use $\tilde{\bm{w}}=(\tilde{\bm{v}},\tilde{\bm{e}},\tilde{\bm{x}})^{\top}\in\mathbb{R}^{N}$ to denote the FEM coefficient vector of $\tilde{w}=(\tilde{v},\tilde{e},\tilde{x})$.

Making the above substitutions in (2.3) we reach the overdetermined semidiscrete system

[TABLE]

with the appropriate initial condition at $t=0$. The above corresponds to using the reduced basis functions as the trial functions and the FEM basis functions as the test functions. Next we present the derivation of two alternative reduced systems of equations for the reduced state vector $\hat{\bm{w}}=(\hat{\bm{v}},\hat{\bm{e}},\hat{\bm{x}})\in\mathbb{R}^{n}$, $n=n_{v}+n_{e}+n_{x}$. The reason for introducing two alternative systems is that both systems will be employed by the hyper-reduction methods introduced in the following sections.

Starting from (3.4), the first reduced system is derived by inverting the mass matrices and transpose multiplying by the reduced basis matrices

[TABLE]

with $\bm{\Psi}_{v}=\bm{M}_{v}^{-1}\bm{\Phi}_{v}\in\mathbb{R}^{N_{v}\times n_{v}}$ and $\bm{\Psi}_{e}=\bm{M}_{e}^{-1}\bm{\Phi}_{e}\in\mathbb{R}^{N_{e}\times n_{e}}$. This corresponds to using the functions represented by the columns of $\bm{\Psi}_{v}$ and $\bm{\Psi}_{e}$ as the test functions. More precisely, the components of force vectors $\bm{F}_{v}^{\psi}\in\mathbb{R}^{n_{v}}$ and $\bm{F}_{e}^{\psi}\in\mathbb{R}^{n_{e}}$ are given by the integrals

[TABLE]

where functions $\psi^{v}_{i}\in\mathcal{V}$ and $\psi^{e}_{i}\in\mathcal{E}$ have the $i$th column of matrices $\bm{\Psi}_{v}$ and $\bm{\Psi}_{e}$ respectively as their FEM coefficient vectors.

The second reduced system is derived by directly transpose multiplying by the reduced basis matrices in (3.4)

[TABLE]

with reduced mass matrices

[TABLE]

The above corresponds to using the reduced basis functions as both trial and test functions. More specifically, the reduced mass matrices are given by the integrals

[TABLE]

and the force vectors $\bm{F}_{v}^{\phi}\in\mathbb{R}^{n_{v}}$ and $\bm{F}_{e}^{\phi}\in\mathbb{R}^{n_{e}}$ by

[TABLE]

Although the size of the reduced systems (3.5) and (3.7) has been successfully reduced to $n\ll N$, the nonlinearity of the force functions in the right hand side of the equations means that they still depend on the full state vector $\tilde{\bm{w}}\in\mathbb{R}^{N}$. As a result, evaluating the force functions requires lifting the reduced state $\hat{\bm{w}}\in\mathbb{R}^{n}$ to $\tilde{\bm{w}}$ at every timestep of the integration. The cost of this lifting operation dominates the computational cost of the simulation and prevents its acceleration. To address this issue, a hyper-reduction method is employed to derive a hyper-reduced model that can be simulated without requiring lifting to the full state vector $\tilde{\bm{w}}$.

In this work we employ the EQP hyper-reduction method, which uses reduced quadrature rules to evaluate the nonlinear force functions with reduced cost [50, 51]. Quadrature methods such as EQP are naturally well suited to FEM spatial discretization procedures, such as the one employed in this work. In that respect, the present work can be viewed as a continuation of [22, 23], where interpolation methods were used to derive hyper-reduced Lagrangian hydrodynamics models.

3.1 Empirical quadrature procedure

In this section we describe the EQP hyper-reduction method using the model equations (3.5). The computation of the force vectors $\bm{F}_{v}^{\psi}$ and $\bm{F}_{e}^{\psi}$ of (3.5) requires the spatial integrations

[TABLE]

with the integrands $\bm{g}_{v}^{\psi}$ and $\bm{g}_{e}^{\psi}$ taking values in $\mathbb{R}^{n_{v}}$ and $\mathbb{R}^{n_{e}}$ respectively, with components

[TABLE]

These integrations can be approximated using numerical quadrature rules employing the underlying FEM basis functions [57, 60]. Let $\{x^{v}_{j},\rho^{v}_{j}\}_{j=0}^{J_{v}-1}$ and $\{x^{e}_{j},\rho^{e}_{j}\}_{j=0}^{J_{e}-1}$ denote the full quadrature rules ($x_{j}$ denoting points, $\rho_{j}\geq 0$ weights) for velocity and energy, with $J_{v}$ and $J_{e}$ their respective sizes. We can approximate the integrations by

[TABLE]

The computational cost of the above summations is proportional to the sizes of the full quadrature rules ($J_{v}$ and $J_{e}$), since the integrands must be evaluated at every point of the corresponding quadrature rule. In addition, since the integrands depend on the state variables they must be updated at every stage of the time integration.

To remove the dependence of the computational cost on the full sizes $J_{v}$ and $J_{e}$ we employ the EQP hyper-reduction method [50, 51]. EQP is used to generate quadrature rules of reduced size that can approximate the evaluation of the force functions with reduced computational cost while maintaining a reasonable level of accuracy. More specifically, EQP seeks to construct new quadrature rules $\{x^{v}_{j},\tilde{\rho}^{v}_{j}\}_{j=0}^{J_{v}-1}$ and $\{x^{e}_{j},\tilde{\rho}^{e}_{j}\}_{j=0}^{J_{e}-1}$ where only $\hat{J}_{v}\ll J_{v}$ and $\hat{J}_{e}\ll J_{e}$ of the new weights $\tilde{\rho}^{v}_{j}$ and $\tilde{\rho}^{e}_{j}$ are nonzero. In this way, computing the sums in (3.13) requires the evaluation of the integrands at only $\hat{J}_{v}$ and $\hat{J}_{e}$ sampled quadrature points respectively, which are the points corresponding to nonzero weights.

Note that EQP does not attempt to identify new quadrature points; rather, it samples a subset of the original quadrature points with appropriately modified weights. In that sense EQP can be thought of as a quadrature analog of interpolation hyper-reduction methods based on gappy POD [61]. The problem of constructing the reduced quadrature rules is posed as a linear optimization problem [50]. To outline the employed procedure we use the construction of a reduced quadrature rule for the approximation of $\bm{F}_{v}^{\psi}$ as an example. The construction of a reduced rule for $\bm{F}_{e}^{\psi}$ can be performed analogously.

Starting with the full quadrature rule $\{x^{v}_{j},\rho^{v}_{j}\}_{j=0}^{J_{v}-1}$ we want to find weights $\{\tilde{\rho}^{v}_{j}\}_{j=0}^{J_{v}-1}$ that minimize the 1-norm $\sum_{j=0}^{J_{v}-1}\tilde{\rho}^{v}_{j}$ subject to the nonnegativity constraints $\tilde{\rho}^{v}_{j}\geq 0$, $j\in\{0,\ldots,J_{v}-1\}$, and a total of $N_{c}^{v}\ll J_{v}$ accuracy constraints

[TABLE]

where $0\leq k<N_{t}$ indexes the time snapshots $\tilde{w}(t_{k})$ used to form the constraints, $0\leq i<n_{v}$ indexes the components of vector $\bm{g}_{v}^{\psi}$, and $0\leq s<N_{c}^{v}$ is the overall index of the accuracy constraints with associated error thresholds $\epsilon_{s}>0$. This implies a total of $N_{c}^{v}=n_{v}N_{t}$ constraints, corresponding to $n_{v}$ constraints (one for each component of $\bm{g}_{v}^{\psi}$) per included snapshot of the state variables $\tilde{w}(t_{k})$. In addition to selecting the snapshots $\tilde{w}(t_{k})$ that will be used to form the accuracy constraints, selecting appropriate values for the error thresholds $\epsilon_{s}$ is also important in constructing an effective reduced quadrature rule. In practice one seeks to strike a balance between sparsity (low number of nonzero weights, $\hat{J}_{v}$) and accuracy (low error thresholds $\epsilon_{s}$).

The linear optimization problem of forming the reduced quadrature rule can be recast as an equivalent nonnegative least squares (NNLS) problem [52, 62]. Details on the steps we take to form the NNLS problem and solve it are given in A.

3.2 Basic EQP hyper-reduction

In this section we employ the EQP method as outlined above to derive a hyper-reduced version of model (3.5). In particular, we use the reduced quadrature rules $\{x^{v}_{j},\tilde{\rho}^{v}_{j}\}_{j=0}^{J_{v}-1}$ and $\{x^{e}_{j},\tilde{\rho}^{e}_{j}\}_{j=0}^{J_{e}-1}$ to define the respective reduced force functions $\hat{\bm{F}}_{v}^{\psi}:\mathbb{R}^{n}\to\mathbb{R}^{n_{v}}$ and $\hat{\bm{F}}_{e}^{\psi}:\mathbb{R}^{n}\times\mathbb{R}^{n_{v}}\to\mathbb{R}^{n_{e}}$

[TABLE]

and semidiscrete hyper-reduced model

[TABLE]

The evaluation of the reduced force functions requires only $\hat{J}_{v}\ll J_{v}$ and $\hat{J}_{e}\ll J_{e}$ evaluations of the respective integrands. As a result, it requires knowledge only of the degrees of freedom of $\tilde{w}$ sampled by the reduced quadrature rules. To make this reduction in computational cost explicit in our notation, the reduced force functions $\hat{\bm{F}}_{v}^{\psi}$ and $\hat{\bm{F}}_{e}^{\psi}$ are written as functions of the reduced state $\hat{\bm{w}}$, instead of the lifted state $\tilde{\bm{w}}$. For the evolution of $\hat{\bm{x}}$ in (3.15), the quantities $\bm{\Phi}_{x}^{\top}\bm{v}_{os}\in\mathbb{R}^{n_{x}}$ and $\bm{\Phi}_{x}^{\top}\bm{\Phi}_{v}\in\mathbb{R}^{n_{x}\times n_{v}}$ can be precomputed, leading to an evolution equation that depends on reduced sizes only.

3.3 Temporal discretization

Applying the RK2A scheme to the hyper-reduced system (3.15), we arrive at the discrete hyper-reduced system

[TABLE]

with $\bar{\hat{\bm{v}}}_{k+1/2}=(\hat{\bm{v}}_{k}+\hat{\bm{v}}_{k+1})/2$. In the above, the timestep control algorithm is enforced based on the lifted state vector $\tilde{\bm{w}}$, using again only the degrees of freedom sampled by the reduced quadrature rules. In general, this leads to timesteps that are different from those used in the full order model, where the timestep control algorithm is enforced using the full state vector $\bm{w}$.

3.4 Time windowing

For unsteady problems dominated by advection, like the Euler equations considered in the present work, the dimension of the linear subspaces required to capture the relevant flow features generally grows with the simulation time. This is both because new dynamical behavior can emerge as the final simulation time is increased, and because dynamically similar flow features can appear in different regions of the spatial domain. As a result, attempting to build a single reduced model for the whole simulation time horizon may require a large reduced basis dimension to achieve the desired degree of accuracy. In turn, this can potentially eliminate any acceleration of the simulation offered by the employed reduction method.

To address this issue we use the method of time windowing [20, 21, 22, 23]. With time windowing, we divide the considered simulation time interval $[t_{0},t_{f}]$ into $N_{w}$ subintervals (windows) $[t_{i},t_{i+1}]$, $0\leq i<N_{w}$, and build a reduced model for each time window. More specifically, for each window we collect the solution snapshots that correspond to the prescribed time interval and use them to build reduced basis matrices and reduced quadrature rules that are local to the chosen time interval. This allows us to build reduced bases and quadrature rules that are of modest size and offer high accuracy in approximating the full solution for each time subinterval.

3.5 Simulation stages

The overall simulation process consists of two main stages: the offline and online stages. In the offline stage we perform numerical simulations using the full model (2.8) and collect the desired solution snapshots. Using those snapshots, we derive the reduced basis matrices and reduced quadrature rules. When multiple time windows are used, this process is carried out for each time window independently as described earlier. These steps complete the construction of the hyper-reduced model.

In the online stage we use the derived hyper-reduced model to perform numerical simulations that depend on reduced sizes only, thereby offering a reduction in computational cost. In the present work we focus on reproductive simulations, where the derived reduced model is used to reproduce simulation data that has been used for its training. This is in contrast to predictive simulations, where the reduced model is tested in predicting previously unseen simulation data.

4 Energy conservative EQP

In this section we develop a modified version of the basic EQP method that enforces the additional constraint of total energy conservation in the discrete hyper-reduced model. To do so, we use the semidiscrete system (3.7). Our goal is to derive reduced force functions $\hat{\bm{F}}_{v}^{\phi}:\mathbb{R}^{n}\to\mathbb{R}^{n_{v}}$ and $\hat{\bm{F}}_{e}^{\phi}:\mathbb{R}^{n}\times\mathbb{R}^{n_{v}}\to\mathbb{R}^{n_{e}}$ to approximate the force function evaluations

[TABLE]

with reduced computational cost. Using the reduced force functions, we can then derive the semidiscrete hyper-reduced model

[TABLE]

and apply the RK2A scheme to reach the discrete hyper-reduced system

[TABLE]

In the following sections we define the reduced force functions $\hat{\bm{F}}_{v}^{\phi}$ and $\hat{\bm{F}}_{e}^{\phi}$. We begin by identifying sufficient conditions for total energy conservation in the strong sense. Then we describe the implementation of the resulting hyper-reduction scheme used to incorporate these conditions.

4.1 Energy conservation

The total energy of the system is the sum of its internal and kinetic energy

[TABLE]

which can be represented using the FEM coefficient vectors

[TABLE]

The semidiscrete full model (2.3) conserves the total energy by design: $\frac{d}{dt}TE(\bm{w}(t))=0$ for all $t\geq 0$. In addition, the RK2A discrete scheme (2.8) conserves the discrete total energy: $TE(\bm{w}_{k+1})=TE(\bm{w}_{k})$ for all $k\geq 0$ [57]. In the following theorem we derive sufficient conditions that ensure that the conservation of discrete total energy is preserved in the discrete hyper-reduced model (4.2).

Theorem 1

Assume that the following three conditions are satisfied.

$\bm{v}_{os}=\bm{0}_{\mathcal{V}}$** 2. 2.

$\bm{1}_{\mathcal{E}}=\bm{\Phi}_{e}\hat{\bm{1}}_{e}$* for some $\hat{\bm{1}}_{e}\in\mathbb{R}^{n_{e}}$* 3. 3.

$\hat{\bm{v}}^{\top}\hat{\bm{F}}_{v}^{\phi}(\hat{\bm{w}})=\hat{\bm{1}}_{e}^{\top}\hat{\bm{F}}_{e}^{\phi}(\hat{\bm{w}},\hat{\bm{v}})$* for all $\hat{\bm{w}}\in\mathbb{R}^{n}$.*

Then the discrete hyper-reduced model (4.2) conserves the discrete total energy; namely,

[TABLE]

for all $k\geq 0$.

Proof 1

It holds true that

[TABLE]

The three conditions presented in Theorem 1 are motivated by the physics and structure of the hyper-reduced model. The first two conditions ensure that the kinetic energy $KE$ and internal energy $IE$ respectively can be represented in terms of the reduced state vectors. The third condition is the main structure preserving property coupling the velocity and energy forcing terms, motivated by the fact that the fluid flow governed by the discrete hyper-reduced model (4.2) is driven by pressure-volume work.

4.2 Implementation

To implement the energy conservative EQP method we begin with the semidiscrete model (3.7). We use zero offset vectors for position, velocity and energy in (3.1), which ensures we satisfy condition 1 of Theorem 1. To meet condition 2 we enrich the energy basis matrix $\bm{\Phi}_{e}$ by adding the energy unit $\bm{1}_{\mathcal{E}}$ as a column vector. Moreover, we orthonormalize the basis matrices $\bm{\Phi}_{v}$ and $\bm{\Phi}_{e}$ such that the reduced mass matrices (3.8) reduce to identity matrices.

To satisfy condition 3 we are going to derive one combined reduced quadrature rule to be used for both the velocity $\bm{F}_{v}^{\phi}$ and energy $\bm{F}_{e}^{\phi}$ force functions. More specifically, we are going to form a reduced quadrature rule for the force function $\bm{F}^{\phi}:\mathbb{R}^{N}\to\mathbb{R}^{n_{v}\times n_{e}}$

[TABLE]

where $\bm{G}$ denotes the integrand taking values in $\mathbb{R}^{n_{v}\times n_{e}}$, with entries

[TABLE]

With respect to $\bm{F}^{\phi}$, the velocity and energy force functions are given by

[TABLE]

We denote by $\{x_{j},\rho_{j}\}_{j=0}^{J-1}$ the full quadrature rule that can be used to approximate the integration (4.3)

[TABLE]

Starting with this rule, we employ the EQP method to derive a reduced quadrature rule $\{x_{j},\tilde{\rho}_{j}\}_{j=0}^{J-1}$ with $\hat{J}\ll J$ nonzero weights, and use it to define the reduced force function $\hat{\bm{F}}^{\phi}:\mathbb{R}^{n}\to\mathbb{R}^{n_{v}\times n_{e}}$

[TABLE]

The evaluation of $\bm{F}_{v}^{\phi}$ and $\bm{F}_{e}^{\phi}$ is then approximated with reduced computational cost by

[TABLE]

It follows that for all $\hat{\bm{w}}\in\mathbb{R}^{n}$

[TABLE]

satisfying condition 3 of Theorem 1. Therefore, we conclude that the discrete hyper-reduced model (4.2) conserves the discrete total energy for all time indices $k\geq 0$. Additional details on the implementation of the energy conservative EQP method are given in A.

5 Numerical experiments

In this section we present numerical results for four problem cases used to compare the performance of the basic and energy conservative EQP methods. Using each of the two EQP methods we simulate the following problems: the 3D Sedov blast problem, the 2D Gresho vortex problem, the 3D triple point problem, and the 3D Taylor-Green vortex problem. They are all problems that have been used in the past as benchmarks for testing the performance of different discretization and reduction methods [57, 22, 61].

5.1 Problems

In the following paragraphs we provide an overview of each problem case and its simulation parameters. The spatial domain in each problem is discretized using an initially uniform, cartesian mesh of hexahedral elements. Table 1 presents discretization parameters used for each of the simulated problems, including the number of basis functions used for each finite element space and the corresponding polynomial basis orders.

Sedov blast

The Sedov blast problem is a three dimensional generalization of the two dimensional problem considered in [57]. The problem simulates the sudden release of energy from a localized source, leading to the formation of a shock wave propagating in the spatial domain and converting the initial internal energy into kinetic one. Its original goal was to test the ability of different artificial viscosity tensors to maintain radial symmetry of the propagating shock wave on an irregular discretization mesh.

The computational domain is $\Omega_{0}=[0,1]^{3}$. We use a constant adiabatic index $\gamma=1.4$, initial density $\rho=1$, and zero initial velocity. The initial energy has magnitude $|e|=0.25$ and is all concentrated at the domain origin. We employ one of the artificial viscosity stress tensors developed in [57] that preserves the propagating shock’s radial symmetry. The simulation time interval is $[0,0.3]$.

Gresho vortex

The Gresho vortex problem is a two dimensional problem originally posed using the incompressible Euler equations [22]. Here we use the steady state solution of the original incompressible problem to build a manufactured solution of the compressible equations. We do that by substituting the original steady state solution into the compressible Euler equations and deriving the source terms required to maintain that steady state [57, 22]. The main challenge posed by this problem case is the large deformations induced on the spatial grid as time evolves.

The computational domain is $\Omega_{0}=[-0.5,0.5]^{2}$. For this problem we use $\gamma=5/3$, initial density $\rho=1$, initial pressure

[TABLE]

where $(r,\phi)$ denote the polar coordinates, and zero artificial viscosity $\sigma_{a}=0$. Additionally, we use zero initial radial velocity $v_{r}=0$ and the initial angular velocity

[TABLE]

Based on the above, the initial energy is determined by the thermodynamic equation of state (2.2). The simulation time interval is $[0,0.4]$.

Triple point problem

The triple point problem is a three dimensional Riemann problem involving two materials in three different states. The material interaction leads to the formation of shock waves and generation of vorticity [57, 22].

The computational domain is $\Omega_{0}=[0,7]\times[0,3]\times[0,1.5]$. We use the adiabatic index

[TABLE]

zero initial velocity $v=0$, initial density

[TABLE]

and initial pressure

[TABLE]

The initial energy is computed from the thermodynamic equation of state (2.2). Given the formation of shock waves, we also use an artificial viscosity tensor [57]. The simulated time interval is $[0,0.8]$.

Taylor-Green vortex

The Taylor-Green vortex is a three dimensional generalization of a two dimensional steady state solution to the incompressible Euler equations [57, 22]. Similarly to the Gresho vortex case, we build a manufactured solution of the compressible equations by deriving the source terms necessary for maintaining the original incompressible steady state solution.

The computational domain is $\Omega_{0}=[0,1]^{3}$. We use $\gamma=5/3$, initial density $\rho=1$, initial pressure

[TABLE]

and initial velocity

[TABLE]

The initial energy is computed using the thermodynamic equation of state (2.2). We use no artificial viscosity stress. The simulated time interval is $[0,0.4]$.

5.2 Results

Figure 1 presents the final velocity snapshots obtained using the full model and the two reduced models for each problem case. Table 2 presents the numerical results for the four problem cases using the two EQP hyper-reduction methods compared in this work.

Starting with the relative energy difference obtained for each simulated problem, we see that the conservative EQP (CEQP) method outperforms basic EQP (BEQP) by several orders of magnitude. The difference in performance is about 9 orders of magnitude for the Sedov case, 6 orders of magnitude for Gresho, 8 orders of magnitude for triple point and 3 orders of magnitude for Taylor-Green. In addition, the energy difference achieved with CEQP is near machine precision for all four simulated problems. These results prove that the changes made in the derivation of the reduced quadrature rule do indeed lead to high accuracy in the conservation of energy, as expected from the theoretical analysis.

In terms of the differences in relative errors achieved for the position, velocity and energy solutions, the two methods perform similarly in Sedov and triple point but demonstrate more pronounced discrepancies in Gresho and Taylor-Green. More specifically, for the Sedov case both methods lead to similar errors, with CEQP outperforming BEQP in position and velocity and conversely in energy. Similarly, the two methods yield similar results for triple point, with BEQP performing better by about an order of magnitude in all three variables. The situation is different for Taylor-Green and Gresho. For Taylor-Green, CEQP underperforms by about 2-3 orders of magnitude in the three solution errors, while for Gresho by about 3-4 orders of magnitude. Despite the more pronounced differences in accuracy in these problem cases, the errors obtained by CEQP are still reasonably low.

Finally, we arrive at the comparison of the simulation speedup (acceleration) values obtained with each EQP method. We see that the two methods yield similar speedup results, although BEQP is outperforming CEQP in all four problem cases. The moderate differences in speedup demonstrate that the changes introduced into the derivation of the reduced quadrature rule for CEQP has not hurt that aspect of the method’s performance.

5.3 Discussion

Among the presented numerical results, the discrepancy in solution accuracy between the basic and energy conservative EQP methods is largest in the Gresho problem case. In addition, Gresho is the problem case generating the largest values of solution errors for CEQP, despite the fact that its energy conservation performance is on par with the other experiments. As of this writing, it is not well understood why that is the case, and whether the elevated errors are due to the numerical implementation of the method or due to some intrinsic challenge posed by the Gresho problem.

The reported speedup factors are relatively small compared to those reported in [22], where interpolation methods are used for hyper-reduction. This is partly because our results are generated using a numerical implementation of the basic and energy conservative EQP methods that is not fully optimized yet. More specifically, the used implementation performs some of the lifting of the reduced state vectors to full state vectors using all degrees of freedom, instead of using only the degrees of freedom sampled by the reduced quadrature rules. Independently of those optimizations, we also expect that the speedup values will increase with increasing values for the mesh refinement parameter $m$ and the FEM basis polynomial order $k$.

A direction of future research is comparing the performance of different NNLS algorithms, especially as it relates to the sparsity and accuracy of the derived reduced quadrature rules. As shown in Table 2, the Lawson-Hanson algorithm used in this work leads to sparse quadrature rules. Nevertheless, it is not clear how different methods would compare against our reported results. In addition, we intend to investigate why the $LQ$ preconditioning of the NNLS constraints matrix (ensuring orthogonality of constraints) leads to improved performance of the Lawson-Hanson algorithm. Although observed empirically, to our knowledge it is not fully understood why that happens [63].

In the present work we focus on reproductive simulations, where each reduced model is tested using initial conditions that have been used in its training. The developed methods readily generalize to predictive simulations, where out of sample initial conditions are used to test the model’s performance. We intend to pursue this in the near future.

Conceptually, the model reduction approach employed in this work builds upon the finite element framework. Starting with the FEM discretization of the problem, with the FEM basis functions acting as the trial and test functions, we proceed to replace these functions by ones derived from simulation data. One important difference between the two model formulations is that the FEM model employs trial and test functions that are local in space, meaning that their support is confined within a small number of elements. On the contrary, the reduced models employ trial and/or test functions with a support that is generally extended across a large number of elements, possibly including the whole spatial domain. Although this allows us to use a smaller number of degrees of freedom to represent the considered dynamics, it can also lead to some of the extrapolation problems encountered by reduction methods to date, since the accurate representation and prediction of out of sample dynamics cannot be guaranteed.

There is ongoing work related to deriving data based reduction methods that use spatially local reduced basis functions, bringing them closer to a data based analog of the finite element framework [64, 65, 66, 67, 68]. It is interesting to explore the applicability of such methods to the problem setting considered in the present work, especially since quadrature hyper-reduction methods are naturally well suited to the finite element framework.

6 Conclusion

We developed reduced models for nonlinear Lagrangian hydrodynamics problems using quadrature hyper-reduction to approximate the evaluation of the nonlinear terms with reduced computational cost. The hyper-reduction was performed using the EQP method, which is well suited to the finite element discretization of the considered nonlinear conservation laws.

Starting with an application of the basic EQP method, we developed a variant of the method that preserves the energy conservation of the full model in the reduction process. In addition to proving exact energy conservation in the resulting reduced model, we performed numerical simulations of four different benchmark problems demonstrating the conservation of energy to near machine precision. Our reported results show that the additional benefit of energy conservation is obtained with only a moderate reduction in the accuracy or speedup of the reduced simulations compared to the basic EQP method.

Acknowledgments

This work was performed under the auspices of the U.S. Department of Energy (DOE) by Lawrence Livermore National Laboratory (LLNL) under Contract No. DE-AC52–07NA27344. Part of this work was performed while C. Vales was a research intern at LLNL supported by the Mathematical Sciences Graduate Internship (MSGI) program administered by the U.S. National Science Foundation. C. Vales also acknowledges support from the U.S. Department of Energy under grant DE-SC0025101. Y. Choi was supported for this work by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, as part of the CHaRMNET Mathematical Multifaceted Integrated Capability Center (MMICC) program, under Award Number DE-SC0023164. IM release: LLNL-JRNL-2010206.

Appendix A Implementation details

In this section we offer additional details on the implementation of the basic and conservative EQP methods developed in the main text.

A.1 Basic EQP

We consider the nonnegative least squares (NNLS) formulation used for the basic EQP method developed in Section 3.1, as well as the preconditioning operations performed before solving the formed NNLS problem.

NNLS formulation

To solve the linear optimization problem for the reduced quadrature rule described in Section 3.1, we reformulate it as an NNLS problem for the vector of quadrature weights $\tilde{\bm{\rho}}^{v}\in\mathbb{R}^{J_{v}}$ [52, 62]. Given the known vector of full quadrature weights $\bm{\rho}^{v}\in\mathbb{R}^{J_{v}}$, we form the accuracy constraints matrix $\bm{C}_{v}\in\mathbb{R}^{N_{c}^{v}\times J_{v}}$, $N_{c}^{v}\ll J_{v}$, with entries

[TABLE]

with row index $0\leq s=i+kn_{v}<N_{c}^{v}$. Denoting $\bm{b}=\bm{C}_{v}\bm{\rho}^{v}\in\mathbb{R}^{N_{c}^{v}}$, the vector $\tilde{\bm{\rho}}^{v}$ is found by solving the NNLS problem

[TABLE]

More specifically, for a given vector of error thresholds $\bm{\epsilon}\in\mathbb{R}^{N_{c}}$, we use the iterative Lawson-Hanson algorithm to solve the following problem for $\bm{r}\geq\bm{0}$

[TABLE]

where $\bm{c}_{v}^{s}\in\mathbb{R}^{J_{v}}$ denotes the $s$th row of matrix $\bm{C}_{v}$, and $b_{s}$, $\epsilon_{s}$ the $s$th respective entries of vectors $\bm{b}$ and $\bm{\epsilon}$ [69, 70].

In the NNLS formulation, (A.1) is used to enforce the accuracy constraints of the linear optimization problem. The nonnegativity constraints are enforced implicitly by the employed algorithm. Finally, compared to the original form of the problem, the NNLS formulation does not explicitly promote sparsity by minimizing the 1-norm of the solution vector $\tilde{\bm{\rho}}^{v}$. Nevertheless, several numerical experiments—including the ones we report in this work—show that sparse solutions are obtained by the Lawson-Hanson algorithm [52, 62, 61].

Preconditioning

Before solving (A.1), we perform two conditioning operations on the constraints matrix $\bm{C}_{v}$, both to reduce the number of iterations required by the Lawson-Hanson algorithm to converge and to increase the derived solution’s accuracy [63]. First, we rescale each row of $\bm{C}_{v}$ so that its entries have absolute value less than or equal to one. Second, we use the reduced LQ decomposition of $\bm{C}_{v}=\bm{L}_{v}\bm{Q}_{v}$, where $\bm{L}_{v}\in\mathbb{R}^{N_{c}^{v}\times N_{c}^{v}}$ is lower triangular and $\bm{Q}_{v}\in\mathbb{R}^{N_{c}^{v}\times J_{v}}$ has orthonormal rows, to transform the NNLS problem (A.1) to

[TABLE]

where $\tilde{\bm{b}}=\bm{Q}_{v}\bm{\rho}^{v}\in\mathbb{R}^{N_{c}^{v}}$, under the assumption that $\bm{L}_{v}$ is nonsingular. In addition, we modify the error thresholds $\epsilon_{s}$ accordingly to ensure that the solution obtained from (A.3) will also satisfy (A.2) with the original error parameters. We define the new error thresholds $\tilde{\epsilon}_{s}$ by

[TABLE]

which is sufficient to ensure that the transformed problem’s solution satisfies (A.2) [63]. The transformed NNLS problem reads

[TABLE]

where $\bm{q}_{v}^{s}\in\mathbb{R}^{J_{v}}$ denotes the $s$th row of matrix $\bm{Q}_{v}$.

A.2 Conservative EQP

We consider the implementation of the NNLS problem for the energy conservative EQP method developed in Section 4, as well as an additional basis enrichment operation performed to ensure conservation of the discrete total energy for simulations consisting of multiple time windows.

NNLS formulation

Following the approach outlined in Section 3.1 and A.1, the reduced quadrature rule is computed by solving an NNLS problem with a total of $N_{c}=(n_{v}+n_{e})N_{t}$ accuracy constraints, where $N_{t}$ denotes the number of snapshots of $\tilde{w}$ used to form the constraints. The constraints are used to control the accuracy of a single reduced quadrature rule that will be used for the approximation of both force vectors $\bm{F}_{v}^{\phi}$ and $\bm{F}_{e}^{\phi}$ introduced in (3.10)

[TABLE]

where the respective integrands $\bm{g}_{v}^{\phi}$ and $\bm{g}_{e}^{\phi}$ take values in $\mathbb{R}^{n_{v}}$ and $\mathbb{R}^{n_{e}}$ and have components

[TABLE]

The accuracy constraints are encoded in the matrix $\bm{C}\in\mathbb{R}^{N_{c}\times J}$, where the first $n_{v}N_{t}$ rows are used to enforce the constraints related to $\bm{F}_{v}^{\phi}$

[TABLE]

with constraint index $s=i+kn_{v}$, and the next $n_{e}N_{t}$ rows the constraints related to $\bm{F}_{e}^{\phi}$

[TABLE]

with constraint index $s=n_{v}N_{t}+i+kn_{e}$. In the above, $\tilde{w}(t_{k})$ denotes the snapshots used to form the constraints.

Reduced bases enrichment

When the simulation consists of more than one time window, an additional modification is made to the reduced bases $\bm{\Phi}_{v}$ and $\bm{\Phi}_{e}$ to ensure that the total energy is conserved across windows. Let $\tilde{\bm{v}}^{(n-1)}$ and $\tilde{\bm{e}}^{(n-1)}$ denote the velocity and energy solution vectors after the final timestep of window $n-1$, $n>1$, has been completed

[TABLE]

where superscripts are used to mark the time window index. Accordingly, we denote by $\tilde{\bm{v}}^{(n)}$ and $\tilde{\bm{e}}^{(n)}$ the corresponding vectors after the change from window $n-1$ to window $n$ has been performed

[TABLE]

Note that no timestepping has been performed after switching from window $n-1$ to $n$, so the two pairs of vectors should be identical. To ensure that this is the case, we extend the basis matrices $\bm{\Phi}_{v}^{(n)}$ and $\bm{\Phi}_{e}^{(n)}$ by adding the respective solution vectors $\tilde{\bm{v}}^{(n-1)}$ and $\tilde{\bm{e}}^{(n-1)}$ as an additional column and reorthonormalize. This guarantees that $\tilde{\bm{v}}^{(n-1)}$ and $\tilde{\bm{e}}^{(n-1)}$ are in the span of the respective basis matrices $\bm{\Phi}_{v}^{(n)}$ and $\bm{\Phi}_{e}^{(n)}$. In turn, this ensures that energy is conserved exactly when changing windows

[TABLE]

No such extension is required for the first window ($n=0$). Since the extension of the basis matrices $\bm{\Phi}_{v}^{(n)}$ and $\bm{\Phi}_{e}^{(n)}$ requires knowledge of the state vectors $\tilde{\bm{v}}^{(n-1)}$ and $\tilde{\bm{e}}^{(n-1)}$, this operation must be performed during the online stage at the time of the window change. To extend and incrementally reorthonormalize the basis matrices efficiently, we employ the modified Gram-Schmidt algorithm with double orthogonalization [71].