Improve support for vector valued kernels

[arieder]

Hello,

I've been experimenting with your IFGF code, trying to replace the fast multipole component of bempp-cl with it. In order to be able to compute the double layer and adjoint double layer operator, I need to be able to compute the gradient of the Helmholtz potential, i.e., something vector valued.
Is there a good way to do this, or should I resort to a separate IFGF operator for each component? (Just for reference: exafmm, the FMM library used by bempp always computes the function value and the gradient in a single pass).

thank you for your input,
Alex

[maltezfaria]

Hi Alex,

It is possible to do that without having to manually interpolate each component separately, but I am afraid things are a bit hacky at this stage. Here is how you could do it (you need to be on the lfaria/dev branch for this to work):

using IFGF, LinearAlgebra, StaticArrays
import IFGF: wavenumber

# random points on a cube
const Point3D = SVector{3,Float64}
m,n = 100_000, 100_000
X,Y = rand(Point3D,m), rand(Point3D,n)

# define a the kernel matrix
struct HelmholtzGradientMatrix <: AbstractMatrix{Adjoint{ComplexF64, SVector{3, ComplexF64}}}
    X::Vector{Point3D}
    Y::Vector{Point3D}
    k::Float64
end

# indicate that this is an ocillatory kernel with wavenumber `k`
wavenumber(A::HelmholtzGradientMatrix) = A.k

# functor interface
function (K::HelmholtzGradientMatrix)(x,y)
    k = K.k
    r = x - y
    d = sqrt(transpose(r)*r)
    return (!iszero(d))*(1/(4π)/d^2 * exp(im*k*d)*(-im*k+1/d)*r')
end

function IFGF.centered_factor(K::HelmholtzGradientMatrix,x,Y)
    k = K.k
    yc = IFGF.center(Y)
    r = x - yc
    d = sqrt(transpose(r)*r)
    return exp(im*k*d)/d
end

# abstract matrix interface
Base.size(A::HelmholtzGradientMatrix) = length(X), length(Y)
Base.getindex(A::HelmholtzGradientMatrix,i::Int,j::Int) = A(A.X[i],A.Y[j])

# create the abstract matrix
k = 8π
A = HelmholtzGradientMatrix(X,Y,k)

L = assemble_ifgf(A,X,Y;order=(5,5,5))

gb = Base.summarysize(L) / 1e9
println("Size of L: $gb gigabytes")

x     = randn(SVector{3,ComplexF64},n)
I     = rand(1:m,100)
exact = [sum(A[i,j]*x[j] for j in 1:n) for i in I] # matrix-vector product
y     = L*x
er    = norm(y[I]-exact) / norm(exact) # relative error
println("approximate relative error: $er")

[maltezfaria]

But I agree with you that it would be nice to have an easy way to "assemble" the four classical integral operators for some standard PDE's such as Helmholtz as part of the API! I am gonna leave this open until such an API is implemented 😅

(I am in middle of writing a paper on the adaptive IFGF algorithm implemented here, so hopefully I will find some free time soon to finish this package and register it as well!)

[mipals]

Hi again, coming over from #32

So, if I understand this correctly then due to the fact that the IFGF requires an actual function we must apply it on the each of the components of the gradient?

Anyhow, when I try to run the example that you posted in this thread I get the following error when I come to the assembly

julia> L = assemble_ifgf(A,X,Y;order=(5,5,5))
ERROR: kernel type Adjoint{ComplexF64, SVector{3, ComplexF64}} not recognized
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] _density_type_from_kernel_type(T::Type)
   @ IFGF ~/.julia/packages/IFGF/8R4yQ/src/utils.jl:114
 [3] assemble_ifgf(kernel::HelmholtzGradientMatrix, Xpoints::Vector{SVector{3, Float64}}, Ypoints::Vector{SVector{3, Float64}}; nmax::Int64, adm::Function, tol::Float64, order::Tuple{Int64, Int64, Int64}, Δs::Float64, threads::Bool)
   @ IFGF ~/.julia/packages/IFGF/8R4yQ/src/ifgfoperator.jl:96
 [4] top-level scope
   @ REPL[28]:1

[maltezfaria]

Hi Mikkel,

Sorry, I need to set aside some time to polish and publish this package. Are you on the lfaria/dev branch? You need the commit 13efffc for this to work I think.

[mipals]

Ahhh, I completely missed that I needed the other branch. It works on here 👍