Inconsistent results for different Helmholtz potential assemblers

[antbonomo]

Using bempp-cl 0.2.0, I obtain different answers for the radiated sound pressure level depending on the assembler used (FMM, OpenCL, or Numba). Below is some code that illustrates the issue. Only the Numba-assembled potentials give pressure levels in good agreement with the analytical solution.

# %% [markdown]
# # Pulsating Sphere (simple direct formulation)
# 
# Compare potential assemblers supported by bempp-cl version 0.2.0
# %% [markdown]
# ## Import Modules

# %%
import bempp.api
import bempp.core
import numpy as np
import pandas as pd
bempp.core.opencl_kernels.show_available_platforms_and_devices()
bempp.core.opencl_kernels.default_device()

# %% [markdown]
# ## Define Wavenumber and Acoustic Properties

# %%
k = 1
c = 343
a = 1
rho = 1.21
Wdot = 1
omega = c*k

# %% [markdown]
# ## Mesh Unit Sphere

# %%
grid = bempp.api.shapes.sphere(h=0.1)

# %% [markdown]
# ## Define piecewise constant basis functions

# %%
piecewise_const_space = bempp.api.function_space(grid, "DP", 0)

# %% [markdown]
# ## Define velocity boundary condition

# %%
@bempp.api.complex_callable
def f(x, n, domain_index, result):
    result[0] = Wdot

vn = bempp.api.GridFunction(piecewise_const_space, fun=f)

# %% [markdown]
# ## Define boundary operators

# %%
I = bempp.api.operators.boundary.sparse.identity(piecewise_const_space, piecewise_const_space, piecewise_const_space)
Mk = bempp.api.operators.boundary.helmholtz.double_layer(piecewise_const_space, piecewise_const_space, piecewise_const_space, k)
Lk = bempp.api.operators.boundary.helmholtz.single_layer(piecewise_const_space, piecewise_const_space, piecewise_const_space, k)

# %% [markdown]
# ## Calculate RHS and LHS

# %%
lhs = Mk-0.5*I
rhs = Lk*vn

# %% [markdown]
# ## Solve

# %%
phi, info = bempp.api.linalg.gmres(lhs, rhs, tol=1e-6)

# %% [markdown]
# ## Calculate sound pressure level at 1000 m

# %%
# Fast Multipole Method
point = np.vstack([0, 0, 1000])
Mk_pot = bempp.api.operators.potential.helmholtz.double_layer(piecewise_const_space, point, k, assembler="fmm")
Lk_pot = bempp.api.operators.potential.helmholtz.single_layer(piecewise_const_space, point, k, assembler="fmm")
phi_field = Mk_pot.evaluate(phi)-Lk_pot.evaluate(vn)
p_field = 1j*rho*omega*phi_field
SPL_FMM = 20*np.log10(np.abs(p_field)/1e-6)

# %%
# OpenCL
point = np.vstack([0, 0, 1000])
Mk_pot = bempp.api.operators.potential.helmholtz.double_layer(piecewise_const_space, point, k, assembler="dense", device_interface="opencl")
Lk_pot = bempp.api.operators.potential.helmholtz.single_layer(piecewise_const_space, point, k, assembler="dense", device_interface="opencl")
phi_field = Mk_pot.evaluate(phi)-Lk_pot.evaluate(vn)
p_field = 1j*rho*omega*phi_field
SPL_OpenCL = 20*np.log10(np.abs(p_field)/1e-6)

# %%
# Numba
point = np.vstack([0, 0, 1000])
Mk_pot = bempp.api.operators.potential.helmholtz.double_layer(piecewise_const_space, point, k, assembler="dense", device_interface="numba")
Lk_pot = bempp.api.operators.potential.helmholtz.single_layer(piecewise_const_space, point, k, assembler="dense", device_interface="numba")
phi_field = Mk_pot.evaluate(phi)-Lk_pot.evaluate(vn)
p_field = 1j*rho*omega*phi_field
SPL_Numba = 20*np.log10(np.abs(p_field)/1e-6)

# %%
R = np.linalg.norm(point)
p_field = rho*c*Wdot*k*a**2*(k*a-1j)/(((k*a)**2+1)*R)*np.exp(1j*k*(R-a))
SPL_analytical = 20*np.log10(np.abs(p_field)/1e-6)

# %%
row_labels = ['Analytical', 'FMM', 'OpenCL', 'Numba']
SPL = np.array([SPL_analytical, SPL_FMM, SPL_OpenCL, SPL_Numba], dtype=float)
df = pd.DataFrame(SPL, index=row_labels, columns=['SPL at 1000 m [dB re 1 uPa]'])
print(df)

[tbetcke]

Hi. Thank you for the report. There are two different issues here:

1.) You are defining the point p as int64 type (since 1000) is an Integer. If you specify p = np.vstack([0, 0, 1000.0]) OpenCL and Numba give the same results. This is an edge case we should check for. 2.) The FMM gives the wrong results because the expansion order is not sufficient. The Exafmm library that we interface is highly performant for low-frequency problems, but not suitable at high-frequencies. Because of the large distance of your evaluation there are many evaluations between the domain and the eval point, which the FMM does not capture properly. You can see this by setting bempp.api.GLOBAL_PARAMETERS.fmm.dense_evaluation = True. In this mode it uses a direct evaluation of all point to point interactions instead of just approximations for the far-field as the FMM usually does. We have built this in for debugging purposes to check FMM results (it is too slow for larger problems). But using this parameter the result becomes immediately correct.

Thank you for reporting this. I will create a separate bug report for the int64 issue and then close this bug report.

[tbetcke]

Have created bug report #58. Closing this one.

[antbonomo]

@tbetcke Thank you very much for your prompt response. Keep up the great work!