# 3.3.1 DG for the acoustic wave equation¶

We consider the acoustic problem \begin{align} \partial_{t} \mathbf{u} - c \nabla p &= 0 \quad \text{ in } \Omega \times I, \\ \partial_{t} p - c \operatorname{div}(\mathbf{u}) &= 0 \quad \text{ in } \Omega \times I, \\ % p (0,\cdot) &= p(1,\cdot), \quad % p (\cdot,0) = p(\cdot,1), \label{eq:per2b}\\ % \mathbf{q} (0,\cdot) &= \mathbf{q}(1,\cdot), \quad % \mathbf{q} (\cdot,0) = \mathbf{q}(\cdot,1), \label{eq:per2bb}\\ p &= p_0 \!\!\! \quad \text{ on } \Omega \times \{0\},\\ u &= 0 \quad \text{ on } \Omega \times \{0\}. \end{align}

Here $$p$$ is the acoustic pressure (the local deviation from the ambient pressure) and $$\mathbf{u}$$ is the local velocity.

+ suitable boundary conditions

A simple grid:

:

from ngsolve import *
# import netgen.gui
from ngsolve.webgui import Draw
from netgen.geom2d import SplineGeometry
geo = SplineGeometry()
geo.AddRectangle( (-1, -1), (1, 1), bcs = ("bottom", "right", "top", "left"))
mesh = Mesh( geo.GenerateMesh(maxh=0.1))
Draw(mesh)

:

BaseWebGuiScene


Find $$p: [0,T] \to \bigoplus_{T\in\mathcal{T}_h} \mathcal{P}^{k+1}(T)$$ and $$\mathbf{u}: [0,T] \to \bigoplus_{T\in\mathcal{T}_h} [\mathcal{P}^k(T)]^d$$ so that

\begin{align*} (\partial_t \mathbf{u}, v) &= b_h(p,v) & &&& \forall v,\\ (\partial_t p, q) &= & - b_h(q,\mathbf{u}) &&& \forall q\\ \end{align*}

with the centered flux approximation:

$b_h(p,v) = \sum_{T} \int_T \nabla p \cdot v + \int_{\partial T} ( \hat{p} - p) v_n$

Here $$\hat{p}$$ is the centered approximation i.e. $$\hat{p} = \{\!\!\{p\}\!\!\}$$.

:

order = 6
fes_p = L2(mesh, order=order+1, all_dofs_together=True)
fes_u = VectorL2(mesh, order=order, piola=True, all_dofs_together=True)
fes = fes_p*fes_u
gfu = GridFunction(fes)
sceneu = Draw(gfu.components, mesh, "u")
scenep = Draw(gfu.components, mesh, "p", min=-0.02, max=0.08, autoscale=False, order=3)


## What is the flag piola doing?¶

The VectorL2 space uses the following definition of basis functions on the mesh.

Let $$\hat{\varphi}(\hat{x})$$ be a (vectorial) basis function on the reference element $$\hat{T}$$, $$\Phi: \hat{T} \to T$$ with $$x = \Phi(\hat{x})$$ and $$F= D\Phi$$, then

$\begin{split}\varphi(x) := \left\{ \begin{array}{rlll} \hat{\varphi}(\hat{x}), & \texttt{piola=False}, & \texttt{covariant=False}, & \text{(default)} \\ \frac{1}{|\operatorname{det} F|} \cdot F \cdot \hat{\varphi}(\hat{x}), & \texttt{piola=True}, & \texttt{covariant=False}, \\ F^{-T} \cdot \hat{\varphi}(\hat{x}), & \texttt{piola=False}, & \texttt{covariant=True}. \\ \end{array} \right.\end{split}$

### Inverse mass matrix operations:¶

Combining Embedding and the inverse mass matrix operation for fes_p allows to realize the following block inverse operations acting on fes:

$\begin{split}\underbrace{\left( \begin{array}{cc} M_{p}^{-1} & 0 \\ 0 & 0 \end{array} \right)}_{\texttt{invp}} = \underbrace{\left( \begin{array}{c} I \\ 0 \end{array} \right)}_{\texttt{emb_p}} \cdot \underbrace{M_{p}^{-1}}_{\texttt{invmassp}} \cdot \underbrace{\left( \begin{array}{cc} I & 0 \end{array} \right)}_{\texttt{emb_p.T}}\end{split}$
:

pdofs = fes.Range(0);
emb_p = Embedding(fes.ndof, pdofs)
invmassp = fes_p.Mass(1).Inverse()
invp = emb_p @ invmassp @ emb_p.T


Analogously, combining Embedding and the inverse mass matrix operation for fes_u allows to realize the following block inverse operations on fes:

$\begin{split}\underbrace{\left( \begin{array}{cc} 0 & 0 \\ 0 & M_{\mathbf{u}}^{-1} \end{array} \right)}_{\texttt{invu}} = \underbrace{\left( \begin{array}{c} I \\ 0 \end{array} \right)}_{\texttt{emb_u}} \cdot \underbrace{M_{\mathbf{u}}^{-1}}_{\texttt{invmassu}} \cdot \underbrace{\left( \begin{array}{cc} I & 0 \end{array} \right)}_{\texttt{emb_u.T}}\end{split}$
:

udofs = fes.Range(1)
emb_u = Embedding(fes.ndof, udofs)
invmassu = fes_u.Mass(Id(mesh.dim)).Inverse()
invu = emb_u @ invmassu @ emb_u.T


### Time loop¶

Assuming the operators B,BT: fes_p$$\times$$ fes_u $$\to$$(fes_p$$\times$$ fes_u )$$'$$ corresponding to $$B_h((u,p),(v,q)=b_h(p,v)$$ are given, the symplectic Euler time stepping method takes the form:

:

def Run(B, BT, t0 = 0, tend = 0.1, tau = 1e-3, backward = False):
t = 0
while t < (tend-t0) - tau/2:
t += tau
if not backward:
gfu.vec.data += -tau * invp @ BT * gfu.vec
gfu.vec.data += tau * invu @ B * gfu.vec
print("\r t = {:}".format(t0 + t),end="")
else:
gfu.vec.data += -tau * invu @ B * gfu.vec
gfu.vec.data += tau * invp @ BT * gfu.vec
print("\r t = {:}".format(tend - t),end="")
sceneu.Redraw() # blocking=False)
scenep.Redraw() # blocking=False)
print("")


Initial values (density ring):

:

gfu.components.Set (exp(-50*(x**2+y**2))-exp(-100*(x**2+y**2)))
gfu.components.vec[:] = 0
scenep.Redraw()


#### The bilinear form for application on the full space fes¶

:

n = specialcf.normal(mesh.dim)
(p,u),(q,v) = fes.TnT()
B = BilinearForm(fes, nonassemble=True)
B += grad(p)*v * dx + 0.5*(p.Other()-p)*(v*n) * dx(element_boundary=True)
BT = BilinearForm(fes, nonassemble=True)
BT += grad(q)*u * dx + 0.5*(q.Other()-q)*(u*n) * dx(element_boundary=True)

:

%time Run(B.mat, BT.mat, backward=False)
%time Run(B.mat, BT.mat, backward=True)

 t = 0.10000000000000007
CPU times: user 5.79 s, sys: 49.1 ms, total: 5.84 s
Wall time: 1.46 s
t = -6.938893903907228e-17
CPU times: user 5.55 s, sys: 54.6 ms, total: 5.61 s
Wall time: 1.4 s


#### Using the TraceOperator and assembling of $$b_h$$¶

We want to use the TraceOperator again:

:

fes_tr = FacetFESpace(mesh, order=order+1)
traceop = fes_p.TraceOperator(fes_tr, False)


We want to assemble the sub-block matrix and need test/trial functions for single spaces (not the product spaces):

:

p = fes_p.TrialFunction()
v = fes_u.TestFunction()
phat = fes_tr.TrialFunction()


We split the operator to $$b_h$$ into * volume contributions (local) * and couplings between the trace (obtained through the trace op) and the volume:

:

Bel = BilinearForm(trialspace=fes_p, testspace=fes_u)
Bel += grad(p)*v * dx -p*(v*n) * dx(element_boundary=True)
%time Bel.Assemble()

Btr = BilinearForm(trialspace=fes_tr, testspace=fes_u)
Btr += 0.5 * phat * (v*n) * dx(element_boundary=True)
%time Btr.Assemble()

B = emb_u @ (Bel.mat + Btr.mat @ traceop) @ emb_p.T

CPU times: user 124 ms, sys: 7.94 ms, total: 132 ms
Wall time: 131 ms
CPU times: user 45.4 ms, sys: 7.3 ms, total: 52.7 ms
Wall time: 52.4 ms


Version 2 in action (assembled matrices, TraceOperators and Embeddings can do transposed multiply):

:

%time Run(B, B.T, backward=False)
%time Run(B, B.T, backward=True)

 t = 0.10000000000000007
CPU times: user 4.96 s, sys: 1.88 s, total: 6.85 s
Wall time: 1.71 s
t = -6.938893903907228e-17
CPU times: user 4.65 s, sys: 2.02 s, total: 6.68 s
Wall time: 1.67 s


#### Using the TraceOperator and geom_free=True with assembling¶

With $$p(x) = \hat{p}(\hat{x})$$ and $$\mathbf{v}(x) = \frac{1}{|\operatorname{det}(F)|} F \cdot \hat{\mathbf{v}}(\hat{x})$$ (Piola mapped) there holds

\begin{align} \int_T \mathbf{v} \cdot \nabla p = \int_{\hat{T}} \hat{\mathbf{v}} \cdot \hat{\nabla} \hat{p} \end{align}

and similarly for the facet integrals, cf. unit-2.11 and

Integrals are independent of the “physical element” (up to ordering)

$$\leadsto$$ same element matrices for a large class of elements

:

Bel = BilinearForm(trialspace=fes_p, testspace=fes_u, geom_free = True)
Bel += grad(p)*v * dx -p*(v*n) * dx(element_boundary=True)
%time Bel.Assemble()

Btr = BilinearForm(trialspace=fes_tr, testspace=fes_u, geom_free = True)
Btr += 0.5 * phat * (v*n) *dx(element_boundary=True)
%time Btr.Assemble()

B = emb_u @ (Bel.mat + Btr.mat @ traceop) @ emb_p.T

CPU times: user 1.67 ms, sys: 268 µs, total: 1.94 ms
Wall time: 857 µs
CPU times: user 1.75 ms, sys: 0 ns, total: 1.75 ms
Wall time: 1.52 ms


Version 3 in action:

:

%time Run(B, B.T, backward=False)
%time Run(B, B.T, backward=True)

 t = 0.10000000000000007
CPU times: user 1.28 s, sys: 50 ms, total: 1.33 s
Wall time: 330 ms
t = -6.938893903907228e-17
CPU times: user 1.28 s, sys: 39.7 ms, total: 1.31 s
Wall time: 328 ms


References:

• [Hesthaven, Warburton. Nodal Discontinuous Galerkin Methods—Algorithms, Analysis and Applications, 2007]