From 971af123ef03f2c22232a7ab8922def90a12ad19 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Wed, 10 Aug 2016 16:23:13 -0500 Subject: [PATCH 1/2] Add level-to-order lookup from pyfmmlib. --- setup.py | 1 + sumpy/expansion/level_to_order.py | 93 +++++++++++++++++++++++++++++++ test/test_fmm.py | 32 +++++++++++ 3 files changed, 126 insertions(+) create mode 100644 sumpy/expansion/level_to_order.py diff --git a/setup.py b/setup.py index 19924d68..c1e5f698 100644 --- a/setup.py +++ b/setup.py @@ -44,6 +44,7 @@ setup(name="sumpy", "pytools>=2013.5.6", "boxtree>=2013.1", "pytest>=2.3", + "pyfmmlib>=2016.1", "six", # If this causes issues, see: diff --git a/sumpy/expansion/level_to_order.py b/sumpy/expansion/level_to_order.py new file mode 100644 index 00000000..249acee2 --- /dev/null +++ b/sumpy/expansion/level_to_order.py @@ -0,0 +1,93 @@ +from __future__ import division + +__copyright__ = "Copyright (C) 2016 Matt Wala" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +__doc__ = """ +.. autofunction:: h2d_level_to_order_lookup +.. autofunction:: l2d_level_to_order_lookup +""" + +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` + """ + + 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 ValueError( + "h2dterms returned error code {ier}".format(ier=ier)) + orders[level] = nterms + + return orders + + +def l2d_level_to_order_lookup(tree, epsilon): + """ + Compute a mapping from level number to expansion order, + Laplace 2D case. + + This wraps the function *l2dterms* from :mod:`pyfmmlib`. + + :arg tree: An instance of :class:`boxtree.Tree`. + :arg epsilon: Precision + + :return: A :class:`numpy.array` of length `tree.nlevels` + """ + + if tree.dimensions != 2: + raise ValueError("tree must be 2d") + + from pyfmmlib import l2dterms + nterms, ier = l2dterms(epsilon) + if ier != 0: + raise ValueError( + "l2dterms returned error code {ier}".format(ier=ier)) + + orders = np.empty(tree.nlevels, dtype=int) + orders.fill(nterms) + + return orders + + +# vim: fdm=marker diff --git a/test/test_fmm.py b/test/test_fmm.py index ce26622a..fa99c56e 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -36,6 +36,10 @@ from sumpy.expansion.multipole import ( from sumpy.expansion.local import ( VolumeTaylorLocalExpansion, H2DLocalExpansion) +from sumpy.expansion.level_to_order import ( + h2d_level_to_order_lookup, + l2d_level_to_order_lookup) + import pytest import logging @@ -50,6 +54,34 @@ else: faulthandler.enable() +@pytest.mark.parametrize("lookup_func, extra_args", [ + (h2d_level_to_order_lookup, (100,)), + (l2d_level_to_order_lookup, ()), + ]) +def test_level_to_order_lookup(ctx_getter, lookup_func, extra_args): + ctx = ctx_getter() + queue = cl.CommandQueue(ctx) + + nsources = 500 + dtype = np.float64 + epsilon = 1e-9 + + from boxtree.tools import ( + make_normal_particle_array as p_normal) + + sources = p_normal(queue, nsources, 2, dtype, seed=15) + from boxtree import TreeBuilder + tb = TreeBuilder(ctx) + + tree, _ = tb(queue, sources, max_particles_in_box=30, debug=True) + + levels = lookup_func(tree, *(extra_args + (epsilon,))) + + assert len(levels) == tree.nlevels + # Should be monotonically nonincreasing. + assert all(levels[i] >= levels[i+1] for i in range(tree.nlevels - 1)) + + @pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class", [ (LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), (LaplaceKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), -- GitLab From 13fbf2a928589046c152e1994a375b98206160ae Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Wed, 10 Aug 2016 16:39:19 -0500 Subject: [PATCH 2/2] ValueError => RuntimeError. --- sumpy/expansion/level_to_order.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sumpy/expansion/level_to_order.py b/sumpy/expansion/level_to_order.py index 249acee2..b0e4ef26 100644 --- a/sumpy/expansion/level_to_order.py +++ b/sumpy/expansion/level_to_order.py @@ -55,7 +55,7 @@ def h2d_level_to_order_lookup(tree, helmholtz_k, epsilon): for level in range(tree.nlevels): nterms, ier = h2dterms(bbox_area / 2 ** level, helmholtz_k, epsilon) if ier != 0: - raise ValueError( + raise RuntimeError( "h2dterms returned error code {ier}".format(ier=ier)) orders[level] = nterms @@ -81,7 +81,7 @@ def l2d_level_to_order_lookup(tree, epsilon): from pyfmmlib import l2dterms nterms, ier = l2dterms(epsilon) if ier != 0: - raise ValueError( + raise RuntimeError( "l2dterms returned error code {ier}".format(ier=ier)) orders = np.empty(tree.nlevels, dtype=int) -- GitLab