From c04082c578a74678f9e38bda3817b226e9f52f1b Mon Sep 17 00:00:00 2001
From: Matt Wala <wala1@illinois.edu>
Date: Fri, 19 Jul 2019 16:41:27 -0500
Subject: [PATCH] Delete broken/obsolete h2dterms and l2dterms wrappers, and
replace with FMMLib-based expansion order finder
Closes #39
- [ ] Point requirements.txt back to pyfmmlib master after pyfmmlib!15 is merged
---
requirements.txt | 2 +-
sumpy/expansion/level_to_order.py | 90 ++++++++++++++-----------------
test/test_misc.py | 20 ++++++-
3 files changed, 59 insertions(+), 53 deletions(-)
diff --git a/requirements.txt b/requirements.txt
index 8e48bc1c..ccf135bd 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -5,4 +5,4 @@ git+https://github.com/inducer/islpy
git+https://github.com/inducer/pyopencl
git+https://gitlab.tiker.net/inducer/boxtree
git+https://github.com/inducer/loopy
-git+https://github.com/inducer/pyfmmlib
+git+https://gitlab.tiker.net/inducer/pyfmmlib@l3dterms-and-h3dterms
diff --git a/sumpy/expansion/level_to_order.py b/sumpy/expansion/level_to_order.py
index 8f012e55..5e29b861 100644
--- a/sumpy/expansion/level_to_order.py
+++ b/sumpy/expansion/level_to_order.py
@@ -23,73 +23,61 @@ THE SOFTWARE.
"""
__doc__ = """
-.. autofunction:: h2d_level_to_order_lookup
-.. autofunction:: l2d_level_to_order_lookup
-
+.. autoclass:: FMMLibExpansionOrderFinder
.. autoclass:: SimpleExpansionOrderFinder
"""
import numpy as np
-def h2d_level_to_order_lookup(tree, helmholtz_k, epsilon):
- """
- Compute a mapping from level number to expansion order,
- Helmholtz 2D case.
-
- This wraps the function *h2dterms* from :mod:`pyfmmlib`.
-
- :arg tree: An instance of :class:`boxtree.Tree`.
- :arg helmholtz_k: Helmholtz parameter
- :arg epsilon: Precision
-
- :return: A :class:`numpy.array` of length `tree.nlevels`
+class FMMLibExpansionOrderFinder(object):
+ r"""Return expansion orders that meet the tolerance for a given level
+ using routines wrapped from :mod:`pyfmmlib`.
"""
- if tree.dimensions != 2:
- raise ValueError("tree must be 2d")
-
- orders = np.empty(tree.nlevels, dtype=int)
- bbox_area = np.max(
- tree.bounding_box[1] - tree.bounding_box[0]) ** 2
-
- from pyfmmlib import h2dterms
- for level in range(tree.nlevels):
- nterms, ier = h2dterms(bbox_area / 2 ** level, helmholtz_k, epsilon)
- if ier != 0:
- raise RuntimeError(
- "h2dterms returned error code {ier}".format(ier=ier))
- orders[level] = nterms
-
- return orders
+ def __init__(self, tol, extra_order=0):
+ """
+ :arg tol: error tolerance
+ :arg extra_order: order increase to accommodate, say, the taking of
+ derivatives of the FMM expansions.
+ """
+ self.tol = tol
+ self.extra_order = extra_order
+ def __call__(self, kernel, kernel_args, tree, level):
+ import pyfmmlib
-def l2d_level_to_order_lookup(tree, epsilon):
- """
- Compute a mapping from level number to expansion order,
- Laplace 2D case.
+ from sumpy.kernel import LaplaceKernel, HelmholtzKernel
- This wraps the function *l2dterms* from :mod:`pyfmmlib`.
+ assert isinstance(kernel, (LaplaceKernel, HelmholtzKernel))
+ assert tree.dimensions in (2, 3)
- :arg tree: An instance of :class:`boxtree.Tree`.
- :arg epsilon: Precision
+ if isinstance(kernel, LaplaceKernel):
+ if tree.dimensions == 2:
+ nterms, ier = pyfmmlib.l2dterms(self.tol)
+ if ier:
+ raise RuntimeError("l2dterms returned error code '%d'" % ier)
- :return: A :class:`numpy.array` of length `tree.nlevels`
- """
+ elif tree.dimensions == 3:
+ nterms, ier = pyfmmlib.l3dterms(self.tol)
+ if ier:
+ raise RuntimeError("l3dterms returned error code '%d'" % ier)
- if tree.dimensions != 2:
- raise ValueError("tree must be 2d")
+ elif isinstance(kernel, HelmholtzKernel):
+ helmholtz_k = dict(kernel_args)[kernel.helmholtz_k_name]
+ size = tree.root_extent / 2 ** level
- from pyfmmlib import l2dterms
- nterms, ier = l2dterms(epsilon)
- if ier != 0:
- raise RuntimeError(
- "l2dterms returned error code {ier}".format(ier=ier))
+ if tree.dimensions == 2:
+ nterms, ier = pyfmmlib.h2dterms(size, helmholtz_k, self.tol)
+ if ier:
+ raise RuntimeError("h2dterms returned error code '%d'" % ier)
- orders = np.empty(tree.nlevels, dtype=int)
- orders.fill(nterms)
+ elif tree.dimensions == 3:
+ nterms, ier = pyfmmlib.h3dterms(size, helmholtz_k, self.tol)
+ if ier:
+ raise RuntimeError("h3dterms returned error code '%d'" % ier)
- return orders
+ return nterms + self.extra_order
class SimpleExpansionOrderFinder(object):
@@ -114,7 +102,7 @@ class SimpleExpansionOrderFinder(object):
extra_order=1):
"""
:arg extra_order: order increase to accommodate, say, the taking of
- oderivatives f the FMM expansions.
+ derivatives of the FMM expansions.
"""
self.tol = tol
diff --git a/test/test_misc.py b/test/test_misc.py
index 33e8535a..8866c289 100644
--- a/test/test_misc.py
+++ b/test/test_misc.py
@@ -124,7 +124,9 @@ class FakeTree:
self.stick_out_factor = stick_out_factor
-@pytest.mark.parametrize("knl", [LaplaceKernel(2), HelmholtzKernel(2)])
+@pytest.mark.parametrize("knl", [
+ LaplaceKernel(2), HelmholtzKernel(2),
+ LaplaceKernel(3), HelmholtzKernel(3)])
def test_order_finder(knl):
from sumpy.expansion.level_to_order import SimpleExpansionOrderFinder
@@ -137,6 +139,22 @@ def test_order_finder(knl):
print(orders)
+@pytest.mark.parametrize("knl", [
+ LaplaceKernel(2), HelmholtzKernel(2),
+ LaplaceKernel(3), HelmholtzKernel(3)])
+def test_fmmlib_order_finder(knl):
+ pytest.importorskip("pyfmmlib")
+ from sumpy.expansion.level_to_order import FMMLibExpansionOrderFinder
+
+ ofind = FMMLibExpansionOrderFinder(1e-5)
+
+ tree = FakeTree(knl.dim, 200, 0.5)
+ orders = [
+ ofind(knl, frozenset([("k", 5)]), tree, level)
+ for level in range(30)]
+ print(orders)
+
+
# {{{ expansion toys p2e2e2p test cases
def approx_convergence_factor(orders, errors):
--
GitLab