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

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:

[1]:

from ngsolve import *
import netgen.gui
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)

 Generate Mesh from spline geometry
 Boundary mesh done, np = 80
 CalcLocalH: 80 Points 0 Elements 0 Surface Elements
 Meshing domain 1 / 1
 load internal triangle rules
 Surface meshing done
 Edgeswapping, topological
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Edgeswapping, metric
 Smoothing
 Split improve
 Combine improve
 Smoothing
 Update mesh topology
 Update clusters

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\}\!\!\}\).

[2]:

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 = FESpace( [fes_p,fes_u] )
gfu = GridFunction(fes)
Draw(gfu.components[1], mesh, "u")
Draw(gfu.components[0], mesh, "p", min=-0.02, max=0.08, autoscale=False, sd=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}\]
[3]:

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}\]
[4]:

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:

[5]:

def Run(B, BT, t0 = 0, tend = 0.1, tau = 1e-3, backward = False):
    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="")
            Redraw(blocking=False)
        print("")

Initial values (density ring):

[6]:

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

setvalues element 930/930

Version 1:

The bilinear form for application on the full space fes

[7]:

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)

[8]:

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

 t = 0.10000000000000007
CPU times: user 5.45 s, sys: 125 ms, total: 5.57 s
Wall time: 1.39 s
 t = -6.938893903907228e-17
CPU times: user 5.31 s, sys: 96 ms, total: 5.41 s
Wall time: 1.35 s

Version 2:

Using the TraceOperator and assembling of \(b_h\)

We want to use the TraceOperator again:

[9]:

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):

[10]:

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:

[11]:

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

assemble mixed bilinearform
CPU times: user 122 ms, sys: 8.2 ms, total: 130 ms
Wall time: 128 ms
assemble mixed bilinearform
CPU times: user 40.7 ms, sys: 12.1 ms, total: 52.9 ms
Wall time: 51.4 ms

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

[12]:

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

 t = 0.10000000000000007
CPU times: user 4.75 s, sys: 2.11 s, total: 6.86 s
Wall time: 1.71 s
 t = -6.938893903907228e-17
CPU times: user 4.67 s, sys: 2.02 s, total: 6.69 s
Wall time: 1.66 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}{|\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

[13]:

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.33 ms, sys: 349 µs, total: 1.68 ms
Wall time: 877 µs
CPU times: user 2.24 ms, sys: 0 ns, total: 2.24 ms
Wall time: 1.77 ms

Version 3 in action:

[14]:

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

 t = 0.10000000000000007
CPU times: user 1.24 s, sys: 109 ms, total: 1.35 s
Wall time: 334 ms
 t = -6.938893903907228e-17
CPU times: user 1.13 s, sys: 43.5 ms, total: 1.18 s
Wall time: 293 ms

References: