inducer / sumpy

Symbolic code generators for multipole and local expansions and translations
32 stars 13 forks source link

test_translations fails when run with 3D kernels #107

Closed isuruf closed 2 years ago

isuruf commented 2 years ago

Currently only 2D kernels are tested and when 3D kernels are tested the convergence verifier fails.

inducer commented 2 years ago

Sample convergence data for Laplace? What stage fails?

isuruf commented 2 years ago

Laplace fails at p2m2m2p while Helmholtz fails at p2m2p.

Laplace 3D - VolumeTaylorExpansion

p2m2p
------------------------------
p | error                 
--+------------------------
2 | 2.767760520941984e-05 
3 | 5.776971943322948e-06 
4 | 1.7429094830716893e-07
------------------------------
------------------------------
p2m2m2p
------------------------------
p | error                 
--+------------------------
2 | 1.2065883172728181e-05
3 | 9.31726501127403e-06  
4 | 1.2322093264050907e-06
------------------------------

Helmholtz 3D - VolumeTaylorExpansion

p2m2p
------------------------------
p | error                
--+-----------------------
2 | 0.0002609529429951883
3 | 0.0002779895576949443
4 | 0.0002829015713450385
------------------------------
inducer commented 2 years ago

I can reproduce that. Not super worried about the Laplace result, as 3->4 shows a healthy decrease. (and 2->3 causes the issue)

Order 5 keeps getting better, so that's IMO just the test that needs fixing:

------------------------------
p | error                 
--+------------------------
2 | 1.2065883172727886e-05
3 | 9.317265011274165e-06 
4 | 1.2322093264059911e-06
5 | 5.730334486283711e-08 
------------------------------

Helmholtz seems well and truly stalled though. And that's just expand and evaluate...

inducer commented 2 years ago

I get a sense that the 3D kernel derivatives for Helmholtz must be wrong.

Consider the following code (using the "expansion toys"): https://gist.github.com/inducer/47424a1d9887256505da9e7d9b906e3a

For Laplace 3D, the error plot of the 5-th order multipole looks the way I expect, with O(2*order) lines of zero error radiating out:

grafik

For Helmholtz (also order 5), this is not the case:

grafik

And since the multipole coefficients for both are actually the same, this clearly points to the derivatives to me. It should be possible to check those using, e.g. the CalculusPatch. Could you take a look?

isuruf commented 2 years ago

Yes, you are right. I changed the derivative taker to the sympy based one and the error changes. (Still have the same error as the Laplace one, but that seems like a test issue instead of a sumpy code issue.)

diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index 635f769..0c00b11 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -562,8 +562,8 @@ class HelmholtzKernel(ExpressionKernel):
         """Return a :class:`sumpy.tools.ExprDerivativeTaker` instance that supports
         taking derivatives of the base kernel with respect to dvec.
         """
-        from sumpy.tools import HelmholtzDerivativeTaker
-        return HelmholtzDerivativeTaker(self.get_expression(dvec), dvec, rscale, sac)
+        from sumpy.tools import HelmholtzDerivativeTaker, ExprDerivativeTaker
+        return ExprDerivativeTaker(self.get_expression(dvec), dvec, rscale, sac)

     def get_pde_as_diff_op(self):
         from sumpy.expansion.diff_op import make_identity_diff_op, laplacian