This page was generated from unit-1.7-helmholtz/pml.ipynb.

1.7.1 Perfectly Matched Layer (PML)

The use of (Perfectly Matched Layer) PML is a standard technique to numerically solve for outgoing waves in unbounded domains. Although scattering problems are posed in unbounded domains, by bounding the scatterer and any inhomogeneities within a PML, one is able to truncate to a bounded computational domain. In this tutorial we see how to use PML for - source problems, and - eigenvalue problems (resonances).

import netgen.gui
from netgen.geom2d import SplineGeometry
from ngsolve import *
import matplotlib.pyplot as plt

Source problem

Consider the problem of finding \(u\) satisfying

\[-\Delta u - \omega^2 u = f \qquad \text{ in } \mathbb{R}^2\]

together with the Sommerfeld (outgoing) radiation condition at infinity

\[\lim_{r \to \infty} r^{1/2} \bigg( \frac{\partial u }{ \partial r} - i \omega u \bigg) = 0\]

where \(r\) is the radial coordinate.

We consider an example where \(f\) represents a time-harmonic pulse that is essentially zero except for a small region.

geo = SplineGeometry()
geo.AddCircle( (0,0), 1.4, leftdomain=2, bc="outerbnd")
geo.AddCircle( (0,0), 1, leftdomain=1, rightdomain=2, bc="innerbnd")
geo.SetMaterial(1, "inner")
geo.SetMaterial(2, "pmlregion")
mesh = Mesh(geo.GenerateMesh (maxh=0.1))
mesh.Curve(3)

f = exp(-20**2*((x-0.3)*(x-0.3)+y*y))
Draw(f, mesh, 'f')

Note how geometry is divided into an inner disk and and an outer annulus called pmlregion.

The PML facility in NGSolve implements a complex coordinate transformation on a given mesh region (which in this example is pmlregion). When this complex variable change is applied to the outgoing solution in the PML region, it becomes a a function that decays exponentially as \(r \to \infty\), allowing us to truncate the unbounded domain.

With the following single line, we tell NGSolve to turn on this coordinate transformation.

mesh.SetPML(pml.Radial(rad=1,alpha=1j,origin=(0,0)), "pmlregion")

Then a radial PML is set in the exterior of a disk centered at origin of radius rad. In addition to origin and rad, the keyword argument alpha may be used to set the PML strength, which determines the rate of increase in the imaginary part of the coordinate map as radius increases.

Having set the PML, the rest of the code now looks very much like that in 1.7:

fes = H1(mesh, order=4, complex=True)
u = fes.TrialFunction()
v = fes.TestFunction()

omega = 10

a = BilinearForm(fes)
a += grad(u)*grad(v)*dx - omega**2*u*v*dx
a += -1j*omega*u*v*ds("outerbnd")

b = LinearForm(fes)
b += f * v * dx

a.Assemble()
b.Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse() * b.vec
Draw(gfu)

Note that above we have kept the same Robin boundary conditions as in 1.7. (Since the solution exponentially decays within PML, it is also common practice to put zero Dirichlet boundary conditions at the outer edge of the PML, an option that is also easily implementable in NGSolve using the dirichlet flag.)

Eigenvalue problem

PML can also be used to compute resonances on unbounded domains. Resonances of a scatterer \(S\) are outgoing waves \(u\) satisfying the homogeneous Helmholtz equation

\[\Delta u + \omega^2 u = 0 \qquad \text{ in } \mathbb{R}^2 \setminus S,\]

i.e., \(u\) is an eigenfunction of \(-\Delta\) associated to eigenvalue \(\lambda :=\omega^2.\)

We truncate this unbounded-domain eigenproblem to a bounded-domain eigenproblem using PML. The following example computes the resonances of a cavity opening on one side into free space.

geo   = SplineGeometry()
pnums = [geo.AddPoint(x, y, maxh=maxh)
         for x,y,maxh in
         [(-2, 0, 1), (-1, 0, 1), (-0.2, 0, 0.01), (-0.2, -0.8, 1), (0.2, -0.8, 1), (0.2, 0, 0.01), (1,0,1), (2,0,1), (2,2,1), (0,2,1), (-2,2,1), (1,1,1), (0,1,1), (-1,1,1) ]
        ]
lines = [(0,1,"dir",2,0), (1,2,"dir", 1,0), (2,3,"dir",1,0), (3,4,"dir",1,0), (4,5,"dir",1,0), (5,6,"dir",1,0), (6,7,"dir",2,0)]
curves= [(7,8,9,"outer",2,0), (9,10,0,"outer",2,0), (6,11,12,"inner",1,2), (12,13,1,"inner",1,2) ]

for p1,p2,bc,left,right in lines:
    geo.Append(["line", pnums[p1], pnums[p2]], bc=bc, leftdomain=left, rightdomain=right)
for p1,p2,p3,bc,left,right in curves:
    geo.Append(["spline3", pnums[p1], pnums[p2], pnums[p3]], bc=bc, leftdomain=left, rightdomain=right)

geo.SetMaterial(1, "air")
geo.SetMaterial(2, "pml")

mesh = Mesh(geo.GenerateMesh(maxh=0.1))
mesh.Curve(5)
Draw(mesh)

The top edge of the rectangular cavity opens into air, which is included in the computational domain as a semicircular region (neglecting backward propagating waves), while the remaining edges of the cavity are isolated from rest of space by perfect reflecting boundary conditions.

You can identify the PML region in the Netgen GUI by double clicking on the outer part of the air.

mesh.SetPML(pml.Radial(rad=1, alpha=1j, origin=(0,0)), "pml")

We now make the complex symmetric system with PML incorporated.

fes = H1(mesh, order=4, complex=True, dirichlet="dir")

u = fes.TrialFunction()
v = fes.TestFunction()

a = BilinearForm (fes, symmetric=True)
a += grad(u)*grad(v)*dx

m = BilinearForm (fes, symmetric=True)
m += u*v*dx

a.Assemble()
m.Assemble()

For solving the eigenproblem, we use an Arnoldi eigensolver, to which a GridFunction with many vectors (allocated via keyword argument multidim) is input, together with a and b. A shift argument to ArnoldiSolver indicates the approximate location of expected eigenvalues.

u = GridFunction(fes, multidim=50, name='resonances')

with TaskManager():
    lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(),
                        u.vecs, shift=400)
Draw(u)

print ("lam: ", lam)

lam:   (362.23,-28.5602)
 (454.049,-51.7395)
 (297.682,-12.5071)
 (259.413,-3.09602)
 (386.691,-125.274)
 (358.628,-125.478)
 (568.415,-2.47227)
 (472.258,-135.782)
 (255.401,-73.613)
 (608.213,-9.97169)
 (575.05,-82.9482)
 (542.504,-133.232)
 (426.086,-184.537)
 (268.798,-150.773)
 (674.984,-22.8565)
 (552.75,-164.806)
 (314.768,-194.426)
 (167.983,-40.878)
 (623.255,-136.703)
 (368.114,-224.738)
 (305.068,-230.26)
 (489.467,-224.234)
 (353.99,-238.438)
 (290.999,-234.224)
 (312.734,-242.057)
 (329.786,-243.983)
 (343.038,-247.339)
 (261.546,-217.501)
 (421.729,-239.072)
 (270.84,-232.624)
 (387.479,-249.659)
 (367.419,-254.404)
 (250.621,-218.838)
 (405.149,-247.777)
 (217.62,-195.691)
 (234.068,-215.784)
 (181.44,-157.408)
 (408.408,-255.361)
 (427.924,-257.003)
 (213.52,-202.725)
 (107.767,-18.0416)
 (479.072,-245.929)
 (461.101,-254.232)
 (447.994,-258.94)
 (194.706,-194.463)
 (174.014,-184.743)
 (160.944,-160.08)
 (769.324,-41.2526)
 (510.609,-238.802)
 (491.857,-263.065)

The more "confined" modes have resonance values \(\omega\) that are closer to the real axis. Here are the computed \(\omega\)-values plotted in the complex plane.

lamr = [sqrt(l).real for l in lam]
lami = [sqrt(l).imag for l in lam]
plt.plot(lamr, lami, ".")
plt.grid(True)
plt.title('Computed $\omega$ values')
plt.show()