6.4. Cartesian time-domain PMLs in ngsolve#
Using netgen.geom2d the PML domain can be generated automatically
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import *
geo = unit_square
geo.SetMaterial(1,'inner')
normals = geo.CreatePML(1)['normals']
print(normals)
{'pml_bottom': <netgen.libngpy._meshing.Vec2d object at 0x73df04318070>, 'pml_right': <netgen.libngpy._meshing.Vec2d object at 0x73df043180b0>, 'pml_top': <netgen.libngpy._meshing.Vec2d object at 0x73df043181b0>, 'pml_left': <netgen.libngpy._meshing.Vec2d object at 0x73df043180f0>}
The normals in the relevant exterior domains are returned as constant vectors
mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))
print(mesh.GetMaterials())
Draw(mesh)
('inner', 'pml_bottom', 'pml_corner', 'pml_right', 'pml_corner', 'pml_top', 'pml_corner', 'pml_left', 'pml_corner')
BaseWebGuiScene
On this domain we want to approximate the problem to find \(p,\hat p\in H1, v,\hat v\in (L^2)^2\)
\begin{aligned} \partial_t \int_\Omega v\cdot w &= -\int_\Omega \nabla p\cdot w-\alpha\int_{\Omega_{tb}\cup\Omega_{lr}} v\cdot n n^\top w+\alpha\int_{\Omega_{tb}\cup\Omega_{lr}} v\cdot(\mathbf I-n n^\top) w-\alpha\int_{\Omega_{tb}\cup\Omega_{lr}}\hat v\cdot\left(\mathbf I- nn^\top\right) w +\int_{\Omega_{\mathrm{int}}}f w\ \partial_t \int_\Omega pq &= \int_\Omega v\cdot \nabla q-2\alpha\int_{\Omega_c}pq+\alpha\int_{\Omega_c}\hat p q\ \partial_t \int_\Omega \hat p\hat q &=\alpha\int_{\Omega_c}p\hat q\ \partial_t \int_\Omega \hat v\cdot \hat w &=-\alpha\int_{\Omega_{tb}\cup\Omega_{lr}}\hat v\cdot\left(\mathbf I- nn^\top\right)\hat w+\alpha\int_{\Omega_{tb}\cup\Omega_{lr}} v\cdot\left(\mathbf I- nn^\top\right)\hat w \end{aligned} for all \(q,\hat q\in H^1, w,\hat w \in (L^2)^2\), where
We set up the appropriate spaces
order = 2
fes_p = H1(mesh,order = order)
fes_v = VectorL2(mesh,order = order-1)
fes = fes_p*fes_p*fes_v*fes_v
p,phat, v,vhat = fes.TrialFunction()
q,qhat, w,what = fes.TestFunction()
p1,q1 = fes_p.TnT()
v1,w1 = fes_v.TnT()
and mass and gradient matrices
B = BilinearForm(grad(p1)*w1*dx).Assemble()
mass_p = BilinearForm(p1*q1*dx).Assemble()
mass_v = BilinearForm(v1*w1*dx).Assemble()
to set up the damping matrices we need the normal vectors as coefficient function
nvec = { mat : ((normals[mat][0], normals[mat][1]) if mat in normals else (0,0)) for mat in mesh.GetMaterials() }
cfn = CF( [CF(nvec[mat]) for mat in mesh.GetMaterials()])
Draw(cfn,mesh,vectors=True)
BaseWebGuiScene
we define the damping matrices
pml1d = mesh.Materials('pml_left|pml_right|pml_top|pml_bottom')
dampp_lrtb = BilinearForm (p1*q1*dx(definedon=pml1d)).Assemble().mat
dampp_corner = BilinearForm(p1*q1*dx("pml_corner")).Assemble().mat
dampv_n = BilinearForm(v1*(OuterProduct(cfn,cfn)*w1)*dx(definedon=pml1d)).Assemble().mat
dampv_t = BilinearForm(v1*((Id(2)-OuterProduct(cfn,cfn))*w1)*dx(definedon=pml1d)).Assemble().mat
used dof inconsistency
(silence this warning by setting BilinearForm(...check_unused=False) )
used dof inconsistency
(silence this warning by setting BilinearForm(...check_unused=False) )
used dof inconsistency
(silence this warning by setting BilinearForm(...check_unused=False) )
used dof inconsistency
(silence this warning by setting BilinearForm(...check_unused=False) )
To define the big matrices we need a lot of embeddings. Note that we only use vector operations (instead of applying mass matrices) for the auxiliar variables
emb_p, emb_phat, emb_v, emb_vhat = fes.embeddings
fullB = emb_v @ B.mat @ emb_p.T
dampingv = emb_v @ (dampv_n-dampv_t) @ emb_v.T - emb_v @ dampv_t @ emb_vhat.T + emb_vhat @ emb_v.T + emb_vhat @ emb_vhat.T
dampingp = emb_p @ (dampp_lrtb+2*dampp_corner) @ emb_p.T - emb_p @ dampp_corner @ emb_phat.T + emb_phat @ emb_p.T
invmassp = mass_p.mat.Inverse()
invmassv = mass_v.mat.Inverse()
invp = emb_p @ invmassp @ emb_p.T + emb_phat @ emb_phat.T
invv = emb_v @ invmassv @ emb_v.T + emb_vhat @ emb_vhat.T
define a GridFunction and set initial values
gfu = GridFunction(fes)
gfu.vec[:] = 0
gfu.components[0].Set(exp(-100*((x-0.4)**2+(y-0.4)**2)))
scene = Draw(gfu.components[0],autoscale=False,min=0,max=0.5)
tmp = gfu.vec.CreateVector()
tau = 0.001
alpha = 5
tend = 2
t = 0
i = 0
while t < tend:
tmp.data = -fullB.T * gfu.vec
tmp.data -= alpha * dampingp * gfu.vec
gfu.vec.data += tau * invp * tmp
tmp.data = fullB * gfu.vec
tmp.data -= alpha * dampingv * gfu.vec
gfu.vec.data += tau * invv * tmp
t += tau
i += 1
if i%10 == 0:
#print("t = {}".format(t),end = '\r')
scene.Redraw()
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Input In [11], in <cell line: 13>()
25 t += tau
26 i += 1
---> 27 if i%10 == 0:
28 #print("t = {}".format(t),end = '\r')
29 scene.Redraw()
KeyboardInterrupt: