Burton-Miller formulation is not correct

Applications
1

I found that the Burton-Miller formulation yields higher errors compared to the standard boundary element equation in my test case. Can someone explain why this might be happening?

Reproduction code:

import bempp.api
import numpy as np

bempp.api.enable_console_logging("debug")
bempp.api.BOUNDARY_OPERATOR_DEVICE_TYPE = "gpu"
bempp.api.POTENTIAL_OPERATOR_DEVICE_TYPE = "gpu"

x0 = np.array([0.0, 0.0, 0.0])
rerr_lst = []
wave_number = 3


@bempp.api.complex_callable
def get_dirichlet(x, n, domain_index, result):
    r_vec = x - x0
    r = np.dot(r_vec, r_vec) ** 0.5
    k = wave_number
    result[0] = np.exp(1j * k * r) / (4 * np.pi * r)


@bempp.api.complex_callable
def get_neumann(x, n, domain_index, result):
    r_vec = x - x0
    r = np.dot(r_vec, r_vec) ** 0.5
    k = wave_number
    ikr = 1j * k * r
    result[0] = (
        -np.exp(ikr) / (4 * np.pi * r * r * r) * (1 - ikr) * ((x - x0) * n).sum()
    )


wavelength = 2 * np.pi / abs(wave_number)
h = wavelength / 6
grid = bempp.api.shapes.ellipsoid(3, 1, 1, h=h * 0.5)

space = bempp.api.function_space(grid, "P", 1)

identity = bempp.api.operators.boundary.sparse.identity(space, space, space)
slp = bempp.api.operators.boundary.helmholtz.single_layer(
    space, space, space, wave_number
)
dlp = bempp.api.operators.boundary.helmholtz.double_layer(
    space, space, space, wave_number
)
hyp = bempp.api.operators.boundary.helmholtz.hypersingular(
    space, space, space, wave_number
)
adlp = bempp.api.operators.boundary.helmholtz.adjoint_double_layer(
    space, space, space, wave_number
)

LHS = 0.5 * identity - dlp
RHS = -slp

neumann_fun = bempp.api.GridFunction(space, fun=get_neumann)

RHS = RHS * neumann_fun

dirichlet_fun, info = bempp.api.linalg.gmres(LHS, RHS, tol=1e-5)
dirichlet_fun_gt = bempp.api.GridFunction(space, fun=get_dirichlet)

rerr = np.linalg.norm(
    dirichlet_fun.coefficients - dirichlet_fun_gt.coefficients
) / np.linalg.norm(dirichlet_fun_gt.coefficients)

rerr_lst.append(rerr)

beta = 1j / wave_number
LHS = 0.5 * identity - dlp - beta * hyp
RHS = (-slp - beta * (adlp + 0.5 * identity)) * neumann_fun

dirichlet_fun, info = bempp.api.linalg.gmres(LHS, RHS, tol=1e-5)
dirichlet_fun_gt = bempp.api.GridFunction(space, fun=get_dirichlet)


rerr = np.linalg.norm(
    dirichlet_fun.coefficients - dirichlet_fun_gt.coefficients
) / np.linalg.norm(dirichlet_fun_gt.coefficients)

rerr_lst.append(rerr)
print(rerr_lst)

The output is

[0.00011713474983450146, 1.5905369315737763]

The error is too high, I even wonder if the Bempp implementation is wrong?

2

I found the problem. There is a negative sign in front of the hypersingular operator defined in Bempp, so negating the hypersingular operator solves the issue.