OSRC preconditioned high-frequency scattering

Download this tutorial as a Jupyter notebook


A typical application for BEM is the high-frequency scattering of an incoming wave \displaystyle u^\text{inc} from a sound-hard obstacle \displaystyle \Omega.

Consider the exterior Helmholtz problem:

\displaystyle \Delta u^\text{+} + k^2u^\text{+} = 0,\quad\text{in }\Omega^\text{+}\\[2mm] \frac{\partial u^\text{+}}{\partial\nu} = -\frac{\partial u^\text{inc}}{\partial\nu},\quad\text{on }\Gamma\\[2mm] \lim_{|\mathbf{x}|\rightarrow+\infty}|\mathbf{x}|\left(\nabla u^\text{+}\cdot \frac{\mathbf{x}}{|\mathbf{x}|}-\mathrm{i}ku^\text{+}\right)=0

Here, \displaystyle u^+ is the scattered field, and hence the total field is \displaystyle u = u^\text{+}+u^{inc}. The scattered field \displaystyle u^\text{+} can be
represented as

\displaystyle u^{+} = \mathcal{K}\phi,

where \displaystyle \mathcal{K} is the double layer potential operator defined by

\displaystyle \left[\mathcal{K}\phi\right](\mathbf{x}) = \int_{\Gamma}g(\mathbf{x},\mathbf{y})\phi(\mathbf{y})\mathrm{d}\mathbf{y}

with \displaystyle g(\mathbf{x},\mathbf{y}) = \frac{\mathrm{e}^{\mathrm{i}k|\mathbf{x}-\mathbf{y}|}}{4\pi|\mathbf{x}-\mathbf{y}|}. Here, \displaystyle \phi=\gamma_0u is the Dirichlet trace of the total field.

The Burton–Miller formulation to compute \displaystyle \phi is given by

\displaystyle \left(\tfrac{1}{2}\mathsf{Id}-\mathsf{K}+\eta\mathsf{W}\right)\phi = \gamma_0 u^\text{inc}+\eta \gamma_1u^\text{inc}

for some \displaystyle \eta\neq 0. Here, \displaystyle \mathsf{Id}, \displaystyle \mathsf{K} and \displaystyle \mathsf{W} are the identity, double layer and hypersingular boundary operators.

A perfect choice for \displaystyle \eta is the exterior Neumann-to-Dirichlet map \displaystyle \mathsf{NtD} with the property

\displaystyle \tfrac{1}{2}\mathsf{Id}-\mathsf{K}-\mathsf{NtD}\circ\mathsf{W} = \mathsf{Id}.

However, this is a complicated non-local operator. The idea of OSRC is to use the pseudodifferential operator approximation

\displaystyle \mathsf{NtD}\approx \frac{1}{\mathrm{i}k}\left(1+\frac{\Delta_{\Gamma}}{(k+\mathrm{i}\epsilon)^2}\right)^{-\frac{1}{2}}

for a regularization parameter \displaystyle \epsilon>0 and to localize the squareroot operator by a Padé approximation. Details of this OSRC preconditioning are given in Antoine & Darbase (2007).

Bempp implements approximate Dirichlet-to-Neumann and Neumann-to-Dirichlet maps. This notebook demonstrates how to use these for high-frequency scattering computations.


In the following we implement an OSRC Burton–Miller formulation for high-frequency scattering. We start with the usual imports.

import bempp.api
import numpy as np

The following defines the wavenumber k and the Dirichlet and Neumann data of the incident plane wave.

= 5

def dirichlet_fun(x, n, domain_index, result):
result[0] = np.exp(1j * k * x[0])

def neumann_fun(x, n, domain_index, result):
result[0] = 1j * k * n[0] * np.exp(1j * k * x[0])

For this example we will use an elongated ellipsoid. The element size is chosen to roughly correspond to 10 elements per wavelength. The function space consists of continuous, piecewise linear basis functions. The OSRC preconditioner only works for continuous function spaces as it assembles a Laplace–Beltrami operator.

= 2*np.pi/(10 * k)
grid = bempp.api.shapes.ellipsoid(3, 1, 1, h=h)
space = bempp.api.function_space(grid, “P”, 1)
print(“The space has {0} dofs”.format(space.global_dof_count))
The space has 2477 dofs

We can now form the Burton–Miller operator.

identity = bempp.api.operators.boundary.sparse.identity(
space, space, space)
dlp = bempp.api.operators.boundary.helmholtz.double_layer(
space, space, space, k)
hyp = bempp.api.operators.boundary.helmholtz.hypersingular(
space, space, space, k, use_slp=True)
ntd = bempp.api.operators.boundary.helmholtz.osrc_ntd(space, k)

burton_miller = .5 * identity  dlp  ntd * hyp

In the above operator assembly we were allowed to directly multiply operators. Bempp automatically assembles the correct mass matrix transformations for this product to make sense. For the assembly of the hypersingular operator we used the parameter use_slp=True. This will cause the hypersingular operatorto be assembled in terms of an underlying Helmholtz single layer operator. This approach takes more memory but the H-matrix assembly is significantly faster.

We next assemble the right-hand side.

dirichlet_grid_fun = bempp.api.GridFunction(space, fun=dirichlet_fun)
neumann_grid_fun = bempp.api.GridFunction(space, fun=neumann_fun)
rhs_fun = dirichlet_grid_fun  ntd * neumann_grid_fun

We can now solve the Burton–Miller formulation using GMRES. We use a strong form discretisation: this automatically performs mass matrix preconditioning.

total_field, info, it_count = bempp.api.linalg.gmres(
burton_miller, rhs_fun, use_strong_form=True,

print(“The linear system was solved in {0} iterations”.format(it_count))

The linear system was solved in 6 iterations

We now want to plot the radar cross section in the \displaystyle z=0 plane. To compute it we use the far-field operators implemented in Bempp.

from matplotlib import pylab as plt

theta = np.linspace(np.pi/2, np.pi/2, 400)
points = np.array([np.cos(theta), np.sin(theta), np.zeros(len(theta))])

dlp_far_field = bempp.api.operators.far_field.helmholtz.double_layer(
space, points, k)
far_field = dlp_far_field * total_field
max_incident = np.abs(dirichlet_grid_fun.coefficients).max()
radiation_pattern = (np.abs(far_field/max_incident)**2).ravel()
db_pattern = 10 * np.log10(4 * np.pi * radiation_pattern)

plt.polar(theta, db_pattern)
plt.ylim(db_pattern.min(), db_pattern.max())
plt.title(‘RCS (DB)’)