This page was generated from unit-3.5.2-wavedg/wavedg.ipynb.

3.5.2 DG for the acoustic wave equation

We consider the acoustic problem $$\begin{array}{rl}{\partial }_t\mathbf{u}-c\mathrm{\nabla }p& =0\text{ in }\mathrm{\Omega }\times I,\\ {\partial }_{t}p-c\mathrm{div}(\mathbf{u})& =0\text{ in }\mathrm{\Omega }\times I,\\ p& =p_0\text{ on }\mathrm{\Omega }\times \{0\},\\ u& =0\text{ on }\mathrm{\Omega }\times \{0\}.\end{array}$$

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:

#imports, geometry and mesh:
from netgen.occ import *
from netgen.webgui import Draw as DrawGeo
geo = OCCGeometry(unit_square_shape.Scale((0,0,0),2).Move((-1,-1,0)), dim=2)
from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(geo.GenerateMesh(maxh=0.1))
Draw(mesh)

BaseWebGuiScene

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

$$\begin{array}{rlrlrl}({\partial }_t\mathbf{u},v)& =b_h(p,v)& & & & \mathrm{\forall }v,\\ ({\partial }_{t}p,q)& =& -b_h(q,\mathbf{u})& & & \mathrm{\forall }q\end{array}$$

with the centered flux approximation:

$$b_h(p,v)=\sum _T{\int }_T\mathrm{\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)

What is the flag piola doing?

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

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

$$\begin{array}{r}\phi (x):=\{\begin{array}{rlll}\hat{\phi }(\hat{x}),& \mathtt{\text{piola=False}},& \mathtt{\text{covariant=False}},& \text{(default)}\\ \frac{1}{|\det F|}\cdot F\cdot \hat{\phi }(\hat{x}),& \mathtt{\text{piola=True}},& \mathtt{\text{covariant=False}},\\ F^{-T}\cdot \hat{\phi }(\hat{x}),& \mathtt{\text{piola=False}},& \mathtt{\text{covariant=True}}.\end{array}\end{array}$$

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{array}{r}\underset{\mathtt{\text{invp}}}{\underbrace{\left(\begin{array}{cc}{M}_p^{-1}& 0\\ 0& 0\end{array}\right)}}=\underset{\mathtt{\text{emb\_p}}}{\underbrace{\left(\begin{array}{c}I\\ 0\end{array}\right)}}\cdot \underset{\mathtt{\text{invmassp}}}{\underbrace{M_p^{-1}}}\cdot \underset{\mathtt{\text{emb\_p.T}}}{\underbrace{\left(\begin{array}{cc}I& 0\end{array}\right)}}\end{array}$$

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{array}{r}\underset{\mathtt{\text{invu}}}{\underbrace{\left(\begin{array}{cc}0& 0\\ 0& M_{\mathbf{u}}^{-1}\end{array}\right)}}=\underset{\mathtt{\text{emb\_u}}}{\underbrace{\left(\begin{array}{c}I\\ 0\end{array}\right)}}\cdot \underset{\mathtt{\text{invmassu}}}{\underbrace{M_{\mathbf{u}}^{-1}}}\cdot \underset{\mathtt{\text{emb\_u.T}}}{\underbrace{\left(\begin{array}{cc}I& 0\end{array}\right)}}\end{array}$$

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 )${}^{\prime }$ 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.25, tau = 1e-3,
        backward = False, scenes = []):
    t = 0
    with TaskManager():
        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="")
            for sc in scenes: sc.Redraw() # blocking=False)
        print("")

Initial values (density ring):

gfu.components[0].Set (exp(-50*(x**2+y**2))-exp(-100*(x**2+y**2)))
gfu.components[1].vec[:] = 0
Draw(gfu.components[0], mesh, "p", min=-0.02, max=0.08, autoscale=False, order=3)

BaseWebGuiScene

Version 1:

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)

scenep=Draw(gfu.components[0], mesh, "p", min=-0.02, max=0.08,
            autoscale=False, order=3, deformation=True)
%time Run(B.mat, BT.mat, backward=False, scenes=[scenep])
%time Run(B.mat, BT.mat, backward=True, scenes=[scenep])

 t = 0.25000000000000017
CPU times: user 15 s, sys: 155 ms, total: 15.2 s
Wall time: 3.8 s
 t = -1.6653345369377348e-16
CPU times: user 14.9 s, sys: 177 ms, total: 15 s
Wall time: 3.76 s

Version 2:

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 114 ms, sys: 24.1 ms, total: 138 ms
Wall time: 137 ms
CPU times: user 38.8 ms, sys: 15.9 ms, total: 54.7 ms
Wall time: 54.5 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.25000000000000017
CPU times: user 12.5 s, sys: 5.29 s, total: 17.7 s
Wall time: 4.44 s
 t = -1.6653345369377348e-16
CPU times: user 12.5 s, sys: 5.19 s, total: 17.7 s
Wall time: 4.44 s

Version 3:

Using the TraceOperator and geom_free=True with assembling

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

$$\begin{array}{r}{\int }_T\mathbf{v}\cdot \mathrm{\nabla }p={\int }_{\hat{T}}\hat{\mathbf{v}}\cdot \hat{\mathrm{\nabla }}\hat{p}\end{array}$$

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

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

$⇝$ 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 2.13 ms, sys: 0 ns, total: 2.13 ms
Wall time: 1.41 ms
CPU times: user 1.63 ms, sys: 220 μs, total: 1.85 ms
Wall time: 1.63 ms

Version 3 in action:

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

 t = 0.25000000000000017
CPU times: user 3.14 s, sys: 154 ms, total: 3.3 s
Wall time: 824 ms
 t = -1.6653345369377348e-16
CPU times: user 3.22 s, sys: 78.1 ms, total: 3.29 s
Wall time: 823 ms

References:

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

• [Koutschan, Lehrenfeld, Schöberl. Computer algebra meets finite elements: An efficient implementation for Maxwells equations. , 2012]