Simple monopole source example gets a wrong result

Applications
1

I am learning to use bempp to solve acoustic problems. I test it with a simple monopole source example:

p(x)=\frac{e^{-i k r}}{4 \pi r}

I use the Fredholm integral equation of the first kind to calculate Dirichlet data according to Neumann data:

-\frac{1}{2} g_{\mathrm{D}}(x)+\left(K_{\kappa} g_{\mathrm{D}}\right)(x) = \left(V_{\kappa} g_{\mathrm{N}}\right)(x)\quad \text { for } x \in \partial \Omega

where Neumann data g_{\mathrm{N}}:=\gamma^{1, \mathrm{ext}} u and Dirichlet data g_{\mathrm{D}}:=\gamma^{0, \mathrm{ext}} u. K_{\kappa} and V_{\kappa} represents Double Layer Potential Operator and Single Layer Potential Operator.

I expect that the Dirichlet data calculated by the integral is equal to the monopole Dirichlet data, but the result is different. I can’t understand why this happens. Can anyone point out the problem? I have been troubled for a long time.

Thank you!

I installed the latest bempp from github according to the documentation. Here is my test code:

import bempp.api
import numpy as np
grid = bempp.api.shapes.sphere(h = 0.1)

dp0_space = bempp.api.function_space(grid, "DP", 0)
p1_space = bempp.api.function_space(grid, "P", 1)
k = 10

@bempp.api.complex_callable()
def monopole_neumann_data(x, n, domain_index, result):
    r = (x[0]**2+x[1]**2+x[2]**2)**0.5
    result[0] = (-1j*k - 1/r)*(np.exp(-1j*k*r))/(4*np.pi*r)

@bempp.api.complex_callable()
def monopole_dirichlet_data(x, n, domain_index, result):
    r = (x[0]**2+x[1]**2+x[2]**2)**0.5
    result[0] = (np.exp(-1j*k*r))/(4*np.pi*r)

monopole_neumann_fun = bempp.api.GridFunction(dp0_space, fun=monopole_neumann_data)
monopole_dirichlet_fun = bempp.api.GridFunction(p1_space, fun=monopole_dirichlet_data)

identity = bempp.api.operators.boundary.sparse.identity(
    p1_space, p1_space, p1_space)
double_layer = bempp.api.operators.boundary.helmholtz.double_layer(
    p1_space, p1_space, p1_space, k)
single_layer = bempp.api.operators.boundary.helmholtz.single_layer(
    dp0_space, p1_space, p1_space, k)

left_side = -.5*identity+double_layer
right_side = single_layer*monopole_neumann_fun

output_dirichlet_fun, info = bempp.api.linalg.gmres(left_side, right_side,tol=1e-5)
print(monopole_dirichlet_fun.coefficients[0])
print(output_dirichlet_fun.coefficients[0])

I get the result:

(-0.06738962350437123+0.042509072890037296j)
(0.05712675186005465-0.05505303182725231j)

PS: If I change the formulation to p(x)=\frac{e^{i k r}}{4 \pi r}, I get a result agree with the original data. But this is not reasonable in theory.

2

Hi. The reason is the last sentence you mentioned. It works for e^{ikr}{4\pi r} because this is what we use as Green’s function. Your monopole definition is in our convention an incoming wave. There is nothing wrong with using either a positve or a negative sign in the definition of the outgoing wave. It just needs to be consistent. Most mathematical literature uses our definition. Some engineering literature uses the minus. The difference comes in when you want do a transformation to a time domain. Then you have to take the choice of sign into account.

3

Thank you very much! So if I change the definition, the Hankel function of the second kind should be changed to the Hankel function of the first kind. Similarly, in other correlation formulas, I only need to change the i to -i. Is that right？

4

Hi. You just need to do what you already proposed, namely instead of e^{ik|x-y|} use e^{-ik|x-y|}. Hankel functions are the Green’s functions in 2 dimensions.

Best wishes

Timo