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{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]:
#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)
[1]:
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:
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 = 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{\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
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
:
[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
:
[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.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):
[6]:
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)
[6]:
BaseWebGuiScene
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]:
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 14.8 s, sys: 103 ms, total: 14.9 s
Wall time: 3.74 s
t = -1.6653345369377348e-16
CPU times: user 14.7 s, sys: 125 ms, total: 14.8 s
Wall time: 3.7 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
CPU times: user 123 ms, sys: 16.2 ms, total: 139 ms
Wall time: 138 ms
CPU times: user 42.9 ms, sys: 12 ms, total: 54.9 ms
Wall time: 54.5 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.25000000000000017
CPU times: user 13 s, sys: 5.47 s, total: 18.4 s
Wall time: 4.61 s
t = -1.6653345369377348e-16
CPU times: user 12.4 s, sys: 5.19 s, total: 17.6 s
Wall time: 4.4 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.43 ms, sys: 0 ns, total: 1.43 ms
Wall time: 1.26 ms
CPU times: user 2.11 ms, sys: 0 ns, total: 2.11 ms
Wall time: 1.89 ms
Version 3 in action:
[14]:
%time Run(B, B.T, backward=False)
%time Run(B, B.T, backward=True)
t = 0.25000000000000017
CPU times: user 3.12 s, sys: 125 ms, total: 3.25 s
Wall time: 811 ms
t = -1.6653345369377348e-16
CPU times: user 3.1 s, sys: 143 ms, total: 3.25 s
Wall time: 811 ms
References:
[Hesthaven, Warburton. Nodal Discontinuous Galerkin Methods—Algorithms, Analysis and Applications, 2007]