From 00e0090cac26bb9fff642d105ee33dbb745259d3 Mon Sep 17 00:00:00 2001
From: Isuru Fernando <idf2@illinois.edu>
Date: Thu, 23 Jun 2022 14:28:19 -0500
Subject: [PATCH] 4r/order for biharmonic (#122)
* 4r/order for biharmonic
* Add a test for magnitudes
---
sumpy/fmm.py | 10 ++++--
test/test_fmm.py | 84 +++++++++++++++++++++++++++++++++++++++++++++++-
2 files changed, 91 insertions(+), 3 deletions(-)
diff --git a/sumpy/fmm.py b/sumpy/fmm.py
index 519b1200..42564444 100644
--- a/sumpy/fmm.py
+++ b/sumpy/fmm.py
@@ -323,8 +323,14 @@ class SumpyExpansionWrangler(ExpansionWranglerInterface):
# YALE UNIV NEW HAVEN CT DEPT OF COMPUTER SCIENCE, 1988.
# rscale that we use in sumpy is the inverse of the scaling used in the
# paper and therefore we should use r / order. However empirically
- # we have observed that 2r / order is better for numerical stability.
- return r * 2 / order
+ # we have observed that 2r / order is better for numerical stability
+ # for Laplace and 4r / order for biharmonic kernel.
+ knl = self.tree_indep.get_base_kernel()
+ from sumpy.kernel import BiharmonicKernel
+ if isinstance(knl, BiharmonicKernel):
+ return r * 4 / order
+ else:
+ return r * 2 / order
# {{{ data vector utilities
diff --git a/test/test_fmm.py b/test/test_fmm.py
index cbf66e6d..4d6d0d48 100644
--- a/test/test_fmm.py
+++ b/test/test_fmm.py
@@ -27,7 +27,8 @@ import numpy.linalg as la
import pyopencl as cl
from pyopencl.tools import ( # noqa
pytest_generate_tests_for_pyopencl as pytest_generate_tests)
-from sumpy.kernel import LaplaceKernel, HelmholtzKernel, YukawaKernel
+from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, YukawaKernel,
+ BiharmonicKernel)
from sumpy.expansion.multipole import (
VolumeTaylorMultipoleExpansion,
H2DMultipoleExpansion, Y2DMultipoleExpansion,
@@ -234,6 +235,87 @@ def test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class,
pconv_verifier()
+@pytest.mark.parametrize("knl", [LaplaceKernel(2), BiharmonicKernel(2)])
+def test_coeff_magnitude_rscale(ctx_factory, knl):
+ """Checks that the rscale used keeps the coefficient magnitude
+ difference small
+ """
+ local_expn_class = LinearPDEConformingVolumeTaylorLocalExpansion
+ mpole_expn_class = LinearPDEConformingVolumeTaylorMultipoleExpansion
+
+ ctx = ctx_factory()
+ queue = cl.CommandQueue(ctx)
+
+ nsources = 1000
+ ntargets = 300
+ dtype = np.float64
+
+ from boxtree.tools import (
+ make_normal_particle_array as p_normal)
+
+ sources = p_normal(queue, nsources, knl.dim, dtype, seed=15)
+ offset = np.zeros(knl.dim)
+ offset[0] = 0.1
+
+ targets = (
+ p_normal(queue, ntargets, knl.dim, dtype, seed=18) + offset)
+
+ from boxtree import TreeBuilder
+ tb = TreeBuilder(ctx)
+
+ tree, _ = tb(queue, sources, targets=targets,
+ max_particles_in_box=30, debug=True)
+
+ from boxtree.traversal import FMMTraversalBuilder
+ tbuild = FMMTraversalBuilder(ctx)
+ trav, _ = tbuild(queue, tree, debug=True)
+
+ from pyopencl.clrandom import PhiloxGenerator
+ rng = PhiloxGenerator(ctx, seed=44)
+ weights = rng.uniform(queue, nsources, dtype=np.float64)
+
+ extra_kwargs = {}
+ dtype = np.float64
+ order = 10
+ if isinstance(knl, HelmholtzKernel):
+ extra_kwargs["k"] = 0.05
+ dtype = np.complex128
+
+ elif isinstance(knl, YukawaKernel):
+ extra_kwargs["lam"] = 2
+ dtype = np.complex128
+
+ from functools import partial
+ target_kernels = [knl]
+
+ tree_indep = SumpyTreeIndependentDataForWrangler(
+ ctx,
+ partial(mpole_expn_class, knl),
+ partial(local_expn_class, knl),
+ target_kernels)
+
+ def fmm_level_to_order(kernel, kernel_args, tree, lev):
+ return order
+
+ wrangler = SumpyExpansionWrangler(tree_indep, trav, dtype,
+ fmm_level_to_order=fmm_level_to_order,
+ kernel_extra_kwargs=extra_kwargs)
+
+ weights = wrangler.reorder_sources(weights)
+ (weights,) = wrangler.distribute_source_weights((weights,), None)
+
+ local_result, _ = wrangler.form_locals(
+ trav.level_start_target_or_target_parent_box_nrs,
+ trav.target_or_target_parent_boxes,
+ trav.from_sep_bigger_starts,
+ trav.from_sep_bigger_lists,
+ (weights,))
+
+ result = np.abs(wrangler.local_expansions_view(local_result, 5)[1][0])
+
+ assert np.max(result) / np.min(result) < 10**6
+
+
def test_unified_single_and_double(ctx_factory):
"""
Test that running one FMM for single layer + double layer gives the
--
GitLab