Fix incorrect scaling of Yukawa kernel

[inducer]

Noticed by Shawn Lin.

[alexfikl]

There's a test in pytential for the Yukawa-Beltrami operator that would fail with this change (right?)
https://github.com/inducer/pytential/blob/d0321e0681390a4486def0613a59ed7c5c8b3741/test/test_beltrami.py#L131-L140

I haven't looked at this in quite a while, but I'm not seeing why there's an extra -, since it seems to match the Laplace version. Maybe worth modifying the comments in there to say why the - goes away?

EDIT: Well, apparently the tests don't fail, but an update to the comments there would still be nice 😁

[inducer]

@alexfikl Thanks for following up. I agree that this change needs to be backed up better, and I think what I'd actually like is a test for the scaling. What we'd really want to test is $\mathcal L G = \delta$, but I don't know how to write a sufficiently general-purpose test for that. Testing a jump relation for the double layer (possibly via QBX) would be another option, though that'd be much easier to do in pytential than sumpy. (Or we could allow the sumpy tests to depend on pytential?) There's precedent for such a test, so this could be a relatively simple copy-paste job.

What do you guys think?

@ShawnL00, your input/help would be very welcome here.

[alexfikl]

@alexfikl Thanks for following up. I agree that this change needs to be backed up better, and I think what I'd actually like is a test for the scaling. What we'd really want to test is L G = δ , but I don't know how to write a sufficiently general-purpose test for that.

This might be silly, but I'm still confused why that Beltrami test in pytential is passing. From what I recall, it tries to solve the eigenvalue problem $\Delta_\Sigma \sigma - k^2 \sigma = \lambda \sigma$. Why does flipping the sign here still work, since neither the exact eigenvalue or eigenfunction in the test changes? I would expect at least one of them to flip sign..

Testing a jump relation for the double layer (possibly via QBX) would be another option, though that'd be much easier to do in pytential than sumpy. (Or we could allow the sumpy tests to depend on pytential?) There's precedent for such a test, so this could be a relatively simple copy-paste job.

I'm ok with putting more jump tests in pytential for these things. However, from what I recall, the 3D jump tests in there are already quite slow, so it might add quite a bit to the runtime :\

[ShawnL00]

What we'd really want to test is L G = δ , but I don't know how to write a sufficiently general-purpose test for that. Testing a jump relation for the double layer (possibly via QBX) would be another option, though that'd be much easier to do in pytential than sumpy. (Or we could allow the sumpy tests to depend on pytential?)

Jump test is definitely a good choice. Perhaps we could also explicitly include the leading singularity for $G$ and its related operator $L$. Then, we know what's the expected jump relation.
To make sure we have $L G = \delta$ and not $L G = -\delta$, I do not think analytic derivation / reference is a bad option.

[inducer]

I do not think analytic derivation / reference is a bad option.

I'm more after something that can be verified with code though.

explicitly include the leading singularity for G and its related operator L

• kernels know their PDE (get_pde_as_diff_op)
• How would we check the leading singularity?

[ShawnL00]

• How would we check the leading singularity?

I am not sure if this is a good method: to verify the leading singularity, we can define a new kernel by subtracting out the leading singularity, evaluate the double layer along some rays that pass through the boundaries (assume the kernel of the double layer potential is not singular), and then check the continuity conditions on the boundary.

• Or we could use QBX to compute two-sided limits.

[inducer]

I'm still confused why that Beltrami test in pytential is passing

I agree; an eigenvalue test should be affected by a change in scaling. Possibly unrelated: Do you recall the source for the eigenvalues in the test? Also: given that it seems to be an eigenvalue test with a known eigenfunction, how come it does a solve (as opposed to a matvec)?

[inducer]

I can report at least one bit of happy news: Without this PR,

pycl test_layer_pot.py 'test_3d_jump_relations(_acf, "sp")' 

fails in pytential, changed to use the Yukawa kernel. It passes with this PR. So @ShawnL00 is definitely onto something.

[ShawnL00]

If I understand correctly, suppose $\lambda$ and $\phi$ are eigenvalues and eigenfunctions for the $-\Delta_{\Gamma}$ operator. We generate the right-hand side for the Yukawa-Beltrami case by setting $f := (-\Delta_{\Gamma} + k^2)\phi = (\lambda + k^2)\phi$. Then, we solve a preconditioned linear system to obtain the density, and $\phi$ is approximated by a single layer potential of the Yukawa kernel. The sign and scaling changes in the Yukawa kernel should not impact the final approximation.