5.1. Discontinuous Galerkin methods for the wave equation

We want to solve the first order system

(5.1)\[\begin{split}\begin{aligned} \partial_tv+\nabla p &= f,&\text{in }&[0,T_0]\times \Omega,\\ \partial_t p+\nabla\cdot v &=0,&\text{in }&[0,T_0]\times \Omega,\\ p(0,\cdot)&=p_0,&\text{in }&\Omega,\\ v(0,\cdot)&=v_0,&\text{in }&\Omega,\\ v\cdot n&=0,&\text{in }&[0,T_0]\times \partial\Omega,\\ \end{aligned}\end{split}\]

numerically, using The Leap Frog time stepping. Since the computational effort depends mainly on the efficiency of the factorization of the mass matrices and the application of their inverses it would be benificial to use piecewise polynomial spaces, with basis functions supported only on single elements.

Let \(\mathcal T_h\) be a decomposition of the domain \(\Omega\) into triangles/tetrahedra such that

\[\bigcup_{T\in\mathcal T_h}\bar T = \bar \Omega.\]

Then we want to use the discrete spaces

(5.2)\[\begin{split}\begin{aligned} \mathcal W_h&:=\left\{p\in L^2(\Omega): p|_T\in\mathcal P^k(T)\, \forall T\in\mathcal T_h\right\},\\ \mathcal V_h&:=\left\{v\in L^2(\Omega)^2: v|_T\in\mathcal P^{k-1}(T)^2\, \forall T\in\mathcal T_h\right\}, \end{aligned}\end{split}\]

where \(\mathcal P^k(T)=\mathrm{span}\{p:T\to\mathbb R: (x,y)\mapsto x^ly^m, l+m\leq k\}\) for a fixed polynomial order \(k\in\mathbb R\).

However, since \(\mathcal W_h\subsetneq H^1(\Omega)\) and \(\mathcal V_h\subsetneq H(\mathrm{div})(\Omega)\) the usual derivation of weak formulations fails. A naive attempt to just insert a basis of the discrete spaces (5.2) into a weak formulation like (4.12) would yield a familiy of decoupled local problems on the single elements \(T\) with homogenuous boundary conditions but not an approxmation of the global problem.

This problem is fixed by modifying the weak forms \(b,\tilde b\) of the gradient and divergence, which for \(p\in H^1(\Omega),v\in H(\mathrm{div})(\Omega)\) and homogeneous Neumann boundary conditions should satisfy

(5.3)\[\begin{split}\begin{aligned} b(p,v) &= \int_\Omega \nabla p\cdot v,\\ \tilde b(v,p) &= \int_\Omega \nabla\cdot v p=\int_\Omega v\cdot\nabla p = b(p,v). \end{aligned}\end{split}\]

To motivate the following we consider the one-dimensional case first.

Numerical fluxes in one dimension

Consider a one-dimensional mesh \(\mathcal T = \{T_j = (a_{j-1},a_{j}), j = 1,\ldots N\}\) for \(a_0<a_1<\ldots<a_N\). Since \(p\in \mathcal W_h\) is still piecewise smooth we may apply the gradient on each element. To make up for the lost information across the element boundaries we add a linear combination of the inner and outer boundary values at \(a_j,a_{j-1}\) with respect to each element, i.e.

(5.4)\[b(p,v) := \sum_{j=1}^N\int_{T_j} \nabla p\cdot v- \left(\alpha p_+(a_{j-1})+\beta p_-(a_{j-1})\right) v(a_{j-1})+\left(\alpha p_-(a_{j})+\beta p_+(a_{j})\right) v(a_j)\]

with \(p_\pm(x)=\lim_{\varepsilon\to 0} p(x\pm\varepsilon)\), where we test with a smooth function \(v\). We want a consistent bilinear form, i.e., if we plug in \(p\in H^1(\Omega)\) we want to obtain the weak gradient from (5.3). It immediately follows that for this hold we need to set \(\beta = -\alpha\). Using the notation \(\{p\}:=\frac{1}{2}(p_++p_-)\) we may rewrite (5.4) with \(\beta=-\alpha\) by

(5.6)\[b(p,v) = \sum_{j=1}^N\int_{T_j} \nabla p\cdot v- 2\alpha \left(\{p\}(a_{j-1})-p_+(a_{j-1})\right) v(a_{j-1})+ 2\alpha \left(\{p\}(a_{j})-p_-(a_{j})\right) v(a_{j})\]

Now note that the definition of \(b\) in (5.4) is also valid for \(v\in L^2(\Omega)^2\) (i.e., also for \(v\in\mathcal V_h\)) if we define it locally on each element. To obtain an energy preserving method we need \(b(p,v)=\tilde b(v,p)\). Moreover we also want consistency of the \(\tilde b\) bilinear form, i.e., if we insert \(v\in H^1(\Omega)\) to obtain the weak divergence from (5.3). Applying integration by parts in (5.6)

(5.6)\[\begin{split}\begin{aligned} b(p,v) &= \sum_{j=1}^N\int_{T_j} \nabla p\cdot v- 2\alpha \left(\{p\}(a_{j-1})-p_+(a_{j-1})\right) v(a_{j-1})+ 2\alpha \left(\{p\}(a_{j})-p_-(a_{j})\right) v(a_{j})\\ &= \sum_{j=1}^N-\int_{T_j} p\nabla\cdot v- \left(2\alpha \{p\}(a_{j-1})+(1-2\alpha)p_+(a_{j-1})\right)v(a_{j-1})+ \left(2\alpha \{p\}(a_{j})-(2\alpha-1)p_-(a_{j})\right) v(a_{j}) \end{aligned}\end{split}\]

which will be consistent for \(\alpha=1/2\).

Jumps and averages

For the discontinuous Galerkin (DG) spaces (5.2) in higher dimensions we may define the average and jump operators as follows:

\[\begin{split}\begin{aligned} \{\cdot\}&:\begin{cases} W_h&\to L^2(\mathcal S)\\ p&\mapsto \left(x\mapsto \frac{1}{2}\lim_{\varepsilon\to 0}p(x+\varepsilon n)+p(x-\varepsilon n)\right) \end{cases},\\ [\cdot]&:\begin{cases} V_h&\to L^2(\mathcal S)\\ v&\mapsto \left(x\mapsto \frac{1}{2}\lim_{\varepsilon\to 0}v(x+\varepsilon n)\cdot n-v(x-\varepsilon n)\cdot n\right) \end{cases}, \end{aligned}\end{split}\]

where \(n\) is the outward normal and the skeleton \(\mathcal S\) is defined by

\[\mathcal S:=\bigcup_{T\in\mathcal T}\partial T.\]

Remark 5.1

The average and jump are independent of the choice of the normal vector \(n\).

On the domain boundary the averaging operator can be adapted to fit boundary conditions

DG in higher dimensions

Following the reasoning for the 1d case we obtain the general DG formulation

(5.7)\[\begin{split}\begin{aligned} \int_\Omega\partial_tv w+b(p,w) &= \int_\Omega f w,\\ \int_\Omega \partial_t pq-b(q,v) &=0,\\ p(0,\cdot)&=p_0,\\ v(0,\cdot)&=v_0, \end{aligned}\end{split}\]

for all \(w\in\mathcal V_h,q\in\mathcal W_h\) with

\[b(p,w) = \sum_{T\in\mathcal T} \int_T \nabla p\cdot w + \int_{\partial T}(\{p\}-p)w\]

Theorem 5.1

The DG formulation (4.10) is consistent, i.e., the weak solution \(p,v\) of (4.10) also solves (5.7)