# Helmholtz equation

Copyright (C) 2018-2025 Samuel Groth and Jørgen S. Dokken

Download sources

- `Python script`

- `Jupyter notebook`

This demo illustrates how to:

- Create a complex-valued finite element formulation In the following example, we will consider the Helmholtz equation solved with both a complex valued and a real valued finite element formulation.

In the complex mode, the exact solution is a plane wave propagating at an angle theta to the positive x-axis. Chosen for comparison with results from Ihlenburg’s book Finite Element Analysis of Acoustic Scattering, p138-139. In real mode, the exact solution corresponds to the real part of the plane wave (a sin function which also solves the homogeneous Helmholtz equation).

We start by importing the necessary modules

    from mpi4py import MPI
    from petsc4py import PETSc

    import numpy as np

    import ufl
    from dolfinx.fem import (
        Expression,
        Function,
        assemble_scalar,
        dirichletbc,
        form,
        functionspace,
        locate_dofs_geometrical,
    )
    from dolfinx.fem.petsc import LinearProblem
    from dolfinx.io import XDMFFile
    from dolfinx.mesh import create_unit_square

We define the necessary parameters for the discretized problem

    k0 = 4 * np.pi  # Wavenumber
    deg = 1  # Approximation space polynomial degree
    n_elem = 64  # Number of elements in each direction of the mesh
    A = 1  # Source amplitude

Next, we create the discrete domain, a unit square and set up the discrete function space.

    msh = create_unit_square(MPI.COMM_WORLD, n_elem, n_elem)
    V = functionspace(msh, ("Lagrange", deg))

## Define variational problem.

The Helmholtz equation can be discretized in the same way for both the real and complex valued formulation. However, note that we use `ufl.inner` instead of `ufl.dot` or the `*` operator between the test and trial function, and that the test-function is always the second variable in the operator. The reason for this is that for complex variational forms, one requires a sesquilinear two-form, with the inner product being \((a,b)=\int_\Omega a \cdot \bar{b}~\mathrm{d}x\).

    u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
    a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx - k0**2 * ufl.inner(u, v) * ufl.dx

We solve for plane wave with mixed Dirichlet and Neumann BCs. We use ufl to manufacture an exact solution and corresponding boundary conditions.

    theta = np.pi / 4
    V_exact = functionspace(
        msh, ("Lagrange", deg + 3)
    )  # Exact solution should be in a higher order space
    u_exact = Function(V_exact, name="u_exact")
    u_exact.interpolate(lambda x: A * np.exp(1j * k0 * (np.cos(theta) * x[0] + np.sin(theta) * x[1])))
    x = ufl.SpatialCoordinate(msh)
    n = ufl.FacetNormal(msh)
    g = -ufl.dot(n, ufl.grad(u_exact))
    L = -ufl.inner(g, v) * ufl.ds

    dofs_D = locate_dofs_geometrical(
        V, lambda x: np.logical_or(np.isclose(x[0], 0), np.isclose(x[1], 0))
    )
    u_bc = Function(V)
    u_bc.interpolate(Expression(u_exact, V.element.interpolation_points))
    bcs = [dirichletbc(u_bc, dofs_D)]

In this problem, we rely on PETSc as the linear algebra backend. PETSc can only be configured for a either real of complex valued matrices and vectors. We can check how PETSc is configured by calling

    is_complex_mode = np.issubdtype(PETSc.ScalarType, np.complexfloating)
    PETSc.Sys.Print(f"PETSc is configured in complex mode: {is_complex_mode}")

    uh = Function(V)
    uh.name = "u"
    problem = LinearProblem(
        a,
        L,
        bcs=bcs,
        u=uh,
        petsc_options_prefix="demo_helmholtz_",
        petsc_options={"ksp_type": "preonly", "pc_type": "lu", "ksp_error_if_not_converged": True},
    )
    _ = problem.solve()

Save solution in XDMF format (to be viewed in ParaView, for example)

    with XDMFFile(
        MPI.COMM_WORLD, "out_helmholtz/plane_wave.xdmf", "w", encoding=XDMFFile.Encoding.HDF5
    ) as file:
        file.write_mesh(msh)
        file.write_function(uh)

Calculate \(L_2\) and \(H_0^1\) errors of FEM solution and best approximation. This demonstrates the error bounds given in Ihlenburg. Pollution errors are evident for high wavenumbers.

    diff = uh - u_exact
    H1_diff = msh.comm.allreduce(
        assemble_scalar(form(ufl.inner(ufl.grad(diff), ufl.grad(diff)) * ufl.dx)), op=MPI.SUM
    )
    H1_exact = msh.comm.allreduce(
        assemble_scalar(form(ufl.inner(ufl.grad(u_exact), ufl.grad(u_exact)) * ufl.dx)), op=MPI.SUM
    )
    PETSc.Sys.Print("Relative H1 error of FEM solution:", abs(np.sqrt(H1_diff) / np.sqrt(H1_exact)))

    L2_diff = msh.comm.allreduce(assemble_scalar(form(ufl.inner(diff, diff) * ufl.dx)), op=MPI.SUM)
    L2_exact = msh.comm.allreduce(
        assemble_scalar(form(ufl.inner(u_exact, u_exact) * ufl.dx)), op=MPI.SUM
    )
    PETSc.Sys.Print("Relative L2 error of FEM solution:", abs(np.sqrt(L2_diff) / np.sqrt(L2_exact)))