The Helmholtz equation

The Helmholtz equation#

In this tutorial, we will learn:

How to solve PDEs with complex-valued fields,
How to import and use high-order meshes from Gmsh,
How to use high order discretizations,
How to use UFL expressions.

Problem statement#

We will solve the Helmholtz equation subject to a first order absorbing boundary condition:

\[\begin{split} \begin{align*} \Delta u + k^2 u &= 0 && \text{in } \Omega,\\ \nabla u \cdot \mathbf{n} - \mathrm{j}ku &= g && \text{on } \partial\Omega, \end{align*} \end{split}\]

where \(k\) is a piecewise constant wavenumber, \(\mathrm{j}=\sqrt{-1}\), and \(g\) is the boundary source term computed as

\[g = \nabla u_\text{inc} \cdot \mathbf{n} - \mathrm{j}ku_\text{inc}\]

for an incoming plane wave \(u_\text{inc}\).

from mpi4py import MPI

import numpy as np

import dolfinx.fem.petsc
import ufl

This example is designed to be executed with complex-valued degrees of freedom. To be able to solve this problem, we use the complex build of PETSc.

import sys

from petsc4py import PETSc

if not np.issubdtype(PETSc.ScalarType, np.complexfloating):
    print("This tutorial requires complex number support")
    sys.exit(0)
else:
    print(f"Using {PETSc.ScalarType}.")

Using <class 'numpy.complex128'>.

Defining model parameters#

# wavenumber in free space (air)
k0 = 10 * np.pi

# Corresponding wavelength
lmbda = 2 * np.pi / k0

# Polynomial degree
degree = 6

# Mesh order
mesh_order = 2

Interfacing with GMSH#

We will use Gmsh to generate the computational domain (mesh) for this example. As long as Gmsh has been installed (including its Python API), DOLFINx supports direct input of Gmsh models (generated on one process). DOLFINx will then in turn distribute the mesh over all processes in the communicator passed to dolfinx.io.gmshio.model_to_mesh.

The function generate_mesh creates a Gmsh model on rank 0 of MPI.COMM_WORLD.

from dolfinx.io import gmshio
from mesh_generation import generate_mesh

# MPI communicator
comm = MPI.COMM_WORLD

file_name = "domain.msh"
generate_mesh(file_name, lmbda, order=mesh_order)
mesh, cell_tags, _ = gmshio.read_from_msh(file_name, comm, rank=0, gdim=2)

Info    : Reading 'domain.msh'...
Info    : 15 entities
Info    : 2985 nodes
Info    : 1444 elements
Info    : Done reading 'domain.msh'

Material parameters#

In this problem, the wave number in the different parts of the domain depends on cell markers, inputted through cell_tags. We use the fact that a discontinuous Lagrange space of order 0 (cell-wise constants) has a one-to-one mapping with the cells local to the process.

W = dolfinx.fem.functionspace(mesh, ("DG", 0))
k = dolfinx.fem.Function(W)
k.x.array[:] = k0
k.x.array[cell_tags.find(1)] = 3 * k0

import matplotlib.pyplot as plt
import pyvista

from dolfinx.plot import vtk_mesh

pyvista.start_xvfb()
pyvista.set_plot_theme("paraview")
sargs = dict(
    title_font_size=25,
    label_font_size=20,
    fmt="%.2e",
    color="black",
    position_x=0.1,
    position_y=0.8,
    width=0.8,
    height=0.1,
)


def export_function(grid, name, show_mesh=False, tessellate=False):
    grid.set_active_scalars(name)
    plotter = pyvista.Plotter(window_size=(700, 700))
    t_grid = grid.tessellate() if tessellate else grid
    plotter.add_mesh(t_grid, show_edges=False, scalar_bar_args=sargs)
    if show_mesh:
        V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
        grid_mesh = pyvista.UnstructuredGrid(*vtk_mesh(V))
        plotter.add_mesh(grid_mesh, style="wireframe", line_width=0.1, color="k")
        plotter.view_xy()
    plotter.view_xy()
    plotter.camera.zoom(1.3)
    plotter.export_html(f"./{name}.html")

export_function(grid, "wavenumber", show_mesh=True, tessellate=True)

Boundary source term#

\[g = \nabla u_{inc} \cdot \mathbf{n} - \mathrm{j}ku_{inc}\]

where \(u_{inc} = e^{-jkx}\), the incoming wave, is a plane wave propagating in the \(x\) direction.

Next, we define the boundary source term by using ufl.SpatialCoordinate. When using this function, all quantities using this expression will be evaluated at quadrature points.

n = ufl.FacetNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
uinc = ufl.exp(1j * k * x[0])
g = ufl.dot(ufl.grad(uinc), n) - 1j * k * uinc

Variational form#

Next, we define the variational problem using a 4th order Lagrange space. Note that as we are using complex valued functions, we have to use the appropriate inner product; see DOLFINx tutorial: Complex numbers for more information.

\[ -\int_\Omega \nabla u \cdot \nabla \bar{v} ~ dx + \int_\Omega k^2 u \,\bar{v}~ dx - j\int_{\partial \Omega} ku \bar{v} ~ ds = \int_{\partial \Omega} g \, \bar{v}~ ds \qquad \forall v \in \widehat{V}. \]

import basix.ufl

element = basix.ufl.element("Lagrange", mesh.topology.cell_name(), degree)
V = dolfinx.fem.functionspace(mesh, element)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

a = (
    -ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
    + k**2 * ufl.inner(u, v) * ufl.dx
    - 1j * k * ufl.inner(u, v) * ufl.ds
)
L = ufl.inner(g, v) * ufl.ds

Linear solver#

Next, we will solve the problem using a direct solver (LU).

opt = {"ksp_type": "preonly", "pc_type": "lu"}
problem = dolfinx.fem.petsc.LinearProblem(a, L, petsc_options=opt)
uh = problem.solve()
uh.name = "u"

Postprocessing#

Visualising using PyVista#

topology, cells, geometry = vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cells, geometry)
grid.point_data["Abs(u)"] = np.abs(uh.x.array)

Post-processing with Paraview#

from dolfinx.io import VTXWriter

u_abs = dolfinx.fem.Function(V, dtype=np.float64)
u_abs.x.array[:] = np.abs(uh.x.array)

VTXWriter#

# VTX can write higher order functions
with VTXWriter(comm, "out_high_order.bp", [u_abs], engine="BP4") as f:
    f.write(0.0)