Sharp Penalty Term and Time Step Bounds for the Interior Penalty Discontinuous Galerkin Method for Linear Hyperbolic Problems

TL;DR

This paper derives sharp bounds for the penalty term and time step size in the SIPDG method to ensure stability across various linear hyperbolic problems and complex mesh configurations.

Contribution

It provides the first sharp, element-wise bounds for penalty and time step size applicable to a wide range of hyperbolic PDEs on unstructured meshes.

Findings

Bounds are sharp and numerically validated.

Applicable to diverse hyperbolic equations and mesh types.

Element-wise computation of bounds enhances practical implementation.

Abstract

We present sharp and sufficient bounds for the interior penalty term and time step size to ensure stability of the Symmetric Interior Penalty Discontinuous Galerkin (SIPDG) method combined with an explicit time-stepping scheme. These conditions hold for generic meshes, including unstructured non-conforming heterogeneous meshes of mixed element types, and apply to a large class of linear hyperbolic problems, including the acoustic wave equation, the (an)isotropic elastic wave equations and Maxwell's equations. The penalty term bounds are computed element-wise, while bounds for the time step size are computed at weighted submeshes requiring only a small number of elements and faces. Numerical results illustrate the sharpness of these bounds.

Tables

Table 1. Table 1 . Parameter estimates on regular meshes, using penalty term η ∗ superscript 𝜂 \eta^{*} .

Mesh	p	$c_{min}^{\ast }$	$\Delta t\left({\eta }_{min}^{\ast }\right)$	$\Delta t\left({\eta }^{\ast }\right)$	$\Delta t_{est}\left({\eta }^{\ast }\right)$	$\frac{\Delta t\left({\eta }_{min}^{\ast }\right)}{\Delta t\left({\eta }^{\ast }\right)}$	$\frac{\Delta t\left({\eta }^{\ast }\right)}{\Delta t_{est}\left({\eta }^{\ast }\right)}$
1D	1	1.00	0.5774	0.5774	0.5774	1.00	1.00
	2	1.00	0.2582	0.2582	0.2582	1.00	1.00
	3	1.00	0.1533	0.1533	0.1533	1.00	1.00
square	1	0.25	0.4082	0.2357	0.2019	1.73	1.17
	2	0.33	0.1826	0.1170	0.0956	1.56	1.22
	3	0.38	0.1084	0.0694	0.0554	1.56	1.25
triangular	1	0.67	0.2579	0.2273	0.1948	1.13	1.17
	2	0.69	0.1406	0.1250	0.1048	1.12	1.19
	3	0.70	0.0906	0.0739	0.0621	1.23	1.19
cubic	1	0.14	0.3333	0.1361	0.1172	2.45	1.16
	2	0.20	0.1491	0.0678	0.0554	2.20	1.22
	3	0.23	0.0885	0.0405	0.0322	2.19	1.26
tetrahedral	1	0.38	0.1035	0.0635	0.0560	1.63	1.13
	2	0.44	0.0598	0.0384	0.0336	1.56	1.14
	3	0.48	0.0360	0.0243	0.0212	1.48	1.15

Table 2. Table 2 . Parameter estimates on regular meshes, using penalty term η ∗ ∗ superscript 𝜂 absent \eta^{**} .

Mesh	p	$c_{min}^{\ast \ast }$	$\Delta t\left({\eta }_{min}^{\ast \ast }\right)$	$\Delta t\left({\eta }^{\ast \ast }\right)$	$\Delta t_{est}\left({\eta }^{\ast \ast }\right)$	$\frac{\Delta t\left({\eta }_{min}^{\ast \ast }\right)}{\Delta t\left({\eta }^{\ast \ast }\right)}$	$\frac{\Delta t\left({\eta }^{\ast \ast }\right)}{\Delta t_{est}\left({\eta }^{\ast \ast }\right)}$
1D	1	1.00	0.5774	0.5774	0.5774	1.00	1.00
	2	1.00	0.2582	0.2582	0.2582	1.00	1.00
	3	1.00	0.1533	0.1533	0.1533	1.00	1.00
square	1	1.00	0.4082	0.4082	0.4082	1.00	1.00
	2	1.00	0.1826	0.1826	0.1826	1.00	1.00
	3	1.00	0.1084	0.1084	0.1084	1.00	1.00
triangular	1	1.00	0.2582	0.2582	0.2427	1.00	1.06
	2	0.96	0.1406	0.1399	0.1275	1.01	1.10
	3	0.96	0.0906	0.0896	0.0755	1.01	1.19
cubic	1	1.00	0.3333	0.3333	0.3333	1.00	1.00
	2	1.00	0.1491	0.1491	0.1491	1.00	1.00
	3	1.00	0.0885	0.0885	0.0885	1.00	1.00
tetrahedral	1	0.75	0.1040	0.0918	0.0803	1.13	1.14
	2	0.74	0.0599	0.0510	0.0455	1.17	1.15
	3	0.81	0.0359	0.0320	0.0279	1.12	1.15

Table 3. Table 3 . Parameter estimates on a regular tetrahedral mesh, using the penalty term, here denoted by η 𝜂 \eta , derived in [ 16 ] , and using penalty terms η ∗ superscript 𝜂 \eta^{*} and η ∗ ∗ superscript 𝜂 absent \eta^{**} .

Mesh	p	$\Delta t\left(\eta \right)$	$\Delta t\left({\eta }^{\ast }\right)$	$\Delta t\left({\eta }^{\ast \ast }\right)$	$\frac{\Delta t\left({\eta }^{\ast }\right)}{\Delta t\left(\eta \right)}$	$\frac{\Delta t\left({\eta }^{\ast \ast }\right)}{\Delta t\left(\eta \right)}$
triangular	1	0.2280	0.2273	0.2582	1.00	1.13
	2	0.1002	0.1250	0.1399	1.25	1.40
	3	0.0567	0.0739	0.0896	1.30	1.58
tetrahedral	1	0.0689	0.0635	0.0918	0.92	1.33
	2	0.0327	0.0384	0.0510	1.17	1.56
	3	0.0196	0.0243	0.0320	1.24	1.63

Table 4. Table 4 . Parameter estimates on deformed 2D meshes, using penalty term η ∗ ∗ superscript 𝜂 absent \eta^{**} . Intervals denote the range in which the time step ratios are found.

Mesh	x	p	$\frac{\Delta t\left({\eta }_{min}^{\ast \ast }\right)}{\Delta t\left({\eta }^{\ast \ast }\right)}$	$\frac{\Delta t\left({\eta }^{\ast \ast }\right)}{\Delta t_{est}\left({\eta }^{\ast \ast }\right)}$
triangular$\left[x\right]$	$\{0.5,0.7,0.9\}$	1	[1.04, 1.06]	[1.14, 1.20]
		2	[1.05, 1.09]	[1.18, 1.20]
		3	[1.05, 1.09]	[1.19, 1.21]
rectangular$\left[x\right]$	$\{0.5,0.7,0.9\}$	1	[1.00, 1.00]	[1.00, 1.11]
		2	[1.00, 1.00]	[1.00, 1.28]
		3	[1.00, 1.00]	[1.00, 1.37]
quadrilateral$\left[x\right]$	$\{0.5,0.6,0.7\}$	1	[1.00, 1.05]	[1.00, 1.20]
		2	[1.00, 1.04]	[1.00, 1.26]
		3	[1.00, 1.05]	[1.00, 1.31]

Table 5. Table 5 . Parameter estimates on deformed 3D meshes, using penalty term η ∗ ∗ superscript 𝜂 absent \eta^{**} . Intervals denote the range in which the time step ratios are found.

Mesh	x	p	$\frac{\Delta t\left({\eta }_{min}^{\ast \ast }\right)}{\Delta t\left({\eta }^{\ast \ast }\right)}$	$\frac{\Delta t\left({\eta }^{\ast \ast }\right)}{\Delta t_{est}\left({\eta }^{\ast \ast }\right)}$
tetrahedral$\left[x\right]$	$\{0.5,0.7,0.9\}$	1	[1.06, 1.13]	[1.12, 1.14]
		2	[1.08, 1.17]	[1.14, 1.17]
		3	[1.07, 1.12]	[1.14, 1.15]
cuboid$\left[x\right]$	$\{0.5,0.7,0.9\}$	1	[1.00, 1.00]	[1.00, 1.11]
		2	[1.00, 1.00]	[1.00, 1.28]
		3	[1.00, 1.00]	[1.00, 1.37]
hexahedral$\left[x\right]$	$\{0.5,0.6,0.65\}$	1	[1.00, 1.09]	[1.00, 1.17]
		2	[1.00, 1.07]	[1.00, 1.25]
		3	[1.00, 1.10]	[1.00, 1.28]

Table 6. Table 6 . Parameter estimates on 2D and 3D meshes for the acoustic wave equation with piecewise linear parameters ρ 𝜌 \rho and c 𝑐 c , using penalty term η ∗ ∗ superscript 𝜂 absent \eta^{**} . Intervals denote the range in which the time step ratios are found.

Mesh	$({\rho }_0,c_0)$	p	$\frac{\Delta t\left({\eta }_{min}^{\ast \ast }\right)}{\Delta t\left({\eta }^{\ast \ast }\right)}$	$\frac{\Delta t\left({\eta }^{\ast \ast }\right)}{\Delta t_{est}\left({\eta }^{\ast \ast }\right)}$
triPL$[{\rho }_0,c_0]$	$\{(1,1),(5,5),(10,1),\dots$	1	[1.01,1.09]	[1.04,1.17]
	$(1,10),(10,10)\}$	2	[1.01,1.11]	[1.06,1.19]
		3	[1,01,1.09]	[1.08,1.20]
squarePL$[{\rho }_0,c_0]$	$\{(1,1),(5,5),(10,1),\dots$	1	[1.00,1.00]	[1.00,1.00]
	$(1,10),(10,10)\}$	2	[1.00,1.00]	[1.00,1.00]
		3	[1.00,1.00]	[1.00,1.04]
tetraPL$[{\rho }_0,c_0]$	$\{(1,1),(5,5),(10,1),\dots$	1	[1.12,1.19]	[1.13,1.14]
	$(1,10),(10,10)\}$	2	[1.12,1.19]	[1.13,1.15]
		3	[1.10,1.12]	[1.14,1.15]
cubicPL$[{\rho }_0,c_0]$	$\{(1,1),(5,5),(10,1),\dots$	1	[1.00,1.00]	[1.00,1.00]
	$(1,10),(10,10)\}$	2	[1.00,1.00]	[1.00,1.00]
		3	[1.00,1.00]	[1.00,1.03]

Table 7. Table 7 . Parameter estimates on a 3D electromagnetic mesh, using penalty term η ∗ ∗ superscript 𝜂 absent \eta^{**} . Intervals denote the range in which the time step ratios are found.

Mesh	${\mu }_1^{-1}$	${\mu }_2^{-1}$	p	$\frac{\Delta t\left({\eta }_{min}^{\ast \ast }\right)}{\Delta t\left({\eta }^{\ast \ast }\right)}$	$\frac{\Delta t\left({\eta }^{\ast \ast }\right)}{\Delta t_{est}\left({\eta }^{\ast \ast }\right)}$
tetraEM$[{\mu }_1,{\mu }_2]$	1	$\{1,10,100\}$	1	[1.04, 1.05]	[1.14, 1.15]
			2	[1.04, 1.07]	[1.15, 1.15]
			3	[1.00, 1.04]	[1.15, 1.16]
cubicEM$[{\mu }_1,{\mu }_2]$	1	$\{1,10,100\}$	1	[1.00, 1.07]	[1.29, 1.29]
			2	[1.00, 1.00]	[1.31, 1.35]
			3	[1.00, 1.01]	[1.33, 1.35]

Table 8. Table 8 . Parameter estimates on a 3D isotropic elastic tetrahedral mesh, using penalty term η ∗ ∗ superscript 𝜂 absent \eta^{**} . Intervals denote the range in which the time step ratios are found.

Mesh	$(\lambda ,\mu )$	p	$\frac{\Delta t\left({\eta }_{min}^{\ast \ast }\right)}{\Delta t\left({\eta }^{\ast \ast }\right)}$	$\frac{\Delta t\left({\eta }^{\ast \ast }\right)}{\Delta t_{est}\left({\eta }^{\ast \ast }\right)}$
tetraISO$[\lambda ,\mu ]$	$\{(1,0),(0,1),(10,1),(100,1)\}$	1	[1.00, 1.11]	[1.14, 1.15]
		2	[1.00, 1.14]	[1.14, 1.15]
		3	[1.00, 1.13]	[1.13, 1.15]
cubicSO$[\lambda ,\mu ]$	$\{(1,0),(0,1),(10,1),(100,1)\}$	1	[1.05, 1.20]	[1.24, 1.28]
		2	[1.00, 1.10]	[1.23, 1.31]
		3	[1.00, 1.07]	[1.24, 1.33]

Equations

$$[AB{]}_{i_1,\dots ,i_N,j_1,\dots ,j_M}:=$$

$$[A\cdot B{]}_{i_1,\dots ,i_N,j_1,\dots ,j_M}:=$$

$$[A:B{]}_{i_1,\dots ,i_N,j_1,\dots ,j_M}:=$$

$$A_{i_N,\dots ,i_1}^t:=A_{i_1,\dots ,i_N},$$

$$(A,B):=\sum _{i_1=1}^{n_1}\cdots \sum _{i_N=1}^{n_N}{A}_{i_1,\dots ,i_N}{B}_{i_1,\dots ,i_N}.$$

$$\Vert A{\Vert }^2:=(A,A).$$

$$\rho {\partial }_t^{2}u$$

$$\hat{n}\cdot C:\hat{n}u$$

$$\hat{n}\cdot C:\nabla u$$

$$u{|}_{t=0}$$

$${\partial }_{t}u{|}_{t=0}$$

$$0<{\rho }_{min}$$

$$0<c_{min}\le \frac{{\mathbf{\sigma }}^t:C:\mathbf{\sigma }}{\Vert \mathbf{\sigma }{\Vert }^2}$$

$$C:\mathbf{\sigma }$$

$${\partial }_t^{2}p=\nabla \cdot c^2\nabla p,$$

$$C_{i,1,1,j}:={\delta }_{ij}{c}^2$$

$$ϵ{\partial }_t^{2}E=-\nabla \times ({\mu }^{-1}\nabla \times E)+{\partial }_{t}j,$$

$$\displaystyle C_{i_{D},i_{M},j_{M},j_{D}}:=\mu^{-1}\big{(}\delta_{i_{D},j_{D}}\delta_{i_{M},j_{M}}-\delta_{i_{D},j_{M}}\delta_{i_{M},j_{D}}\big{)}$$

$$C_{i_D,i_M,j_M,j_D}=C_{i_M,i_D,j_M,j_D}=C_{i_D,i_M,j_D,j_M}$$

$$C_{i_D,i_M,j_M,j_D}=\lambda {\delta }_{i_D,i_M}{\delta }_{j_D,j_M}+\mu ({\delta }_{i_D,j_D}{\delta }_{i_M,j_M}+{\delta }_{i_D,j_M}{\delta }_{i_M,j_D})$$

$$U$$

$$⟨{\partial }_t(\rho {\partial }_{t}u),w⟩+a(u,w)$$

$$a(u,w):={\int }_{\Omega }(\nabla u{)}^t:C:\nabla wdx.$$

$$\displaystyle U_{e}:=\big{\{}u:e\rightarrow\mathbb{R}^{m}\;|\;u=\tilde{u}\circ\phi_{e}^{-1},\text{ for some }\tilde{u}\in\tilde{U}_{e}\big{\}}.$$

$$\displaystyle U_{h}:=\big{\{}u\in[L^{2}(\Omega)]^{m}\;\big{|}\;u|_{e}\in U_{e},e\in\mathcal{T}_{h}\big{\}}.$$

$$\displaystyle\{\!\!\{\phi\}\!\!\}\big{|}_{f}$$

$$\displaystyle[\![u]\!]\big{|}_{f}$$

$$u^{\ast }{|}_f$$

$${\nu }_h{|}_{\partial e\cap f}$$

$$\rho {|}_e$$

$$(\rho {\partial }_t^{2}u,w)+a_h(u,w)=$$

$$(\rho u{|}_{t=0},w)=$$

$$(\rho {\partial }_{t}u{|}_{t=0},w)=$$

$$a_h(u,w)$$

Full text

*This work was funded by the Shell Global Solutions International B.V. under contract no. PT45999.

1. Introduction

Realistic wave problems often involve large three-dimensional domains, with complex boundary layers, sharp material interfaces and detailed internal structures. Solving such problems therefore requires a numerical method that is efficient in terms of both memory usage and computation time and is flexible enough to deal with interfaces and internal structures without leading to an unnecessary overhead.

A standard finite difference scheme therefore falls short, since it cannot efficiently deal with complex material interfaces, and since small internal structures impose global restrictions on the grid resolution. Finite element methods overcome these problems, since they can be applied to unstructured meshes. However, the finite element method, combined with an explicit time-stepping scheme, requires solving mass matrix-vector equations during every time step. This significantly increases the computational time when the mass matrix is not (block)-diagonal. To obtain diagonal mass matrices without losing the order of accuracy, several mass-lumping techniques have been developed; see, for example, [13, 6, 7, 15]. However, for higher order elements, these techniques require additional quadrature points and degrees of freedom to maintain the optimal order of accuracy.

An alternative is the discontinuous Galerkin (DG) finite element method. This method is similar to the conforming finite element method but allows its approximation functions to be discontinuous at the element boundaries, which naturally results in a block-diagonal mass matrix. Additional advantages of this method are that it also supports meshes with hanging nodes, and that the extension to arbitrary higher order polynomial basis functions is straightforward and can be adapted elementwise. The downside of the DG method, however, is that the discontinuous basis functions can result in significantly more degrees of freedom.

Still, because of its advantages, numerous DG schemes have already been developed and analyzed for linear wave problems, including the symmetric interior penalty discontinuous Galerkin (SIPDG) method; see, for example, [10, 11, 3]. The advantage of the SIPDG method is that it is based on the second order formulation of the problem, while schemes based on a first order formulation require solving additional variables leading to more memory usage. The SIPDG and several alternatives have also been compared and analyzed in [4], from which it follows that the SIPDG method is one of the most attractive options because of its consistency, stability, and optimal convergence rate.

However, to efficiently apply the SIPDG method with an explicit time-stepping scheme, the interior penalty term needs to be sufficiently large and the time step size needs to be sufficiently small. If the penalty term is set too small or the time step size too large, the SIPDG scheme will become unstable. On the other hand, increasing the penalty term will lead to a more severe time step restriction, and a smaller time step size results in a longer computation time. For this reason, there have been multiple studies on finding suitable choices for these parameters.

In [17, 9], for example, sufficient conditions have been derived for the penalty term, for triangular and tetrahedral meshes, while the results of [9] have been sharpened in [16] for tetrahedral meshes. However, the numerical results in Section 7 illustrate that these estimates are still not always very sharp. In [16] they also studied the effects of the penalty term, element shape, and polynomial order on the CFL condition for tetrahedral elements, although these results may not give sufficient conditions for nonuniform grids. Penalty term conditions for regular square and cubic meshes have been studied in [2, 8, 1], where [8, 1] also studied the CFL conditions for these meshes. However, the analysis in these studies only holds for uniform homogeneous meshes. For generic heterogeneous meshes, sharp parameter conditions have remained an open problem.

In this paper we solve these problems by deriving sufficient conditions for both the penalty term and time step size, which lead to sharp estimates, and which hold for generic meshes, including unstructured nonconforming heterogeneous meshes of mixed element types with different types of boundary conditions. These conditions also apply to a large class of linear wave problems, including the acoustic wave equation, Maxwell’s equations, and (an)isotropic elastic wave equations. We compare our estimates to some of the existing ones and illustrate the sharpness of our parameter estimates with several numerical tests.

The paper is constructed as follows: in Section 2 we introduce some tensor notation, such that we can present the general linear hyperbolic model in Section 3, and present the symmetric interior penalty discontinuous Galerkin method in Section 4. After this, we derive sufficient conditions for the penalty parameter in Section 5 and sufficient conditions for the time step size in Section 6. Finally, we compare and test the sharpness of our estimates in Section 7 and give a summary in Section 8.

2. Some tensor notation

Before we present the linear hyperbolic model, it is useful to define some tensor notation. Let $M$ and $N$ be two nonnegative integers. Also, let $A\in\mathbb{R}^{n_{1}\times\cdots\times n_{N}}$ be a tensor of order $N$ and $B\in\mathbb{R}^{m_{1}\times\cdots\times m_{M}}$ be a tensor of order $M$. We define the tensor product $AB\in\mathbb{R}^{n_{1}\times\cdots\times n_{N}\times m_{1}\times\cdots\times m_{M}}$, which is of order $N+M$, as follows:

[TABLE]

for $(i_{1},\dots,i_{N},j_{1},\dots,j_{M})\in(\mathbb{N}_{n_{1}},\dots,\mathbb{N}_{n_{N}},\mathbb{N}_{m_{1}},\dots,\mathbb{N}_{m_{M}})$, where $\mathbb{N}_{n}:=\{1,\dots,n\}$. Now let $A\in\mathbb{R}^{n_{1}\times\cdots\times n_{N}\times p}$ be a tensor of order $N+1$ and $B\in\mathbb{R}^{p\times m_{1}\times\cdots\times m_{M}}$ be a tensor of order $M+1$. We define the dot product $A\cdot B\in\mathbb{R}^{n_{1}\times\cdots\times n_{N}\times m_{1}\times\cdots\times m_{M}}$, which is of order $N+M$, as follows:

[TABLE]

for $(i_{1},\dots,i_{N},j_{1},\dots,j_{M})\in(\mathbb{N}_{n_{1}},\dots,\mathbb{N}_{n_{N}},\mathbb{N}_{m_{1}},\dots,\mathbb{N}_{m_{M}})$. For the double dot product, let $A\in\mathbb{R}^{n_{1}\times\cdots\times n_{N}\times p_{1}\times p_{2}}$ be a tensor of order $N+2$ and $B\in\mathbb{R}^{p_{2}\times p_{1}\times m_{1}\times\cdots\times m_{M}}$ be a tensor of order $M+2$. We define $A:B\in\mathbb{R}^{n_{1}\times\cdots\times n_{N}\times m_{1}\times\cdots\times m_{M}}$, which is of order $N+M$, as follows:

[TABLE]

for $(i_{1},\dots,i_{N},j_{1},\dots,j_{M})\in(\mathbb{N}_{n_{1}},\dots,\mathbb{N}_{n_{N}},\mathbb{N}_{m_{1}},\dots,\mathbb{N}_{m_{M}})$. Now let $A\in\mathbb{R}^{n_{1}\times\cdots\times n_{N}}$ be a tensor of order $N$ again. We define the transpose $A^{t}\in\mathbb{R}^{n_{N}\times\cdots\times n_{1}}$ as follows:

[TABLE]

for $(i_{N},\dots,i_{1})\in(\mathbb{N}_{n_{N}},\dots,\mathbb{N}_{n_{1}})$. A tensor $A$ is called symmetric if $A=A^{t}$ and we define $\mathbb{R}^{n_{1}\times\cdots\times n_{N}}_{sym}$ to be the set of symmetric tensors in $\mathbb{R}^{n_{1}\times\cdots\times n_{N}}$. Now let $B\in\mathbb{R}^{n_{1}\times\cdots\times n_{N}}$ be a tensor of the same size as $A$. We define the inner product as follows:

[TABLE]

The corresponding norm of a tensor is given by

[TABLE]

In the next section we will present the general linear hyperbolic problem, which we will solve using the SIPDG method.

3. A general linear hyperbolic model

Let $\Omega\subset\mathbb{R}^{d}$ be a $d$-dimensional open domain with a Lipschitz boundary $\partial\Omega$, and let $(0,T)$ be the time domain. Also, let $\{\Gamma_{d},\Gamma_{n}\}$ be a partition of $\partial\Omega$, corresponding to Dirichlet and von Neumann boundary conditions, respectively. We define the following linear hyperbolic problem:

[TABLE]

where $u:\Omega\times(0,T)\rightarrow\mathbb{R}^{m}$ is a vector of $m$ variables that are to be solved, $\nabla$ is the gradient operator, $\rho:\Omega\rightarrow\mathbb{R}^{+}$ is a positive scalar field, $C:\Omega\rightarrow\mathbb{R}^{d\times m\times m\times d}_{sym}$ a fourth order tensor field, $f:\Omega\times(0,T)\rightarrow\mathbb{R}^{m}$ the external volume force, and $\hat{n}:\partial\Omega\rightarrow\mathbb{R}^{d}$ the outward pointing normal unit vector.

We make some assumptions on the material parameters. First, we assume that there exist positive constants $\rho_{min},\rho_{max}$ such that

[TABLE]

We also assume that the tensor field is symmetric, $C=C^{t}$, and that there exist linear subspaces $\Sigma_{0},\Sigma_{1}\subset\mathbb{R}^{d\times m}$, with $\Sigma_{0}+\Sigma_{1}=\mathbb{R}^{d\times m}$ and $\Sigma_{0}\perp\Sigma_{1}$, and constants $c_{min},c_{max}$ such that $C$ is nonnegative and bounded in the following sense:

[TABLE]

By $\Sigma_{0}+\Sigma_{1}=\mathbb{R}^{d\times m}$ we mean that any $\boldsymbol{\sigma}\in\mathbb{R}^{d\times m}$ can be written as $\boldsymbol{\sigma}=\boldsymbol{\sigma}_{0}+\boldsymbol{\sigma}_{1}$ for some $\boldsymbol{\sigma}_{0}\in\Sigma_{0},\boldsymbol{\sigma}_{1}\in\Sigma_{1}$, and by $\Sigma_{0}\perp\Sigma_{1}$ we mean that $(\boldsymbol{\sigma}_{0},\boldsymbol{\sigma}_{1})=0$ for all $\boldsymbol{\sigma}_{0}\in\Sigma_{0},\boldsymbol{\sigma}_{1}\in\Sigma_{1}$.

By choosing the correct tensor and scalar field we can obtain a number of hyperbolic problems, including Maxwell’s equations, the acoustic wave equation and the (an)isotropic elastic wave equations. We illustrate this in the following examples, where we define the tensor and vector fields using Cartesian coordinates.

Example 3.1**.**

Consider the acoustic wave equation written in the following dimensionless form:

[TABLE]

where $p:\Omega\times(0,T)\rightarrow\mathbb{R}$ is the pressure field and $c:\Omega\rightarrow\mathbb{R}^{+}$ the acoustic wave velocity. We can write these equations in the form of (1) by setting $m=1$, $u=p$, $f=0$, $\rho=1$, and

[TABLE]

for $i,j=1,\dots,d$, where $\delta$ is the Kronecker delta function.

Example 3.2**.**

Consider Maxwell’s equations in a domain with zero conductivity written in the following dimensionless form:

[TABLE]

where $E:\Omega\times(0,T)\rightarrow\mathbb{R}^{3}$ is the electric field, $\epsilon:\Omega\rightarrow\mathbb{R}^{+}$ the relative electric permittivity, $\mu:\Omega\rightarrow\mathbb{R}^{+}$ the relative magnetic permeability, and $j:\Omega\times(0,T)\rightarrow\mathbb{R}^{3}$ the applied current density. We can write these equations in the form of (1) by setting $d=3$, $m=3$, $u=E$, $f=\partial_{t}j$, $\rho=\epsilon$, and

[TABLE]

for $i_{D},i_{M},j_{D},j_{M}=1,2,3$, where $\delta$ is the Kronecker delta function.

Example 3.3**.**

Consider the linear anisotropic elastic wave equations. They can immediately be written in the form of (1) with $u:\Omega\times(0,T)\rightarrow\mathbb{R}^{d}$ being the displacement vector and $C:\Omega\rightarrow\mathbb{R}^{d\times d\times d\times d}_{sym}$ being the elasticity tensor, which has the additional symmetries

[TABLE]

for $i_{D},i_{M},j_{D},j_{M}=1,\dots,d$. For the special case of isotropic elasticity, we can write

[TABLE]

for $i_{D},i_{M},j_{D},j_{M}=1,\dots,d$, where $\delta$ is the Kronecker delta function and $\lambda,\mu$ are the Lamé parameters.

In the next section we present the DG method that we use to solve these linear hyperbolic problems.

4. A discontinuous Galerkin method

To explain the DG method, we first present the weak formulation of (1). After that, we introduce the tesselation of the domain, the discrete function spaces, and the trace operators. We then present the symmetric interior penalty discontinuous Galerkin (SIPDG) method and derive some of its properties.

4.1. The weak formulation

Define the following function space:

[TABLE]

Assume that $u_{0}\in U$, $v_{0}\in L^{2}(\Omega)^{m}$, and f\in L^{2}\big{(}0,T;L^{2}(\Omega)^{m}\big{)}. The weak formulation of (1) is finding u\in L^{2}\big{(}0,T;U\big{)}, with \partial_{t}u\in L^{2}\big{(}0,T;L^{2}(\Omega)^{m}\big{)} and \partial_{t}(\rho\partial_{t}u)\in L^{2}\big{(}0,T;U^{-1}\big{)}, such that $u|_{t=0}=u_{0}$, $\partial_{t}u|_{t=0}=v_{0}$, and

[TABLE]

Here $(\cdot,\cdot)$ denotes the inner product of $L^{2}(\Omega)^{m}$, $\langle\cdot,\cdot\rangle$ denotes the pairing between $U^{-1}$ and $U$, and $a(\cdot,\cdot):U\times U\rightarrow\mathbb{R}$ is the semielliptic bilinear operator given by

[TABLE]

Using (3) it can be shown that $U$ is a separable Hilbert space, and from (2) it follows that the norm $\|u\|_{\rho}^{2}:=(\rho u,u)$ is equivalent to the standard $L^{2}(\Omega)^{m}$ inner product. Using these properties, it can be proven, in a way analogous to the proof of [14, Chapter 3, Theorem 8.1] that (4) is well-posed and has a unique solution.

4.2. Tesselation, discrete function space, and trace operators

Let $\mathcal{T}_{h}$ be a set of nonoverlapping open domains in $\mathbb{R}^{d}$, referred to as elements, such that every element $e\in\mathcal{T}_{h}$ fits inside a $d$-dimensional sphere of radius $h$, and such that $\overline{\Omega}:=\bigcup_{e\in\mathcal{T}_{h}}\overline{e}$, where $\overline{e}$ and $\overline{\Omega}$ are the closures of $e$ and $\Omega$, respectively. We call $\mathcal{T}_{h}$ the tesselation of $\Omega$. Using the tesselation we define the set of faces $\mathcal{F}_{h}:=\mathcal{F}_{h,in}\cup\mathcal{F}_{h,b}$ and the union of all faces $\Gamma_{h}:=\bigcup_{f\in\mathcal{F}_{h}}f$, where \mathcal{F}_{h,in}:=\big{\{}\partial e\cap\partial e^{\prime}\big{\}}_{e,e^{\prime}\in\mathcal{T}_{h}} is the set of all internal faces, \mathcal{F}_{h,b}:=\big{\{}\partial e\cap\partial\Omega\big{\}}_{e\in\mathcal{T}_{h}} is the set of all boundary faces, and $\partial e$ denotes the element boundary. Furthermore, we let $\{\mathcal{F}_{h,d},\mathcal{F}_{h,n}\}$ be the partition of $\mathcal{F}_{h,b}$ corresponding to the Dirichlet and Neumann boundary conditions, such that $\bigcup_{f\in\mathcal{F}_{h,d}}f=\Gamma_{d}$ and $\bigcup_{f\in\mathcal{F}_{h,n}}f=\Gamma_{n}$.

We use these sets of elements and faces to construct the discrete function space. To do this, let $e\in\mathcal{T}_{h}$ be an element with $\tilde{e}$ the corresponding reference element, which depends only on the shape of $e$. For every reference element $\tilde{e}$ we define a discrete function space $\tilde{U}_{e}:\tilde{e}\rightarrow\mathbb{R}^{m}$, such that $\tilde{U}_{e}=[\mathcal{P}^{K_{e}}(\tilde{e})]^{m}$, where $\mathcal{P}^{K_{e}}(\tilde{e})$ is a finite set of polynomial functions on $\tilde{e}$, which depends only on the degree $K_{e}$ and $\tilde{e}$. Now let $\phi_{e}:\tilde{e}\rightarrow e$ be an invertible polynomial mapping, such that $|\mathbf{J}_{e}|:=|\mathrm{det}(\nabla\phi_{e})|\geq|\mathbf{J}_{e}|_{min}>0$ in $e$, for some positive constant $|\mathbf{J}_{e}|_{min}$. Using this mapping and the reference function space $\tilde{U}_{e}$, we can construct the function space on the physical element $U_{e}$ as follows:

[TABLE]

This can then be used to construct the discrete finite element space:

[TABLE]

The functions in the finite element space are continuous within every element, but can be discontinuous at the faces between two elements. To construct a discrete version of the bilinear form $a$, which can deal with these discontinuities, we introduce trace operators $\{\!\!\{\cdot\}\!\!\}$ and $[\![\cdot]\!]$ given below:

[TABLE]

where $\mathcal{T}_{f}$ is the set of elements adjacent to $f$, $|\mathcal{T}_{f}|$ is the number of elements adjacent to $f$, and $\hat{n}|_{\partial e}$ is the outward pointing normal vector of element $e$. The first trace operator is the average of traces, while the second operator is known as the jump operator. Using the first trace operator $\{\!\!\{\cdot\}\!\!\}$, we can construct the numerical flux operator $(\cdot)^{*}:U_{h}\rightarrow[L^{2}(\Gamma_{h})]^{m}$, which assigns a unique value for $u\in U_{h}$ at the faces, as follows:

[TABLE]

In order to ensure that the discrete bilinear form $a_{h}$ remains semielliptic, we also introduce penalty terms $\eta_{e}\in\mathbb{R}^{+}$ for every element $e\in\mathcal{T}_{h}$, and a penalty scaling function $\nu_{h}\in\bigotimes_{e\in\mathcal{T}_{h}}L^{\infty}(\partial e)$, with $\nu_{h}>0$, which means $\nu|_{\partial e}:\partial e\rightarrow\mathbb{R}^{+}$ for all $e\in\mathcal{T}_{h}$. The penalty terms $\eta_{e}$ are positive dimensionless constants for which lower bounds will be derived in Section 5. The function $\nu_{h}$ scales with order $h^{-1}$ and is chosen as follows:

[TABLE]

where $\mathcal{F}_{e}$ denotes the faces adjacent to element $e$, $|\mathbf{J}_{e}|:=|\mathrm{det}(\nabla\phi_{e})|$ is the reference-to-physical element scale, and $|\mathbf{J}_{f}|$ is the reference-to-physical face scale. The face scale satisfies $|\mathbf{J}_{f}|=1$ in 1D, $|\mathbf{J}_{f}|=|\partial_{1}\phi_{f}|$ in 2D, and $|\mathbf{J}_{f}|=|\partial_{1}\phi_{f}\times\partial_{2}\phi_{f}|$ in 3D, where $\phi_{f}:\tilde{f}\rightarrow f$ is the reference-to-physical face mapping, and $\partial_{i}\phi_{f}$ is the derivative of $\phi_{f}$ in reference coordinate $i$, assuming Cartesian reference coordinates. In our numerical tests we use this scaling function, although the stability analysis in this paper holds for arbitrary positive functions $\nu_{h}$.

Finally, to ensure that the discrete version of $a$ is well defined we also make the following additional assumptions on the material parameters $\rho$ and $C$:

[TABLE]

where $W^{1,\infty}$ denotes the Sobolev space of differentiable functions with uniformly bounded weak derivatives. These assumptions together with the trace inequality imply that the element traces of $C$ and $\rho$ are well defined and bounded.

We have now introduced the function spaces, operators, and parameter assumptions needed to present the DG method in the next subsection.

4.3. The symmetric interior penalty discontinuous Galerkin method

We present a DG method, which is known as the symmetric interior penalty discontinuous Galkerkin (SIPDG) method. The SIPDG method is solving $u:[0,T]\rightarrow U_{h}$ such that

[TABLE]

where $a_{h}:U_{h}\times U_{h}\rightarrow\mathbb{R}$ is the discrete version of the elliptic operator, given by

[TABLE]

with

[TABLE]

for all $u,w\in U_{h}$. The bilinear form $a^{(C)}_{h}$ is the same as the original elliptic operator $a$ and is the part that remains when both input functions are continuous. The bilinear form $a^{(DG)}_{h}$ can be interpreted as the additional part that results from partial integration of the elliptic operator $a$ when the first input function is discontinuous. Finally, the bilinear form $a^{(IP)}$ is the part that contains the interior penalty terms needed to ensure stability of the scheme.

Using the definition of the numerical flux in (5), we can rewrite the bilinear forms $a_{h}^{(DG)}$ and $a_{h}^{(IP)}$, as follows:

[TABLE]

for all $u,w\in U_{h}$. Here $\epsilon_{f}:=1/2$ for internal faces and $\epsilon_{f}:=1$ for faces in $\mathcal{F}_{h,d}$, and $\eta_{h}\in\bigotimes_{e\in\mathcal{T}_{h}}L^{\infty}(\partial e)$ is defined by $\eta_{h}|_{\partial e}:=\eta_{e}$ for all $e\in\mathcal{T}_{h}$. This scheme conforms with existing SIPDG schemes, except for a possible deviation in the interior penalty part $a_{h}^{(IP)}$. For example, for the acoustic wave equation given in Example 3.1, this scheme is equivalent to the one in [10] when choosing their penalty term $\mathbf{a}$ as $\mathbf{a}|_{f}=\epsilon_{f}\{\!\!\{\eta\nu_{h}c\}\!\!\}|_{f}$ for all $f\in\mathcal{F}_{h}$.

Since the bilinear form is symmetric, $a_{h}(u,w)=a_{h}(w,u)$ for all $u,w\in U_{h}$, we can obtain, by substituting $w=\partial_{t}u$ into (6a), the following energy equation:

[TABLE]

where E_{h}:=\frac{1}{2}\big{\|}\rho^{1/2}w\big{\|}_{0}^{2}+\frac{1}{2}a_{h}(u,u) is the discrete energy, with $\|\cdot\|_{0}$ the $[L^{2}(\Omega)]^{m}$ norm. In the absence of an external force $f$ this implies that the discrete energy is conserved.

However, for this energy to be well defined, in the sense that it is always nonnegative, the discrete bilinear form needs to remain semielliptic: $a_{h}(u,u)\geq 0$ for all $u\in U_{h}$. This then implies that any nonzero discrete eigenmode cannot grow unbounded in the absence of an external force. In case $u$ is a zero discrete eigenmode, we can substitute $w=u$ into (6a) to obtain $\partial_{t}^{2}\|\rho^{1/2}u\|^{2}_{0}=2(f,u)$. This implies that, in the absence of an external force, zero eigenmodes grow at most linearly in time. This behavior can correspond to physical rigid motions, or, when there is a discrepancy between the physical and discrete zero modes, a linear drift of a spurious mode. For acoustic and elastic waves these spurious modes are absent, while for electromagnetic waves, these modes have been analyzed in, for example, [5]. However, even if there are spurious modes, we will not consider their drift as numerical instability, since the numerical error is expected to grow linearly in time anyway due to dispersion errors.

In the next section we will find sufficient lower bounds for the penalty term to make sure $a_{h}$ is semielliptic. In particular, we will show there that $a_{h}$ satisfies a coercivity condition that is commonly used to show optimal convergence in the energy-norm.

5. Sufficient penalty term estimates

In this section we derive a sufficient lower bound for the penalty term and a positive constant $c_{coer}>0$, where $c_{coer}$ is independent of the mesh $\mathcal{T}_{h}$, such that

[TABLE]

where $|\cdot|_{1,h}$ is the discrete seminorm defined by \big{|}u\big{|}^{2}_{1,h}:=\sum_{e\in\mathcal{T}_{h}}\big{|}u\big{|}^{2}_{1,e}, with

[TABLE]

Here $C^{1/2}\in\bigotimes_{e\in\mathcal{T}_{h}}W^{1,\infty}(e)^{d\times m\times m\times d}_{sym}$ is a tensor field such that $C^{1/2}:C^{1/2}=C$. The existence of such a tensor field follows from Lemma A.1. The numerical flux $u^{*}$ is defined in (5), although the stability analysis in this paper holds for arbitrary linear flux operators. Note that (9) is satisfied when

[TABLE]

where $a_{e}(u,w):=a_{e}^{(C)}(u,w)+a_{\partial e}^{(DG)}(u,w)+a_{\partial e}^{(DG)}(w,u)+\eta_{e}a_{\partial e}^{(IP)}(u,w)$. Since we can write $|u|_{1,e}^{2}=a_{e}^{(C)}(u,u)+\eta_{e}a_{\partial e}^{(IP)}(u,u)$, it remains to bound $a_{\partial e}^{(DG)}(u,u)$ in terms of $a_{e}^{(C)}(u,u)$ and $a_{\partial e}^{(IP)}(u,u)$. In order to do this, we first introduce the auxiliary bilinear form $a_{\partial e}^{(C*)}:U_{h}\times U_{h}\rightarrow\mathbb{R}$ defined by

[TABLE]

Note that this operator is similar to $a_{e}^{(C)}$, but integrates over the element boundary instead of the interior. Next, we show that $a_{\partial e}^{(DG)}(u,u)$ can be bounded in terms of $a_{\partial e}^{(C*)}(u,u)$ and $a_{\partial e}^{(IP)}(u,u)$:

Lemma 5.1**.**

Consider an arbitrary element $e\in\mathcal{T}_{h}$, and let $c>0$ be an arbitrary positive constant. Then the following inequality holds:

[TABLE]

Proof.

Take an arbitrary function $u\in U_{h}$. We can write

[TABLE]

Using the Cauchy–Schwarz and the Cauchy inequalities, we can then obtain

[TABLE]

∎

We now construct the following constant:

[TABLE]

where $U_{e}$ is the discrete function space restricted to element $e$. From its definition, we can immediately obtain the inequality $a_{\partial e}^{(C*)}(u,u)\leq\kappa_{e}^{*}a_{e}^{(C)}(u,u)$ for any $u\in U_{h}$. Using this property and Lemma 5.1 we can prove in Theorem 5.6 that $\kappa_{e}^{*}$ is a sufficient lower bound for $\eta_{e}$ to ensure $a_{e}$ to be coercive.

However, before we give this proof, we first show that $\kappa_{e}^{*}$ is well defined and show how it can be computed efficiently. To do this, consider an arbitrary $e\in\mathcal{T}_{h}$, and let $\{w_{i}\}_{i=1}^{n}$ be a set of basis functions spanning $U_{e}$. Using these basis functions we can define positive semidefinite matrices $A_{e},A^{*}_{\partial e}\in\mathbb{R}^{n\times n}_{sym}$ as follows:

[TABLE]

If matrix $A$ had been positive definite, then we could have obtained, using Lemma A.5, the following relation:

[TABLE]

where $\lambda_{max}(A_{e}^{-1}A^{*}_{\partial e})$ denotes the largest eigenvalue of $A_{e}^{-1}A^{*}_{\partial e}$. However, the matrix $A_{e}$ is only positive semidefinite, so we need to obtain some intermediate results before we can show that a similar type of relation still holds. First, we show that the kernel of $A_{e}$ is a subset of the kernel of $A^{*}_{\partial e}$.

Lemma 5.2**.**

Let $e\in\mathcal{T}_{h}$ be an arbitrary element, and let $A_{e},A^{*}_{\partial e}$ be matrices defined as in (12). Then $\mathrm{Ker}(A_{e})\subset\mathrm{Ker}(A^{*}_{\partial e})$.

Proof.

Let $\underline{u}\in\mathrm{Ker}(A_{e})$, and define $u\in U_{e}$ as follows: $u:=\sum_{i=1}^{n}\underline{u}_{i}w_{i}$. Then

[TABLE]

From this it follows that $C:\nabla u=\mathbf{0}$ in $e$. Now let $\tilde{u}:=u\circ\phi_{e}$ be the reference function, with $\phi_{e}$ the reference-to-physical element mapping. Since $\phi_{e}$ is assumed to be a polynomial function, such that $|\mathbf{J}_{e}|:=|\mathrm{det}(\nabla\phi_{e})|\geq|\mathbf{J}_{e}|_{min}>0$ in $e$, for some constant $|\mathbf{J}_{e}|_{min}$, it follows that $\phi_{e}^{-1}\in C^{\infty}(e)^{d}$. Furthermore, since the reference function $\tilde{u}$ is also assumed to be polynomial, this implies that $u=\tilde{u}\circ\phi_{e}^{-1}\in C^{\infty}(e)^{d}$. Because $C|_{e}\in W^{1,\infty}(e)^{d\times m\times m\times d}_{sym}$, it then follows from the trace theorem that $C:\nabla u=\mathbf{0}$ is also satisfied on the boundary $\partial e$. This implies $a^{(C*)}_{\partial e}(u,w)=0$ for any $w\in U_{e}$ and therefore $\underline{u}\in\mathrm{Ker}(A^{*}_{\partial e})$. ∎

Now let $k$ be the rank of $A_{e}$, and let $V_{e}\in\mathbb{R}^{n\times n}$ be a nonsingular matrix such that $V_{e}^{t}A_{e}V_{e}=D_{e}$, where $D_{e}\in\mathbb{R}^{n\times n}_{sym}$ is a diagonal matrix with the last $n-k$ entries being zero. Such a matrix decomposition can be obtained from, for example, a symmetric Gauss elimination procedure or a singular value decomposition. We then use these matrices to construct matrices $\tilde{D}_{e},\tilde{A}^{*}_{\partial e}\in\mathbb{R}^{k\times k}_{sym}$ as follows:

[TABLE]

Using these matrices $\tilde{D}_{e}$ and $\tilde{A}^{*}_{\partial e}$ we can compute $\kappa_{e}^{*}$ and show that it is well defined.

Lemma 5.3**.**

Let $e\in\mathcal{T}_{h}$ be an arbitrary element, and let $\tilde{D}_{e},\tilde{A}^{*}_{\partial e}$ be the matrices defined as in (13). The constant $\kappa_{e}^{*}$ is well defined and satisfies

[TABLE]

where \lambda_{max}\big{(}\tilde{D}_{e}^{-1}\tilde{A}^{*}_{\partial e}\big{)} denotes the largest eigenvalue of $\tilde{D}_{e}^{-1}\tilde{A}^{*}_{\partial e}$.

Proof.

First, consider the decomposition $V_{e}^{t}A_{e}V_{e}=D_{e}$ which was used to construct $\tilde{D}_{e}$ and $\tilde{A}^{*}_{\partial e}$. Since matrix $A_{e}$ has rank $k$ and the last $n-k$ entries of $D_{e}$ are zero, and since $V_{e}$ is nonsingular, this implies that the last $n-k$ columns of $V_{e}$ span the kernel of $A_{e}$. From Lemma 5.2 it follows that these columns are also in the kernel of $A^{*}_{\partial e}$. Now let $\underline{w}\in\mathbb{R}^{n}$, and let $\tilde{}\underline{w}\in\mathbb{R}^{k}$ be the vector composed of the first $k$ entries of $\underline{w}$. We can then obtain the following relation:

[TABLE]

Since $A_{e}$ is positive semidefinite, it also follows that all entries of $\tilde{D}_{e}$ are strictly positive. Furthermore, since $A^{*}_{\partial e}$ is positive semidefinite, the matrix $\tilde{A}^{*}_{\partial e}$ will be positive semidefinite as well. Using these properties, we can prove (14) as follows:

[TABLE]

In the third step we substituted $\underline{u}$ by $V_{e}\underline{w}$, in the fourth step we used (15), and in the last step we used Lemma A.5 combined with the fact that $\tilde{D}_{e}$ is positive definite and $\tilde{A}^{*}_{\partial e}$ is positive semidefinite. ∎

Remark 5.4**.**

A symmetric Gauss elimination procedure or a singular value decomposition algorithm usually does not give the exact decomposition $V_{e}^{t}A_{e}V_{e}=D_{e}$, but only a numerical approximation. The diagonal entries of $D_{e}$ are then considered to be [math] when they are smaller than a given tolerance.

Remark 5.5**.**

The largest eigenvalue \lambda_{max}\big{(}\tilde{D}_{e}^{-1}\tilde{A}^{*}_{\partial e}\big{)} can be efficiently obtained using a power iteration method.

We can now derive the following sufficient estimate for the penalty term.

Theorem 5.6**.**

Let $e\in\mathcal{T}_{h}$ be an arbitrary element, and let $c_{\kappa}\geq 1$ be an arbitrary constant. If $\eta_{e}\geq c_{\kappa}\kappa^{*}_{e}$, then $a_{e}(u,u)\geq 0$ for all $u\in U_{h}$. Moreover, if $c_{\kappa}>1$, then

[TABLE]

where

[TABLE]

Proof.

Take an arbitrary function $u\in U_{h}$ and scalar $x\in[1,c_{\kappa}]$. We can then derive the following inequality:

[TABLE]

In the second line we used Lemma 5.1 with $c=x\kappa^{*}_{e}$, and in the last line we used the definition of $\kappa^{*}_{e}$. Now note that we can write $|u|^{2}_{1,e}=a^{(C)}_{e}(u,u)+c_{\kappa}\kappa^{*}_{e}a^{(DG)}_{\partial e}(u,u)$. Combining this with the inequality above gives

[TABLE]

Taking the supremum over all $x\in[1,c_{\kappa}]$ results in (16). ∎

The penalty term estimate depends on the constant $\kappa^{*}_{e}$. However, this constant does not include any effects of the normal vector on the positivity of the bilinear operator, which may cause the penalty term estimate to be less sharp. Therefore, we consider an additional penalty term estimate which does include the effect of the normal vector, and is shown to be considerably sharper in Section 7. To do this, we first define the tensor field $\mathbf{c}_{\hat{n}}\in\bigotimes_{e\in\mathcal{T}_{h}}L^{\infty}(\partial e)^{m\times m}_{sym}$ as follows:

[TABLE]

where $\hat{n}|_{\partial e}$ is the outward pointing normal vector of element $e$. We also define the following function space:

[TABLE]

Lemma A.2 shows that there exists a pseudoinverse $\mathbf{c}_{\hat{n}}^{-1}\in\bigotimes_{e\in\mathcal{T}_{h}}L^{\infty}(\partial e)^{m\times m}_{sym}$ such that $\mathbf{c}_{\hat{n}}^{-1}\cdot\mathbf{c}_{\hat{n}}\cdot\hat{u}=\mathbf{c}_{\hat{n}}\cdot\mathbf{c}_{\hat{n}}^{-1}\cdot\hat{u}=\hat{u}$ for all $\hat{u}\in\hat{U}_{h}$. We use this tensor field to define an alternative auxiliary bilinear operator $a^{(C**)}_{\partial e}:U_{h}\times U_{h}\rightarrow\mathbb{R}$ as follows:

[TABLE]

The penalty term estimate and the coercivity result are obtained in the same way as before, except that we now use $a^{(C**)}_{\partial e}$ instead of $a^{(C*)}_{\partial e}$. We start again by deriving a bound on $a^{(DG)}_{\partial e}$:

Lemma 5.7**.**

Consider an arbitrary element $e\in\mathcal{T}_{h}$, and let $c>0$ be an arbitrary positive constant. Then the following inequality holds:

[TABLE]

Proof.

From Lemma A.2 it follows that $\mathbf{c}_{\hat{n}}$ and $\mathbf{c}_{\hat{n}}^{-1}$ are positive semidefinite tensor fields, and therefore there exist symmetric positive semidefinite tensor fields $\mathbf{c}_{\hat{n}}^{1/2},\mathbf{c}_{\hat{n}}^{-1/2}\in\bigotimes_{e\in\mathcal{T}_{h}}L^{\infty}(\partial e)^{m\times m}_{sym}$ such that $\mathbf{c}_{\hat{n}}^{1/2}\cdot\mathbf{c}_{\hat{n}}^{1/2}=\mathbf{c}_{\hat{n}}$ and $\mathbf{c}_{\hat{n}}^{-1/2}\cdot\mathbf{c}_{\hat{n}}^{-1/2}=\mathbf{c}_{\hat{n}}^{-1}$, and such that $\mathbf{c}_{\hat{n}}^{-1/2}\cdot\mathbf{c}_{\hat{n}}^{1/2}\cdot\hat{u}=\mathbf{c}_{\hat{n}}^{1/2}\cdot\mathbf{c}_{\hat{n}}^{-1/2}\cdot\hat{u}=\hat{u}$ for all $\hat{u}\in\hat{U}_{h}$.

Now take an arbitrary function $u\in U_{h}$, and define the function $\hat{u}\in\hat{U}_{h}$ as follows:

[TABLE]

We can then write

[TABLE]

Using the Cauchy–Schwarz and the Cauchy inequalities, we can then obtain

[TABLE]

∎

We now use the bilinear operator $a_{\partial e}^{(C**)}$ to construct the following constant:

[TABLE]

The proof of the existence of this constant and the way to compute it is analogous to $\kappa_{e}^{*}$. In a similar way as before we can use this constant to obtain the following sufficient penalty term estimate.

Theorem 5.8**.**

Let $e\in\mathcal{T}_{h}$ be an arbitrary element, and let $c_{\kappa}\geq 1$ be an arbitrary constant. If $\eta_{e}\geq c_{\kappa}\kappa^{**}_{e}$, then $a_{e}(u,u)\geq 0$ for all $u\in U_{h}$. Moreover, if $c_{\kappa}>1$, then

[TABLE]

where

[TABLE]

Proof.

The proof is analogous to that of Theorem 5.6. ∎

We have now derived conditions for the penalty term to ensure that $a_{h}$ is positive semidefinite and showed how the penalty term can be computed. In the next section we will derive time step estimates to ensure that the local time-stepping scheme is stable.

6. Sufficient time step estimates

We start by rewriting the DG method as a linear system of ordinary differential equations. We then show how we can obtain sufficient upper bounds for the spectral radius of $M^{-1}A$, and therefore sufficient lower bounds for the time step size for a large class of explicit time integration schemes, by splitting the mass mass matrix $M$ and stiffness matrix $A$ into multiple parts. Finally, we introduce weighted mesh decompositions to explain how this splitting of matrices can be done efficiently.

6.1. A system of ordinary differential equations

Let $\{w_{i}\}_{i=1}^{N}$ be a linear basis of $U_{h}$ and define, for $u\in U_{h}$, the vector $\underline{u}\in\mathbb{R}^{N}$ such that $u=\sum_{i=1}^{N}\underline{u}_{i}w_{i}$. We can rewrite the DG method, given in (6), as the following system of ordinary differential equations: we solve $\underline{u}:[0,T]\rightarrow\mathbb{R}^{N}$, such that

[TABLE]

where $M_{h},A_{h}\in\mathbb{R}^{N\times N}$ are matrices, $\underline{u}_{0,h}^{*},\underline{v}_{0,h}^{*}\in\mathbb{R}^{N}$ are vectors, and $\underline{f}_{h}^{*}:[0,T]\rightarrow\mathbb{R}^{N}$ is a vector function, defined as follows:

[TABLE]

For a large class of explicit time integrators, including Lax–Wendroff schemes and explicit Runge–Kutta schemes, the time step size condition is of the form

[TABLE]

where $c_{method}>0$ is a constant, depending only on the type of time integration method, and $\lambda_{max}(M^{-1}A)$ is the largest eigenvalue of $M^{-1}A$, which is also known as the spectral radius of $M^{-1}A$. For example, the stability condition for the leap–frog scheme is well known to be (20) with $c_{method}=2$. Because of the form of (20), it remains to find an upper estimate for the spectral radius. In the next section we show how this can be done by splitting the matrices $M$ and $A$ into multiple parts.

6.2. Spectral radius estimates by splitting matrices

In order to obtain a bound for the spectral radius we first introduce the mapping $\mathcal{I}:\mathbb{R}^{N\times N}_{sym}\rightarrow\mathbb{R}^{N\times N}_{sym}$, which maps a symmetric matrix to a diagonal matrix with entries [math] and $1$ to indicate the nonzero rows or columns of the input matrix:

[TABLE]

We also define $\mathcal{I}^{*}(S)$ as the matrix $\mathcal{I}(S)$ with all zero-columns removed. Using these definitions we can formulate the following theorem.

Theorem 6.1**.**

Let $M\in\mathbb{R}^{N\times N}_{sym}$ be a symmetric positive definite matrix and $A\in\mathbb{R}^{N\times N}_{sym}$ a symmetric positive semidefinite matrix. Also let $M_{(i)},A_{(i)}\in\mathbb{R}^{N\times N}_{sym}$ for $i=1,\dots,n$, be symmetric matrices such that

[TABLE]

where $M_{(i)}^{\prime}:=(I_{(i)}^{*})^{t}M_{(i)}I^{*}_{(i)}$, $I_{(i)}:=\mathcal{I}(M_{(i)})$ and $I_{(i)}^{*}:=\mathcal{I}^{*}(M_{(i)})$. Then

[TABLE]

where $A_{(i)}^{\prime}:=(I^{*}_{(i)})^{t}A_{(i)}I^{*}_{(i)}$, and $\lambda_{max}(\cdot)$ denotes the largest eigenvalue in magnitude.

Remark 6.2**.**

The matrices $M^{\prime}_{(i)}$ and $A^{\prime}_{(i)}$ are the submatrices of $M_{(i)}$ and $A_{(i)}$, respectively, obtained by removing all rows and columns corresponding to the zero rows and columns of $M_{(i)}$. The condition $I_{(i)}A_{(i)}I_{(i)}=A_{(i)}$ means that any zero column or row of $M_{(i)}$ is also a zero column or row of $A_{(i)}$, and the condition $M_{(i)}^{\prime}\succ 0$ means that the submatrices of $M$ are positive definite.

Proof.

For any $\underline{u}\in\mathbb{R}^{n}$, define the following set of indices:

[TABLE]

Using Lemma A.5 and Lemma A.3 we can then bound the largest eigenvalue as follows:

[TABLE]

Using Lemma A.5 again, we can obtain, for any $i=1,..,n$, the following:

[TABLE]

where $N_{i}$ is the number of nonzero columns of $M_{(i)}$. Combining these results gives (22). ∎

To apply the above theorem it remains to find a decomposition of the matrices $M$ and $A$ such that (21) is satisfied. For continuous finite elements such a decomposition can be easily obtained from the element matrices,

[TABLE]

where $M_{e}$ and $A_{e}$ are the element matrices corresponding to the mass matrix $M$ and stiffness matrix $A$, respectively. Using Theorem 6.1 we then obtain the following estimate for the spectral radius:

[TABLE]

where $M_{e}^{\prime}:=(I_{e}^{*})^{t}M_{e}I_{e}^{*}$, $A_{e}^{\prime}:=(I_{e}^{*})^{t}A_{e}I_{e}^{*}$ and $I_{e}^{*}:=\mathcal{I}^{*}(M_{e})$. In other words, the largest eigenvalue of the global matrix is bounded by the supremum over all elements of the largest eigenvalue of the element matrix. This result was already mentioned by [12]. For discontinuous elements, however, a suitable decomposition of the matrices is less straightforward due to the face integral terms. In the next subsection we show how we can decompose the matrices for discontinuous elements, using a weighted mesh decomposition.

6.3. A weighted mesh decomposition

We define a weighted submesh $\omega:\mathcal{T}_{h}\cup\mathcal{F}_{h}\rightarrow[0,1]$ to be a function that assigns to every element and face a weight value $\omega_{e}$ and $\omega_{f}$ between [math] and $1$, such that if $\omega_{f}>0$ for a certain face $f$, then $\omega_{e}>0$ for the adjacent elements $e\in\mathcal{T}_{f}$. We call a set of weighted submeshes $\mathcal{W}_{h}$ a weighted mesh decomposition of $\mathcal{T}_{h}$ if the sum of all weighted submeshes adds up to one for every face and element: $\sum_{\omega\in\mathcal{W}_{h}}\omega_{e}=1$ for all $e\in\mathcal{T}_{h}$ and $\sum_{\omega\in\mathcal{W}_{h}}\omega_{f}=1$ for all $f\in\mathcal{F}_{h}$. An illustration of a weighted mesh decomposition is given in Figure 1.

We can use a weighted submesh to construct bilinear forms $(\cdot,\cdot)_{\omega},a_{\omega}:U_{h}\times U_{h}\rightarrow\mathbb{R}$ as follows:

[TABLE]

with

[TABLE]

Note that $(u,w)=\sum_{\omega\in\mathcal{W}_{h}}(u,w)_{\omega}$ and $a_{h}(u,w)=\sum_{\omega\in\mathcal{W}_{h}}a_{\omega}(u,w)$, for all $u,w\in U_{h}$.

For the numerical tests, we will in particular consider a weighted mesh decomposition based on the vertices, as illustrated in Figure 2. The vertex-based mesh decomposition is given by $\mathcal{W}_{h}:=\{\omega^{(q)}\}_{q\in\mathcal{Q}}$, with

[TABLE]

where $|\mathcal{Q}_{e}|$, $|\mathcal{Q}_{f}|$ are the number of vertices adjacent to element $e$ and face $f$, respectively, and $\mathcal{T}_{q}$, $\mathcal{F}_{q}$ are the set of elements and faces adjacent to $q$, respectively.

Now let $\{w_{i}\}_{i=1}^{N}$ be a linear basis of $U_{h}$, such that every basis function is nonzero on only a single element $e_{i}$. We can use a weighted mesh decomposition $\mathcal{W}_{h}$ to decompose the mass matrix and stiffness matrix as follows:

[TABLE]

where $[M_{\omega}]_{ij}:=(w_{i},w_{j})_{\omega}$ and $[A_{\omega}]_{ij}:=a_{\omega}(w_{i},w_{j})$, for $i,j=1,\dots,N$. Using Theorem 6.1 we can immediately obtain the following estimate for the spectral radius and therefore the time step size.

Theorem 6.3**.**

Let $\mathcal{W}_{h}$ be a weighted mesh decomposition. Then

[TABLE]

where $M_{\omega}^{\prime}:=(I_{\omega}^{*})^{t}M_{\omega}I_{\omega}^{*}$, $A_{\omega}^{\prime}:=(I_{\omega}^{*})^{t}A_{\omega}I_{\omega}^{*}$, and $I_{\omega}^{*}:=\mathcal{I}^{*}(M_{\omega})$, and where $\lambda_{max}(\cdot)$ denotes the largest eigenvalue in magnitude.

Remark 6.4**.**

When the weighted submeshes $\omega$ are nonzero for only a few elements and faces, then $M^{\prime}_{\omega}$ and $A^{\prime}_{\omega}$ are relatively small matrices. The largest eigenvalue \lambda_{max}\big{(}(M_{\omega}^{\prime})^{-1}A_{\omega}^{\prime}\big{)} can then be efficiently computed in parallel for each submatrix, using a power iteration method requiring only a relatively small number of iterations.

In the next section we show several numerical results illustrating the sharpness of the penalty term and time step estimates.

7. Numerical results

7.1. Computing the spectral radius for periodic meshes

To test the sharpness of the penalty term estimates and time step estimates we consider a $d$-dimensional cubic domain of the form $(0,N)^{d}$ with periodic boundary conditions. We then create a uniform mesh of $N^{d}$ unit cubes, after which we subdivide every cube into smaller elements and choose basis function sets and material parameters for every subelement. These subelements, basis functions and material parameters are chosen identically for every cube. An illustration of such a mesh is given in Figure 3. The advantage of such a uniform periodic mesh is that we can rather easily obtain the exact spectral radius by using a Fourier analysis, in a way similar to the von Neumann method for finite difference schemes.

To apply a Fourier analysis we first choose a linear basis $\{w_{i}\}_{i=1}^{N^{d}\times M}$ of the discrete function space $U$, such that the linear basis is of the form

[TABLE]

where $\mathbf{k}\in\mathbb{Z}_{N}^{d}$ is the identifier of unit cube $(k_{1},k_{1}+1)\times\cdots\times(k_{d},k_{d}+1)$ and $\{w_{\mathbf{k},i}\}_{i=1}^{M}$ is a linear basis of the discrete space $U$ restricted to this cube. We can then define the submatrices $M_{\mathbf{k},\mathbf{l}},A_{\mathbf{k},\mathbf{l}}\in\mathbb{R}^{M\times M}$ as follows:

[TABLE]

By construction of the mesh, most of these submatrices are identical. Fix any $\mathbf{l}\in\mathbb{Z}_{N}^{d}$. Then the submatrices $M_{\mathbf{k},\mathbf{k}+\mathbf{l}}$ are identical for any $\mathbf{k}\in\mathbb{Z}_{N}^{d}$. The same holds for $A_{\mathbf{k},\mathbf{k}+\mathbf{l}}$. Moreover, by definition of the mass matrix, $M_{\mathbf{k},\mathbf{k}+\mathbf{l}}$ is only nonzero when $\mathbf{l}=\mathbf{0}$, and by construction of the stiffness matrix, $A_{\mathbf{k},\mathbf{k}+\mathbf{l}}$ is only nonzero when $|\mathbf{l}|\leq 1$. Therefore, we only have to consider the submatrices $M_{0}:=M_{\mathbf{k},\mathbf{k}}$, $A_{0}:=A_{\mathbf{k},\mathbf{k}}$, and $A_{i}^{\pm}:=A_{\mathbf{k},\mathbf{k}\pm\mathbf{e}_{i}}$ for $i=1,\dots,d$, where $\mathbf{k}$ is an arbitrary vector in $\mathbb{Z}_{N}^{d}$ and $\mathbf{e}_{i}$ is the unit vector in direction $i$.

Now let $\underline{u}\in\mathbb{R}^{N^{d}\times M}$ be a vector of coefficients, and let $\underline{u}_{\mathbf{k}}\in\mathbb{R}^{M}$ be the vector of coefficients corresponding to cube $\mathbf{k}$. Suppose that $\underline{w}=M^{-1}A\underline{u}$. We can then write

[TABLE]

Define $\underline{u}^{(\mathbf{z})}\in\mathbb{R}^{N^{d}\times M}$ as follows:

[TABLE]

where $\underline{u}_{0}\in\mathbb{R}^{M}$ is an arbitrary vector of coefficients corresponding to a single cube, $\mathbf{i}:=\sqrt{-1}$ is the imaginary number, and $\mathbf{z}\in\mathbb{Z}_{N}^{d}$ is a vector of integers. Then $\underline{w}^{(\mathbf{z})}:=M^{-1}A\underline{u}^{(\mathbf{z})}$ satisfies

[TABLE]

where

[TABLE]

From (25) and (26) it follows that if $(\lambda,\underline{u}_{0})$ is an eigenpair of $Z^{(\mathbf{z})}$, then $(\lambda,\underline{u}^{(\mathbf{z})})$ is an eigenpair of $M^{-1}A$. Since $Z^{(\mathbf{z})}$ has $M$ eigenpairs and since there are $N^{d}$ possible choices for $\mathbf{z}$, every eigenvalue of $M^{-1}A$ is an eigenvalue of $Z^{(\mathbf{z})}$ for some $\mathbf{z}\in\mathbb{Z}_{N}^{d}$. For the time step estimates we are only interested in the largest eigenvalue $\lambda_{max}(M^{-1}A)$, which we can then compute by

[TABLE]

For the numerical tests that we will present here, we have taken $N=2$, since in most cases the largest eigenvalue $\lambda_{max}(M^{-1}A)$ no longer increases significantly for $N>2$.

7.2. Sharpness of the penalty term and time step estimates

For testing the sharpness of our parameter estimates we use polynomial basis functions up to degree $p$ for simplicial elements, and polynomials up to degree $p$ in the direction of each reference coordinate for quadrilateral and hexahedral elements. First, we consider several regular homogeneous meshes for the acoustic wave equation in 1D, 2D and 3D. After that we test on meshes with deformed elements and meshes with piecewise linear parameter fields. We also test on meshes for electromagnetic and elastic wave problems, including heterogeneous meshes with sharp material contrasts and meshes with sharp contrasts in primary and secondary wave velocities.

To test the sharpness of the parameters we first compute the penalty terms as in Theorem 5.6 or Theorem 5.8 with $c_{\kappa}=1$. We will refer to the first penalty terms as $\eta_{e}^{*}$ and to the second as $\eta_{e}^{**}$. We then find the smallest scale $c_{min}\in[0,1]$ such that the stiffness matrix $A$ is still positive semidefinite when using the downscaled penalty terms $\eta_{min,e}:=c_{min}\eta_{e}$. We compute $c_{min}$ accurate to two decimal places using the bisection method.

After we have computed $\eta$ and $c_{min}$, we consider the time step condition for the leap–frog scheme and use this to compute the time step size in the three ways given below:

[TABLE]

Here $A_{min}$ is the stiffness matrix that results from using the downscaled penalty terms $\eta_{min,e}$, and $M^{\prime}_{\omega}$ and $A^{\prime}_{\omega}$ are the submatrices corresponding to the weighted submesh $\omega$. These time step sizes can be interpreted as follows: $\Delta t(\eta_{min})$ is the largest allowed time step size when using the minimum downscaled penalty terms $\eta_{min,e}$, $\Delta t(\eta)$ is the largest allowed time step size when using the penalty term estimates $\eta_{e}$, and $\Delta t_{est}(\eta)$ is the time step estimate when using the penalty term estimates $\eta_{e}$ and a weighted mesh decomposition.

For our time step estimate $\Delta t_{est}(\eta)$ we will use the vertex-based mesh decomposition as given in (23). We will measure the sharpness of the penalty term by $\Delta t(\eta_{min})/\Delta t(\eta)$, and we will measure the sharpness of our time step estimate by $\Delta t(\eta)/\Delta t_{est}(\eta)$.

7.2.1. Regular meshes

For the first tests we consider the acoustic wave equation as given in Example 3.1, with $c=1$. We use meshes of the form described in Section 7.1, with element subdivisions as listed below. An illustration of some of the element subdivisions is given in Figure 4.

• •

1D: mesh constructed from unit intervals.

• •

square: 2D mesh constructed from unit squares.

• •

triangular: 2D mesh constructed from unit squares, with each square subdivided into two triangles.

• •

cubic: 3D mesh constructed from unit cubes.

• •

tetrahedral: 3D mesh constructed from unit cubes, with each cube subdivided into six pyramids, and every pyramid subdivided into four tetrahedra.

The results of the parameter estimates, when using the penalty term $\eta^{*}$, are given in Table 1. The results when using $\eta^{**}$ are given in Table 2. From these tables we can already see that our second penalty term estimate $\eta^{**}$ is in general much sharper than the first estimate $\eta^{*}$. This is true especially for cubes, where $\eta^{*}$ causes a reduction in the largest allowed time step size of more than $2$. Also, for square and tetrahedral meshes the first penalty term estimate causes a reduction in the time step size of more than a factor $1.5$. On the other hand, when using $\eta^{**}$ the largest allowed time step size is never reduced more than a factor $1.2$, and for many of the regular meshes $\Delta t(\eta_{min}^{**})/\Delta t(\eta^{**})$ is even below $1.01$. For the 1D, square, and cubic meshes, this penalty term and corresponding time step estimate even coincide with the analytic results derived in [1]. In general, the time step estimate does not reduce the time step size by a factor more than $1.2$ when using $\eta^{**}$ and not more than $1.3$ when using $\eta^{*}$.

For the case of tetrahedral meshes we can compare our penalty term estimates with the estimate derived in [16]. The penalty term derived there is equivalent to $\eta_{e}=p(p+2)$, with the penalty scaling function given by $\nu_{h}|_{\partial e\cap f}=1/(\epsilon_{f}\min_{e\in{\mathcal{T}_{f}}}d_{i,e})$, where $d_{i,e}$ is the diameter of the inscribed sphere of $e$ and $\epsilon_{f}\in\{1/2,1\}$ is defined as in (8b). Their analysis can be readily extended to triangles by replacing the trace inverse inequality for tetrahedra by the trace inverse inequality for triangles, given in Theorem 3 of [18], which is equivalent to setting $\eta_{e}=p(p+1)$. The results of these estimates are given in Table 3, from which we can see that this penalty term estimate has a similar sharpness as $\eta^{*}$, but is significantly less sharp than $\eta^{**}$, having a time step size more than $1.5$ times smaller than when using $\eta^{**}$ for $p=2,3$ on tetrahedra and $p=3$ on triangles.

Since $\eta^{**}$ is significantly sharper than $\eta^{*}$, we will only use $\eta^{**}$ throughout the following numerical tests.

7.2.2. Meshes with deformed elements

In this subsection we consider the acoustic wave equation with $c=1$ again, but now using deformed elements. For the penalty term we will only use $\eta^{**}$. An overview of the different meshes is listed below, and an illustration of the element subdivisions is given in Figures 5 and 6.

• •

rectangular$[x]$: 2D mesh constructed from unit squares, with each square subdivided into $2\times 2$ rectangles adjacent to a central node at $(x,x)$.

• •

quadrilateral$[x]$: 2D mesh constructed from unit squares, with each square subdivided into $2\times 2$ smaller uniform squares, after which the central node at (0.5,0.5) is moved to $(x,x)$.

• •

triangular$[x]$: 2D mesh constructed from unit squares, with each square subdivided into four triangles adjacent to the central node at $(x,x)$.

• •

cuboid$[x]$: 3D mesh constructed from unit cubes, with each cube subdivided into $2\times 2\times 2$ cuboids adjacent to the central node at $(x,x,x)$.

• •

hexahedral$[x]$: 3D mesh constructed from unit cubes, with each cube subdivided into $2\times 2\times 2$ smaller uniform cubes, after which the central node at $(0.5,0.5,0.5)$ is moved to $(x,x,x)$.

• •

tetrahedral$[x]$: 3D mesh constructed from unit cubes, with each cube subdivided into six pyramids adjacent to the central node at $(x,x,x)$, and with each pyramid subdivided into four tetrahedra.

The results of the parameter estimates are given in Tables 4 and 5. In all cases, the penalty term estimate causes a reduction in the time step size of no more than a factor $1.1$, with respect to the time step size using the minimal downscaled penalty term $\eta_{min}^{**}$. The time step estimates for triangular meshes causes a reduction of no more than $1.25$ and for the tetrahedral meshes a reduction of no more than $1.2$ with respect to the largest allowed time step size. For quadrilateral and hexahedral meshes the time step estimate tends to get less sharp for higher order polynomial basis functions and more strongly deformed meshes, but in all of our tests $\Delta t(\eta^{**})/\Delta t_{est}(\eta^{**})$ remains below $1.4$.

Since it is hard to compute the element and face integrals for general quadrilateral and hexahedral meshes exactly, we approximate them using the Gauss–Legendre quadrature rule with $p+1$ points in every direction, where $p$ is the polynomial order of the basis functions. We do not consider meshes of the type quadrilateral$[x]$ with $x\geq 0.75$ and hexahedral$[x]$ with $x\geq 2/3$ since in those cases one of the elements no longer has a well-defined mapping.

7.2.3. Meshes with piecewise linear parameters

We now consider the acoustic wave equation with piecewise linear parameter fields $\rho$ and $c$ instead of constant parameters. An overview of the different meshes is listed below.

• •

squarePL$[\rho_{0},c_{0}]$: 2D mesh constructed from unit squares, with each square subdivided into $2\times 2$ smaller squares and with piecewise linear parameters $\rho,c$ such that $\rho=c=1$ at $y=0$ and $y=1$ and $\rho=\rho_{0},c=c_{0}$ at $y=0.5$.

• •

triPL$[\rho_{0},c_{0}]$: 2D mesh constructed from unit squares, with each square subdivided into four uniform triangles and with piecewise linear parameters $\rho,c$ such that $\rho=c=1$ at the boundary and $\rho=\rho_{0},c=c_{0}$ at the center of each square.

• •

cubicPL$[\rho_{0},c_{0}]$: 3D mesh constructed from unit cubes, with each cube subdivided into $2\times 2\times 2$ smaller cubes and with piecewise linear parameters $\rho,c$ such that $\rho=c=1$ at $z=0$ and $z=1$ and $\rho=\rho_{0},c=c_{0}$ at $z=0.5$.

• •

tetrahedralPL$[\rho_{0},c_{0}]$: 3D mesh constructed from unit cubes, with each cube subdivided into $24$ uniform tetrahedra and with piecewise linear parameters $\rho,c$ such that $\rho=c=1$ at the boundary and $\rho=\rho_{0},c=c_{0}$ at the center of each cube.

The results of the parameter estimates are given in Table 6. The penalty term estimate cause a reduction in the time step size of no more than a factor $1.2$, with respect to the time step size using the minimal downscaled penalty term $\eta_{min}^{**}$ for triangular and tetrahedral elements and no reduction at all for square and cubic elements. The time step estimates for the triangular and tetrahedral meshes causes a reduction in the time step size of no more than $1.2$ with respect to the largest allowed time step size, while for the square and cubic meshes they often cause no reduction at all or a reduction not larger than $1.05$.

7.2.4. Meshes for electromagnetic and elastic wave problems

In this subsection we consider the electromagnetic wave equations as given in Example 3.2 and the 3D isotropic elastic wave equations as given in Example 3.3. We use the electromagnetic wave equations to test heterogeneous media with sharp contrasts in material parameters, while we use the elastic wave equations to test media with large contrasts in primary and secondary wave velocities. An overview of the different meshes is listed below.

• •

tetraEM$[\mu_{1},\mu_{2}]$: 3D mesh constructed from unit cubes, with each cube subdivided into $24$ smaller uniform tetrahedra, where all tetrahedra below the surface $x=z$ have a relative magnetic permeability of $\mu_{1}$ and all tetrahedra above this surface have a permeability of $\mu_{2}$. For all tetrahedra the relative electric permittivity $\epsilon$ equals $1$.

• •

cubicEM$[\mu_{1},\mu_{2}]$: 3D mesh constructed from unit cubes, with each cube subdivided into $2\times 2\times 2$ smaller uniform cubes, where the bottom four cubes have a relative magnetic permeability of $\mu_{1}$ and the top four cubes a permeability of $\mu_{2}$. For all cubes the relative electric permittivity $\epsilon$ equals $1$.

• •

tetraISO$[\lambda,\mu]$: 3D mesh constructed from unit cubes, with each cube subdivided into $24$ smaller uniform tetrahedra, and with all tetrahedra having Lamé parameters $\lambda$ and $\mu$ and mass density $\rho=1$.

• •

cubicISO$[\lambda,\mu]$: 3D mesh constructed from unit cubes, with each cube having Lamé parameters $\lambda$ and $\mu$ and mass density $\rho=1$.

The results of the parameter estimates are given in Tables 7 and 8. For the heterogeneous electromagnetic problems, the penalty term estimate causes a reduction in the time step size of no more than a factor $1.1$, with respect to the time step size using the minimal downscaled penalty term $\eta_{min}^{**}$. For the isotropic elastic meshes this reduction is no more than a factor $1.2$. The time step estimates for all the tetrahedral meshes causes a reduction in the time step size of no more than $1.25$ with respect to the largest allowed time step size, while for the cubic meshes this reduction remains below a factor $1.4$.

8. Conclusion

We have derived sharp and sufficient conditions for the penalty parameter in Theorem 5.6 and Theorem 5.8 to guarantee stability of the SIPDG method for linear wave problems. In addition, we derived sufficient upper bounds for the spectral radius and the time step size in Theorem 6.3 and Section 6.3, by introducing a weighted mesh decomposition. These conditions hold for generic meshes, including unstructured nonconforming heterogeneous meshes of mixed element types, with different types of boundary conditions. Moreover, the estimates hold for general linear hyperbolic partial differential equations, including the acoustic wave equation, the (an)isotropic elastic wave equations, and Maxwell’s equations. Both the penalty term and time step size can be efficiently computed in parallel using a power iteration method on each element and submesh, respectively.

To test the sharpness of our estimates we have considered several semiuniform meshes made of $N^{d}$ uniform cubes or squares, which are then uniformly subdivided into smaller meshes, including meshes with deformed elements and heterogeneous media. From the results we can see that the penalty term estimate given in Theorem 5.8 is significantly sharper than the estimate in Theorem 5.6 for 2D and 3D meshes. Especially for cubic elements, we see that downscaling the penalty estimate constructed using Theorem 5.6 allows a time step size more than twice as large, while the penalty term in Theorem 5.8 does not allow any downscaling at all. Furthermore, the penalty term estimate from Theorem 5.8 allows a time step size more than a factor $1.5$ times larger than when using the penalty term estimate derived in [16], when using square or cubic basis functions on tetrahedra or cubic basis functions on triangles. For regular square and cubic meshes, this same penalty term estimate, together with the time step estimate from Theorem 6.3 using the weighted mesh decomposition of (23), conforms with the largest allowed time step sizes analytically derived in [1].

In general we see that downscaling the penalty term from Theorem 5.8 does not increase the maximum allowed time step size by more than a factor $1.2$, even when the material parameters are nonconstant within the elements, while our time step estimate, based on a vertex-based weighted mesh decomposition, does not become smaller than a factor $1.2$ compared to the maximum allowed time step size for triangular and tetrahedral meshes, and not smaller than a factor $1.4$ for quadrilateral and hexahedral meshes.

Acknowledgments

We would like to thank Shell Global Solutions International B.V. for supporting this project, and in particular Wim Mulder for giving us useful advice and sharing his practical experiences related to the topic.

Appendix A Linear Algebra

Lemma A.1**.**

Let $C\in\mathbb{R}^{d\times m\times m\times d}_{sym}$ be a symmetric tensor such that

[TABLE]

Then there exists a symmetric tensor field $C^{1/2}\in\mathbb{R}_{sym}^{d\times m\times m\times d}$ such that $C^{1/2}:C^{1/2}=C$.

Proof.

The result follows from the fact that $\boldsymbol{\sigma}\rightarrow C:\boldsymbol{\sigma}$ is a linear, self-adjoint, and positive semidefinite operator on a finite-dimensional subspace $\mathbb{R}^{d\times m}$. ∎

Lemma A.2**.**

Let $\hat{n}\in\mathbb{R}^{d}$ be a vector and $C\in\mathbb{R}^{d\times m\times m\times d}_{sym}$ be a symmetric tensor such that

[TABLE]

Also, define the second order tensor $\mathbf{c}_{\hat{n}}:=\hat{n}\cdot C\cdot\hat{n}$ and the following function space:

[TABLE]

There exists a pseudoinverse $\mathbf{c}_{\hat{n}}^{-1}\in\mathbb{R}^{m\times m}$ such that $\mathbf{c}_{\hat{n}}^{-1}\cdot\mathbf{c}_{\hat{n}}\cdot u=\mathbf{c}_{\hat{n}}\cdot\mathbf{c}_{\hat{n}}^{-1}\cdot u=u$, for all $u\in U$. Moreover, both $\mathbf{c}_{\hat{n}}$ and $\mathbf{c}_{\hat{n}}^{-1}$ are symmetric and positive semidefinite.

Proof.

The fact that $\mathbf{c}_{\hat{n}}$ is symmetric follows from its definition and the fact that $C$ is symmetric. The fact that this tensor is also positive semidefinite follows from (28) and by writing

[TABLE]

To prove that the pseudoinverse with respect to $U$ exists and is symmetric positive semidefinite, it is sufficient to show that

[TABLE]

Since $\mathbf{c}_{\hat{n}}$ is positive semidefinite we then only need to show that

[TABLE]

Now take an arbitrary $u\in U$, and suppose that $u\cdot\mathbf{c}_{\hat{n}}\cdot u=0$. Using Lemma A.1 we can write

[TABLE]

From this it follows that $C:\hat{n}u=\mathbf{0}$. Because of the definition of $U$ we can write $u=\hat{n}\cdot C:\boldsymbol{\sigma}$ for some $\boldsymbol{\sigma}\in\mathbb{R}^{d\times m}$, and therefore $C:\hat{n}(\hat{n}\cdot C:\boldsymbol{\sigma})=\mathbf{0}$. Applying the double dot product with $\boldsymbol{\sigma}^{t}$ on the left side gives

[TABLE]

Therefore, $u=0$, which proves (30). ∎

Lemma A.3**.**

Let $\{a_{i}\}_{i=1}^{n}$ be a set of nonnegative scalars, and let $\{b_{i}\}_{i=1}^{n}$ be a set of positive scalars. Then

[TABLE]

Proof.

[TABLE]

∎

Lemma A.4**.**

Let $A\in\mathbb{R}^{n\times n}_{sym}$ be a symmetric matrix and $M\in\mathbb{R}^{n\times n}_{sym}$ a positive definite matrix. Then there exists a diagonalization $M^{-1}A=VDV^{-1}$, such that $D$ is a real diagonal matrix and $V$ satisfies $V^{t}MV=I$, where $I$ is the identity matrix.

Proof.

Since $M$ is symmetric positive definite, there exists a symmetric positive definite matrix $M^{1/2}$ such that $M^{1/2}M^{1/2}=M$. Define $\tilde{A}:=M^{-1/2}AM^{-1/2}$. Since $A$ is symmetric, $\tilde{A}$ is symmetric as well and can be diagonalised as $\tilde{A}=\tilde{V}D\tilde{V}^{-1}$ with $D$ a diagonal real matrix and $\tilde{V}$ satisfying $\tilde{V}^{t}\tilde{V}=I$. Now define $V:=M^{-1/2}\tilde{V}$. Then $M^{-1}A=V^{-1}DV$ and $V^{t}MV=I$. ∎

Lemma A.5**.**

Let $A\in\mathbb{R}^{n\times n}_{sym}$ be a symmetric matrix and $M\in\mathbb{R}^{n\times n}_{sym}$ a positive definite matrix. Then

[TABLE]

where \lambda_{max}\big{(}M^{-1}A\big{)} denotes the largest eigenvalue of $M^{-1}A$ in magnitude.

Proof.

From Lemma A.4 it follows that there exist a diagonal matrix $D$ and a nonsingular matrix $V$ such that $MVDV^{-1}=A$ and $V^{t}MV=I$. Now consider an arbitrary $\underline{u}\in\mathbb{R}^{n},\underline{u}\neq\underline{0}$. We set $\underline{w}:=V^{-1}\underline{u}$ to obtain

[TABLE]

The lemma follows from this equality and the following relation:

[TABLE]

∎