Helmholtz wrappers's self-interaction handling when source and target points are the same

[yongzhongli]

Hi FMM3D developers,

I am testing hfmm3d in Julia and noticed a discrepancy when setting the target points equal to the source points (a situation that naturally arises in integral equation methods).

According to the documentation:

“When x = xj, the term corresponding to xj is dropped from the sum.”

From my understanding, this means that the Green’s function contribution for coincident source–target pairs is ignored. Is this equivalent to treating the kernel value G(0) as zero when r=0?

And I found when I use hfmm3d with targets identical to sources, the resulting potentials differ from those computed with h3ddir. Should I expect hfmm3d and h3ddir results to match when targets coincide with sources, or is the discrepancy expected?

For integral equation contexts where self-terms must be treated carefully (e.g., singular integrals), is the recommendation to always avoid exact coincident source–target evaluations with hfmm3d?

[lu1and10]

Hi @yongzhongli

In the code, we do treat G(0) as zero, as in the direct evaluation kernels, G(r) is set to zero for all r < thresh,

FMM3D/src/Helmholtz/helmkernels.f

Line 104 in 70455e0

if(d.lt.thresh) goto 1000

The thresh used in hfmm3d is,

FMM3D/src/Helmholtz/hfmm3d.f

Line 778 in 70455e0

thresh = 2.0d0**(-51)*boxsize(0)

While in the julia wrapper, the h3ddir default thresh is set to

FMM3D/julia/src/FMM3D.jl

Line 460 in 70455e0

thresh::Float64=1e-16)

unless you pass a thresh.

hfmm3d should differ from h3ddir by requested tolerance.
It would be great that you could provide us a minimal reproducible code example as your test reproducing the discrepancy you find and we could debug from it.

[yongzhongli]

Hi @lu1and10,

Thanks for your explanation, the code snippet looks good to me, but I don't know why in my case it is not working properly. I can use other reference to confirm the h3ddir gives correct results, but hfmm3d doesn't. If they both set G = 0 when r is small, I wonder why the results should differ.

Here is a minimal reproducible code that can reproduce the discrepancy I found

# Minimal reproducible example for FMM3D issue:
# Case: exact coincidence (targets == sources)
# Kernel convention: G = exp(-i k r) / (4π r)
# -> call FMM3D with zk = -k, then scale by 1/(4π).

using LinearAlgebra, Random
using FMM3D


function unit_sphere_points(M; rng=Random.default_rng())
    X = randn(rng, 3, M)
    @inbounds for j in 1:M
        X[:, j] /= norm(X[:, j])
    end
    return X
end

function eval_both(sources, targets; k=5.0, eps=1e-6, charges=ones(ComplexF64, size(sources,2)))
    zk = complex(-k, 0.0)
    @time fmm = FMM3D.hfmm3d(eps, zk, sources; targets=targets, charges=charges, pgt=1)
    @time dir = FMM3D.h3ddir(zk, sources, targets; charges=charges, pgt=1)
    u_fmm = fmm.pottarg
    u_dir = dir.pottarg
    return u_fmm, u_dir
end

function main(; M=10000, k=2*3.1415, eps=1e-6, seed=1234)
    rng = MersenneTwister(seed)
    sources = unit_sphere_points(M; rng)
    charges = randn(rng, M) .+ 1im * randn(rng, M)
    targets = sources
    u_fmm, u_dir = eval_both(sources, targets; k=k, eps=eps, charges=charges)
    rel_err = norm(u_fmm - u_dir) / norm(u_dir)
    println("Relative error = ", rel_err)
    for j in 1:5
        println("j=$j  u_fmm=$(u_fmm[j])   u_dir=$(u_dir[j])")
    end
end

main()

I found that only if setting eps = 1e-16 in hfmm3d makes its results match h3ddir, but then both take about the same time.

[lu1and10]

Hi @yongzhongli

I think the current Helmholtz FMM implementation requires zk in the first quadrant. You may try zk in the first quadrant(real positive zk in your case) for different tolerances.

The reason you get FMM matches the direct calculation for tolerance 1e-16 is that the code switches to direct calculation internally for that tolerance.