Helmholtz equation

Copyright (C) 2018 Samuel Groth

Helmholtz problem in both complex and real modes 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 Method of Manufactured Solutions is used to produce the exact solution and source term.

import numpy as np
from dolfinx import (Function, FunctionSpace, UnitSquareMesh, fem,
                     has_petsc_complex)
from dolfinx.fem.assemble import assemble_scalar
from dolfinx.io import XDMFFile
from mpi4py import MPI
from ufl import FacetNormal, TestFunction, TrialFunction, dx, grad, inner

# wavenumber
k0 = 4 * np.pi

# approximation space polynomial degree
deg = 1

# number of elements in each direction of mesh
n_elem = 128

mesh = UnitSquareMesh(MPI.COMM_WORLD, n_elem, n_elem)
n = FacetNormal(mesh)

# Source amplitude
if has_petsc_complex:
    A = 1 + 1j
else:
    A = 1


# Test and trial function space
V = FunctionSpace(mesh, ("Lagrange", deg))

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V)
f.interpolate(lambda x: A * k0**2 * np.cos(k0 * x[0]) * np.cos(k0 * x[1]))
a = inner(grad(u), grad(v)) * dx - k0**2 * inner(u, v) * dx
L = inner(f, v) * dx

# Compute solution
uh = fem.Function(V)
uh.name = "u"
problem = fem.LinearProblem(a, L, u=uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
problem.solve()

# Save solution in XDMF format (to be viewed in Paraview, for example)
with XDMFFile(MPI.COMM_WORLD, "plane_wave.xdmf", "w", encoding=XDMFFile.Encoding.HDF5) as file:
    file.write_mesh(mesh)
    file.write_function(uh)

Calculate L2 and H1 errors of FEM solution and best approximation. This demonstrates the error bounds given in Ihlenburg. Pollution errors are evident for high wavenumbers.

# Function space for exact solution - need it to be higher than deg
V_exact = FunctionSpace(mesh, ("Lagrange", deg + 3))
u_exact = Function(V_exact)
u_exact.interpolate(lambda x: A * np.cos(k0 * x[0]) * np.cos(k0 * x[1]))

# H1 errors
diff = uh - u_exact
H1_diff = mesh.mpi_comm().allreduce(assemble_scalar(inner(grad(diff), grad(diff)) * dx), op=MPI.SUM)
H1_exact = mesh.mpi_comm().allreduce(assemble_scalar(inner(grad(u_exact), grad(u_exact)) * dx), op=MPI.SUM)
print("Relative H1 error of FEM solution:", abs(np.sqrt(H1_diff) / np.sqrt(H1_exact)))

# L2 errors
L2_diff = mesh.mpi_comm().allreduce(assemble_scalar(inner(diff, diff) * dx), op=MPI.SUM)
L2_exact = mesh.mpi_comm().allreduce(assemble_scalar(inner(u_exact, u_exact) * dx), op=MPI.SUM)
print("Relative L2 error of FEM solution:", abs(np.sqrt(L2_diff) / np.sqrt(L2_exact)))