Entropy stable discontinuous Galerkin methods for nonlinear conservation laws on networks and multi-dimensional domains

We begin with introducing the two-dimensional shallow water equations

Here, $h$ denotes the water height as measured from the lake or channel bottom. The velocity in the $x$ direction is denoted by $u$ and the velocity in the $y$ direction is denoted by $v$. The gravitational constant is denoted by $g$. In this example, we have the conservative variable $\bm{u}=[h,hu,hv]^{T}$ and the flux functions are $\bm{f}_{1}=[hu,hu^{2}+gh^{2}/2,huv]^{T}$ and $\bm{f}_{2}=[hv,huv,hv^{2}+gh^{2}/2]^{T}$.

We can derive the 1D shallow water equations from the 2D shallow water equations by assuming a rectangular domain with length $L_{x}$ and $L_{y}$ in the $x$ and $y$ directions, respectively. If $L_{y}\ll L_{x}$ and wall boundary conditions are imposed, then we expect $v$ to be small and $h,u$ to be near-constant along the $y$-direction. These simplifications result in the one-dimensional shallow water equations

The mathematical entropy for the shallow water equations corresponds to total energy, and is given by

where $\left\|U\right\|^{2}=u^{2}$ in one dimension and $\left\|U\right\|^{2}=u^{2}+v^{2}$ in two dimensions. The entropy variables $\bm{v}$ in two dimensions are given by

In one dimension, the entropy variables are simply $\bm{v}=\left[v_{1},v_{2}\right]^{T}$. The inverse mapping in 2D is given by

where we can compute $\left\|U\right\|^{2}=v_{2}^{2}+v_{3}^{2}$ in terms of the entropy variables. The inverse mapping in 1D follows by ignoring $hv$ and setting $\left\|U\right\|^{2}=v_{2}^{2}$.

In this paper, we only consider systems where the conservative variables, the entropy potential, and the numerical fluxes can all be transformed between 1D and 2D.

We adopt a notation which distinguishes between discretized and continuous quantities. Unless otherwise specified, continuous vector and matrix quantities are denoted using lower and upper case bold font, respectively. We denote spatially discrete quantities using a bold sans serif font. Finally, the output of continuous functions evaluated over discrete vectors is interpreted as a discrete vector.

For example, if $\bm{\mathsf{x}}$ denotes a vector of point locations, i.e., $(\bm{\mathsf{x}})_{i}=\bm{x}_{i}$, then $u(\bm{\mathsf{x}})$ is interpreted as the vector

Similarly, if $\bm{\mathsf{u}}={u}(\bm{\mathsf{x}})$, then ${f}(\bm{\mathsf{u}})$ corresponds to the vector

Vector-valued functions are treated similarly. For example, given a vector-valued function $\bm{f}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ and a vector of coordinates $\bm{\mathsf{x}}$, $\left(\bm{f}(\bm{\mathsf{x}})\right)_{i}=\bm{f}(\bm{x}_{i})$.

We construct entropy stable numerical schemes for networks based on high order entropy stable DG formulations in 1D and 2D. These formulations ensure entropy stability over individual segments and subdomains of a multi-dimensional network Abc (2000). We begin by introducing some mathematical notation. We denote the reference element by $\widehat{D}$ with boundary $\partial{\widehat{D}}$. In 1D, the reference element is the interval $[-1,1]$ and in 2D, the reference element is the bi-unit right triangle. We construct entropy conservative schemes on multiple elements, where the domain $\Omega$ is broken up into $K$ non-overlapping elements $D^{k}$. Each element can be represented as the affine mapping $\Phi_{k}$ of the reference interval $\widehat{D}$. Because this mapping is affine, $J^{k}$ (the determinant of the Jacobian of $\Phi_{k}$) is constant over each element. We use $\widehat{n}_{i}$ to represent the $i$th component of the outward normal vector scaled by the face Jacobian on the boundary of the reference element.

The solution is approximated over the reference element by polynomials of total degree $N$

We denote the dimension of $P^{N}$ as $N_{p}=\rm{dim}(P^{N}(\widehat{D}))$. Moreover, let $\{\phi_{i}\}_{i=1}^{N_{p}}$ denote a basis for $P^{N}$, such that for $u(\bm{x})\in P^{N}(\widehat{D})$, there exist coefficients $u_{i}$ such that

We also assume the use of volume and surface quadrature rules $\left\{\widehat{\bm{x}}_{i},w_{i}\right\}_{i=1}^{N_{q}}$, $\left\{\widehat{\bm{x}}^{f}_{i},w^{f}_{i}\right\}_{i=1}^{N^{f}_{q}}$. We denote the number of volume and surface quadrature nodes by $N_{q}$ and $N^{f}_{q}$ respectively, and assume that the volume and surface rules are exact for polynomials of degree $2N-1$ and $2N$, respectively. We furthermore assume the volume quadrature is sufficiently accurate such that the mass matrix is positive-definite.

We now introduce quadrature-based operators. Let $\bm{\mathsf{W}}$ and $\bm{\mathsf{W}}_{f}$ denote diagonal matrices whose entries are volume and surface quadrature weights. We then define the volume and surface quadrature interpolation matrices $\bm{\mathsf{V}}_{q}$ and $\bm{\mathsf{V}}_{f}$ as:

The matrix $\bm{\mathsf{V}}_{q}$ maps coefficients of $\bm{\mathsf{u}}=\left[u_{1},u_{2},\ldots,u_{N_{p}}\right]$ in terms of polynomial basis to evaluations of $u(\bm{x})$ at volume quadrature points and, similarly, the matrix $\bm{\mathsf{V}}_{f}$ interpolates $\bm{\mathsf{u}}$ to surface quadrature points.

We now define $\bm{\mathsf{D}}_{i}$ as the differentiation matrix with respect to the $i$th coordinate. $\bm{\mathsf{D}}_{i}$ is defined implicitly by:

Here, $\bm{\mathsf{D}}_{i}$ maps the basis coefficients of some polynomial $u\in P^{N}$ to coefficients of its $i$th directional derivative with respect to the reference coordinate $\widehat{\bm{x}}_{i}$.

With the matrix $\bm{\mathsf{V}}_{q}$, we can now introduce the element mass matrix whose entries are the evaluations of inner products of different basis functions with quadrature points:

We can define the quadrature-based $L^{2}$ projection matrix $\bm{\mathsf{P}}_{q}$, by inverting the mass matrix:

The matrix $\bm{\mathsf{P}}_{q}$ maps a function in terms of its evaluations at quadrature points to its coefficients of the $L^{2}$ projection in the basis $\phi_{i}(\bm{\widehat{x}})$.

In d-dimensions, define the following matrices

With the above definitions, we have can show that Abc (2000)

By combining the projection matrix $\bm{\mathsf{P}}_{q}$ with the matrix $\bm{\mathsf{\widehat{Q}}}_{i}$, we can construct a nodal differentiation operator at quadrature points Abc (2000):

We also define the the matrix $\bm{\mathsf{E}}$, which extrapolates volume quadrature nodes to surface quadrature nodes, as

Then we have the following generalized summation-by-parts (SBP) property:

Entropy stable formulations for nonlinear conservation laws usually introduce numerical flux terms which couple together all degrees of freedom on neighboring elements Abc (2000). To avoid this, we introduce the hybridized operator $\bm{\mathsf{Q}}_{i,h}$, which is given explicitly as

This operator is designed to be applied to vectors of solution values at both volume and surface quadrature nodes and mimics the structure of boundary terms used in hybridized DG methods Abc (2000). When used in a DG formulation, it allows one to construct entropy stable formulations using more standard DG numerical fluxes. We have the following theorem:

We also construct differentiation and boundary matrices $\bm{\mathsf{Q}}^{k}_{i,h},\bm{\mathsf{B}}^{k}_{i}$ on the physical element $D^{k}$ through the chain rule and linear combinations of differentiation matrices on the reference element. These will be used to construct entropy conservative and entropy stable schemes in the following section. Let $g^{k}_{ij}=J^{k}\frac{\partial x_{i}}{\partial\hat{x}_{j}}$ denote geometric terms on $D^{k}$. Then, physical SBP operators can be constructed by taking a linear combination of the reference SBP operators

It was shown in Abc (2000) that these operators satisfy a physical SBP property

where $\bm{\mathsf{B}}^{k}_{i}$ is a diagonal matrix containing the $i$th component of the outward unit normal on $D^{k}$ scaled by quadrature weights and surface Jacobian factors. We note that this construction uses the fact that $g^{k}_{ij}$ is constant over $D^{k}$ for affinely mapped elements. We also introduce the physical mass matrix $\bm{\mathsf{M}}^{k}=J^{k}\bm{\mathsf{M}}$, which is scaled by a Jacobian factor. Finally, we note that physical SBP operators on curved elements can be constructed using a “split form” version of (3.15) chan2019discretely.

Finally, we introduce $\bm{v}_{h}$ as the $L^{2}$ projection of the entropy variables and $\bm{\tilde{u}}$ as the evaluations of the conservative variables in terms of the $L^{2}$ projected entropy variables

Here $\bm{u}_{q}$ and $\bm{v}_{q}$ denote the conservative variables and entropy variables evaluated at the volume quadrature points. The vector $\bm{\tilde{v}}$ denotes the evaluations of the $L^{2}$ projection of the entropy variables at both volume and surface quadrature points, while $\bm{\tilde{u}}$ denotes the evaluations of the conservative variables in terms of the projected entropy variables $\bm{u}(\Pi_{N}\bm{v})$, where $\Pi_{N}$ denotes the $L^{2}$ projection operator.

In this section, we introduce entropy conservative numerical fluxes and discrete formulations Abc (2000, 2000, 2000, 2000). To construct entropy stable schemes in $d$ dimensions, we utilize entropy conservative fluxes as defined in Abc (2000):

The flux $\bm{f}_{i,S}$ can be used to construct entropy conservative and entropy stable finite volume methods. This numerical flux can also be used to construct discretely entropy stable DG schemes using an approach referred to as flux differencing Abc (2000, 2000, 2000, 2000).

Using flux differencing Abc (2000, 2000), we can approximate the derivative of $\bm{f}_{i}(\bm{u}(x))$ using the differentiation matrices $\bm{\mathsf{Q}}^{k}_{i,h}$ and $\bm{f}_{i,S}$. We define a flux matrix $\bm{\mathsf{F}}_{i}$ by evaluating $\bm{f}_{i,S}$ using solution values at quadrature points:

Then, $\int\frac{\partial\bm{f}_{i}(\bm{u})}{\partial x_{i}}$ can be approximated by the term $2(\bm{\mathsf{Q}}^{k}_{i,h}\circ\bm{\mathsf{F}}_{i})\bm{1}$, where $\circ$ denotes the Hadamard product between two matrices. A discrete entropy conservative formulation is then given as follows on an element $D^{k}$:

where $\bm{\mathsf{S}}$ denotes source terms depending on the bottom geometries.

Entropy stable schemes can be constructed from an entropy conservative scheme by adding appropriate penalization term that dissipate entropy at element interfaces Abc (2000). This modification converts schemes which satisfy a global entropy equality into schemes which satisfy a global entropy inequality.

We first introduce notation for 1D and 2D fluxes. Let $\bm{f}_{i,S}$ denote the 2D flux, where the subscript $i$ denotes the $i$th component of the 2D fluxes, (where $i=1,2$ correspond to the $x$ and $y$ direction respectively). We denote the 1D flux vector by $\bm{f}_{S}$, without the subscript for the coordinate index. We use $\bm{f}_{J,1D}$ and $\bm{f}_{J,2D}$ to denote the numerical flux at a multi-dimensional junction on 1D and 2D sides respectively. We use $\bm{f}_{J,i}$ to denote the fluxes at a 1D-1D junction, where the subscript $i$ denotes the numerical flux for the $i$th 1D channel. We use $n_{1D}$ and $n_{2D}$ to denote the sign of the outward surface normal on the 1D and 2D domains at a 1D-2D junction.

Finally, we restrict ourselves to velocity-based systems, where in $d$ dimensions, the conservative variables contain a velocity or momentum vector with $d$ components. We also assume that 1D system can be derived from 2D system by ignoring or projecting certain velocity components. We note our approach is not directly applicable to systems which do not satisfy these assumptions, such as Burger’s equation.

For flows over networks of 1D channels, one approach to model junction behavior is to explicitly discretize the junction using a 2D mesh Abc (2000, 2000). Since each junction can be accurately modeled using a relatively small number of elements, this approach inherits most of the computational advantages of a purely 1D network simulation. We illustrate the construction of a 1D-2D coupling scheme using the simplified setting shown in Figure 2. We assume that the 2D domain has an interface which is perpendicular to the 1D domain. This property can be enforced when constructing meshes for 2D domains. We couple 1D and 2D domains together by treating the 1D domain as a channel with the same width of the 2D domain.

In both 1D and 2D, each cell communicates through fluxes (calculated in terms of “interior” and “exterior” values of the conservative variables) at the interfaces, and we treat the coupling interface similarly. We first compute the fluxes at the interface of the coupling in the 2D domain with the following steps. Let $\bm{u}_{2D}$ denote the vector of the conservative variables in the “interior” of 2D cell, and let $\bm{u}^{+}_{2D}$ denotes the “exterior” values of the conservative variables from the neighbouring cell. The conservative variables in 1D have fewer fields than the variables in 2D. To resolve this mismatch, we introduce a matrix $\bm{R}$ to transform between 1D and 2D variables. For shallow water equations in 1D and 2D, we have

where $[n_{1},n_{2}]^{T}$ is the unit outward normal vector for the 1D domain with respect to the 2D coordinates. Then we can use $\bm{R}$ to transform between 1D and 2D variables as follows:

Note that $\bm{R}^{T}\bm{u}_{2D}$ constructs the $1D$ momentum from the normal component of the 2D momentum. Similarly, we can apply this $\bm{R}$ operator to the entropy conservative fluxes and entropy variables in both 1D and 2D.

Let $\bm{u}_{f,2D}$ and $\bm{u}_{f,1D}$ denote the vector of values of the conservative variables at surface quadrature points in the 2D and 1D domains on general elements and let $\bm{u}_{J,2D}$ and $\bm{u}_{J,1D}$ denote the vector of values of the conservative variables at surface quadrature points on the 2D and 1D side of a 1D-2D interface. Then the flux on the 2D side of the 1D-2D interface $\bm{f}_{J,2D}$ is defined as

We next apply the transformation matrix $\bm{R}^{T}$ to the 2D fluxes to reduce them down to 1D dimension. We define the surface flux for 1D side of the 1D-2D interface as follows:

where $\bm{\mathsf{w}}_{J,f}$ is the vector containing the quadrature weights of the surface quadrature points on the entire 1D-2D interface. Note that operation $\bm{R}^{T}\bm{f}_{J,2D}$ is performed point-wise at each surface quadrature point.

We now present a proof of entropy conservation for our 1D-2D coupling scheme. Let $K_{1D}$ and $K_{2D}$ denote the number of elements in the 1D and 2D domains respectively. Let $J^{k}_{1D}$ and $J^{k}_{2D}$ be the Jacobian on the $k$th elements in the 1D and 2D domains respectively. We first introduce the following lemma:

The proof of the lemma is restatement of results from Abc (2000). Next, assume a simple 1D-2D coupling as shown in Figure 2, where we assume the 1D and 2D domain are channels of the same width denoted by $A$. The 1D-2D junction terms can be expressed as

where $\bm{f}_{J,2D}(\bm{\mathsf{\tilde{u}}}_{J,1D},\bm{\mathsf{\tilde{u}}}_{J,2D})$ and $\bm{f}_{J,1D}(\bm{\mathsf{\tilde{u}}}_{J,1D},\bm{\mathsf{\tilde{u}}}_{J,2D})$ are defined in (4.5) and (4.4) respectively. Together, these 1D-2D discretizations imply entropy conservation:

We note that a 2D domain can be coupled to multiple 1D domains. For example, in Figure 3, we calculate the flux for Channel 1 with fluxes from the top two cells and the flux for Channel 2 with the bottom six cells of the 1D-2D interface.

In this section, we introduce entropy stable methods to join multiple (2 or more) 1D domains together without explicitly representing the junction as a 2D domain. This approach simplifies meshing and further reduces the overall computational cost, though it may not accurately reproduce multi-dimensional effects in certain flow setting Abc (2000, 2000).

We first present a general method to couple an arbitrary number of 1D domains at a junction. The key is to distribute shared fluxes between each domain connected to the junction in a conservative fashion. Let $I_{J}$ be a set of all domains at the junction and let $|I_{J}|$ denotes the dimension of the set. Let $A_{i}$ denote the width of the $i$th domain at a junction. Let $n_{J,i}=\pm{1}$ denote the outward surface normal for channel $i\in I_{J}$, as shown in Figure 4. We introduce on channel $i\in I_{J}$ the values of solution $\bm{u}_{f,i}$ and junction flux $\bm{f}_{J,i}$

where we have the weighting coefficients $c_{ij}$. We require these coefficients to satisfy

We introduce the convention that non-zero values of $c_{ii}$ represent “partial” wall boundary conditions (which we describe in more detail in Section 4.3.2), where the solution of the a domain partially interacts with itself.

Let $K_{i}$ denote the number of elements in the $i$th channel and let $J^{k}_{i}$ be the Jacobian of the $k$th element of the $i$th channel. With Lemma 4.1, we have the following theorem: