3.3. Complex-valued waves

In NGSolve finite element spaces can be built over the complex field and the resulting complex matrix systems can be solved. This tutorial shows how to compute the solution of the Helmholtz equation with impedance boundary conditions in complex arithmetic. The boundary value problem is to find \(u\) satisfying

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

together with the impedance (outgoing) boundary condition

\[ \frac{\partial u }{ \partial n} - i \omega u = 0 \quad \text{ on } \partial \Omega \]

where \(i =\) 1j is the imaginary unit.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *

air: Face = Circle((0.5, 0.5), 0.8).Face()
air.edges.name = 'outer'
scatterer = MoveTo(0.7, 0.3).Rectangle(0.05, 0.4).Face()
scatterer.edges.name = 'scat'
geo = OCCGeometry(air - scatterer, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.05))
Draw(mesh);

3.3.1. Declare a complex-valued finite element space

fes = H1(mesh, order=5, complex=True)
u, v = fes.TnT()

omega = 10
pulse = 5e4*exp(-(40**2)*((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5)))
Draw(10**-5*pulse, mesh, order=3);

3.3.2. Forming the system

The weak form for \(u \in H^1\):

\[ \int_\Omega\big[ \nabla u \cdot \nabla \bar v - \omega^2 u \bar v \big] \, dx - i \,\omega\, \int_{\partial \Omega} u \bar v \, ds = \int_{\Omega} f \bar v \]

for all \(v\) in \(H^1\).

# Forms

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


f = LinearForm(pulse * v * dx).Assemble();

3.3.3. Solve

gfu = GridFunction(fes, name="u")
gfu.vec.data = a.mat.Inverse() * f.vec
Draw(gfu, mesh, min=-10, max=10, order=3, animate_complex=True);

Open controls and explore webgui’s further viewing options such as viewing real and imaginary parts, viewing absolute value, and viewing with deformation turned on.

3.3.4. Add a PML

Consider the same problem as before

\[ -\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 \]

define now an outer region called for the PML

outer = Circle((0.5, 0.5), 1).Face()
outer.edges.name = 'outer'
air: Face = Circle((0.5, 0.5), 0.8).Face()
#air.edges.name = 'outer'
scatterer = MoveTo(0.7, 0.3).Rectangle(0.05, 0.4).Face()
scatterer.edges.name = 'scat'

pmlregion = outer - air
pmlregion.faces.name = 'pmlregion'
geo = OCCGeometry(Glue([air-scatterer, pmlregion]), dim=2)

mesh = Mesh(geo.GenerateMesh(maxh=0.05))
Draw(mesh);

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=0.8, alpha=1j,origin=(0.5,0.5)), "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 Unit 1.7:

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


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

f = LinearForm(pulse * v * dx).Assemble()

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

settings = {}
Draw(gfu, animate_complex=True, deformation=False, scale = 0.01);