Partial differential equations » Solving the Helmholtz equation

Resolution of the complex-valued Helmholtz equation using Rodin.

Image

Introduction

The Helmholtz equation models time-harmonic wave phenomena in acoustics, electromagnetics, and vibration analysis. In this example we solve the complex-valued Helmholtz equation with Dirichlet boundary conditions:

\[ \left\{ \begin{aligned} -\Delta u - k^2 u &= f && \mathrm{in} \ \Omega \\ u &= u_{\mathrm{ex}} && \mathrm{on} \ \partial \Omega \end{aligned} \right. \]

where $ k = 2\pi / \lambda $ is the wave number, $ \lambda $ is the wavelength, $ f $ is a known source term, and $ u_{\mathrm{ex}} $ is a known exact solution imposed on the boundary.

This example demonstrates several important Rodin features:

  • Complex-valued finite element spaces (P1<Complex>)
  • Extracting real and imaginary parts of complex GridFunction objects
  • The Helmholtz bilinear form combining stiffness and mass terms

Weak Formulation

The weak formulation of the Helmholtz equation is:

\[ \text{Find } u \in H^1(\Omega) \text{ s.t. } \forall v \in H^1_0(\Omega), \quad \int_{\Omega} \nabla u \cdot \nabla v \ dx - k^2 \int_{\Omega} u \, v \ dx = \int_{\Omega} f \, v \ dx \]

Note the key difference from the Poisson equation: the presence of the mass term $ -k^2 \int u v \, dx $ , which makes the bilinear form indefinite for sufficiently large $ k $ .

Implementation

Includes and Namespaces

We include the standard Rodin modules plus the math constants for $ \pi $ :

#include <Rodin/Solver.h>
#include <Rodin/Geometry.h>
#include <Rodin/Assembly.h>
#include <Rodin/Variational.h>

using namespace Rodin;
using namespace Rodin::Geometry;
using namespace Rodin::Variational;

Mesh and Parameters

We use a fine mesh ( $ 64 \times 64 $ ) to adequately resolve the wave, following the rule of thumb that the mesh size $ h $ should satisfy $ k h \lesssim 1 $ . The mesh is created using Mesh::UniformGrid() and scaled to the unit domain via Mesh::scale():

static constexpr Real waveLength = 0.5;
const Real waveNumber = 2 * Math::Constants::pi() / waveLength;

Mesh mesh;
mesh = mesh.UniformGrid(Polytope::Type::Triangle, { 64, 64 });
mesh.getConnectivity().compute(1, 2);
mesh.scale(1.0 / 63.0);

Complex Finite Element Space

The Helmholtz equation naturally involves complex-valued fields. Rodin supports this through template-parameterized finite element spaces. Using P1<Complex> creates a space where each DOF is a complex number:

P1<Complex> vh(mesh);

TrialFunction u(vh);
TestFunction  v(vh);

Manufactured Exact Solution

For verification we use a manufactured exact solution and compute the corresponding right-hand side $ f $ :

// Exact solution: u_ex(x,y) = exp(i k (x cos(theta) + y sin(theta)))
ComplexFunction u_ex =
  [&](const Geometry::Point& p) -> Complex
  {
    const Real x = p.x(), y = p.y();
    return std::exp(Complex(0, waveNumber * (x * cosTheta + y * sinTheta)));
  };

// Right-hand side computed from u_ex
ComplexFunction f =
  [&](const Geometry::Point& p) -> Complex
  {
    const Real x = p.x(), y = p.y();
    return /* ... manufactured source term ... */;
  };

Problem Definition

The Problem is assembled using the same form language as for real-valued PDEs. The key terms are:

Problem helmholtz(u, v);
helmholtz = Integral(Grad(u), Grad(v))                   // stiffness
          - waveNumber * waveNumber * Integral(u, v)      // mass (shift)
          - Integral(f, v)                                // load
          + DirichletBC(u, u_ex);                         // boundary data

Here:

  • Integral(Grad(u), Grad(v)) is the stiffness term $ \int_\Omega \nabla u \cdot \nabla v \, dx $
  • waveNumber * waveNumber * Integral(u, v) is the scaled mass term $ k^2 \int_\Omega u v \, dx $
  • Integral(f, v) is the right-hand side
  • DirichletBC(u, u_ex) imposes the exact solution on the boundary

Solving the System

Because the system is complex-valued and symmetric but not positive-definite (the bilinear form is indefinite), we use a direct solver:

Solver::SparseLU(helmholtz).solve();

Extracting Real and Imaginary Parts

Since standard visualization tools typically work with real-valued data, we extract the real and imaginary parts of the solution into separate real-valued GridFunction objects using Re() and Im():

P1<Real> rh(mesh);
GridFunction uRe(rh), uIm(rh);

uRe = Re(u.getSolution());
uIm = Im(u.getSolution());

uRe.save("uRe.gf", IO::FileFormat::MFEM);
uIm.save("uIm.gf", IO::FileFormat::MFEM);

Full Source Code

The complete source code is:

See Also

« Solving the elasticity equation | Partial differential equations | Solving the Navier-Stokes equations »