diff --git a/sumpy/codegen.py b/sumpy/codegen.py
index 948512dfe884c35e0f421c4a4314796e38b0ab7b..bef7fc6bf93a171a96626f9db88c312c23167bf2 100644
--- a/sumpy/codegen.py
+++ b/sumpy/codegen.py
@@ -581,6 +581,28 @@ class BigIntegerKiller(CSECachingMapperMixin, IdentityMapper):
# }}}
+# {{{ convert 123000000j to 123000000 * 1j
+
+class ComplexRewriter(CSECachingMapperMixin, IdentityMapper):
+
+ def __init__(self, float_type=np.float32):
+ IdentityMapper.__init__(self)
+ self.float_type = float_type
+
+ def map_constant(self, expr):
+ """Convert complex values not within complex64 to a product for loopy
+ """
+ if not isinstance(expr, complex):
+ return IdentityMapper.map_constant(self, expr)
+
+ if complex(self.float_type(expr.imag)) == expr.imag:
+ return IdentityMapper.map_constant(self, expr)
+
+ return expr.real + prim.Product((expr.imag, 1j))
+
+ map_common_subexpression_uncached = IdentityMapper.map_common_subexpression
+
+
# {{{ vector component rewriter
INDEXED_VAR_RE = re.compile("^([a-zA-Z_]+)([0-9]+)$")
@@ -684,6 +706,7 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[],
ssg = SumSignGrouper()
fck = FractionKiller()
bik = BigIntegerKiller()
+ cmr = ComplexRewriter()
def convert_expr(name, expr):
logger.debug("generate expression for: %s" % name)
@@ -694,6 +717,7 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[],
expr = fck(expr)
expr = ssg(expr)
expr = bik(expr)
+ expr = cmr(expr)
#expr = cse_tag(expr)
for m in pymbolic_expr_maps:
expr = m(expr)
@@ -705,7 +729,7 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[],
result = [
lp.Assignment(id=None,
assignee=name, expression=convert_expr(name, expr),
- temp_var_type=lp.auto)
+ temp_var_type=lp.Optional(None))
for name, expr in assignments]
logger.info("loopy instruction generation: done")
diff --git a/sumpy/e2p.py b/sumpy/e2p.py
index d0d05cbf8c45e8d9f1c4552f784950bdfcab4bfb..7b0072ad55708ba205ceed3448914ad0d8eda281 100644
--- a/sumpy/e2p.py
+++ b/sumpy/e2p.py
@@ -119,7 +119,7 @@ class E2PBase(KernelCacheWrapper):
assignee="kernel_scaling",
expression=sympy_conv(
self.expansion.kernel.get_global_scaling_const()),
- temp_var_type=lp.auto)]
+ temp_var_type=lp.Optional(None))]
def get_cache_key(self):
return (type(self).__name__, self.expansion, tuple(self.kernels))
diff --git a/sumpy/expansion/level_to_order.py b/sumpy/expansion/level_to_order.py
index 8f012e55827f7676e2429f76502f5146c085387c..5e29b86138ec705ef57747929dead3403160b17a 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/sumpy/kernel.py b/sumpy/kernel.py
index 99222f684340fe1f36d888355a24cbffe29573c0..aebaa9fe76da1088f33dc38b3b950591732d3287 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -431,13 +431,14 @@ class BiharmonicKernel(ExpressionKernel):
init_arg_names = ("dim",)
def __init__(self, dim=None):
+ # See https://arxiv.org/abs/1202.1811
r = pymbolic_real_norm_2(make_sym_vector("d", dim))
if dim == 2:
- expr = r**2 * var("log")(r)
- scaling = 1/(8*var("pi"))
+ expr = r**2 * (var("log")(r) - 1)
+ scaling = -1/(8*var("pi"))
elif dim == 3:
expr = r
- scaling = 1 # FIXME: Unknown
+ scaling = 1/(8*var("pi"))
else:
raise RuntimeError("unsupported dimensionality")
diff --git a/sumpy/p2p.py b/sumpy/p2p.py
index 4c4eb945bb90d2d71bcb5ed3327d744e4ad26391..e3b457dd52e885bf66462e475b3b58edb78bd098 100644
--- a/sumpy/p2p.py
+++ b/sumpy/p2p.py
@@ -125,7 +125,7 @@ class P2PBase(KernelComputation, KernelCacheWrapper):
return [lp.Assignment(id=None,
assignee="pair_result_%d" % i, expression=expr,
- temp_var_type=lp.auto)
+ temp_var_type=lp.Optional(None))
for i, expr in enumerate(exprs)]
def get_default_src_tgt_arguments(self):
diff --git a/sumpy/qbx.py b/sumpy/qbx.py
index 0d4359b6566e146d66f717d25c92569803893d0d..9708764c01d0448d7c7e8314efb461989c838acd 100644
--- a/sumpy/qbx.py
+++ b/sumpy/qbx.py
@@ -143,7 +143,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper):
return [lp.Assignment(id=None,
assignee="pair_result_%d" % i, expression=expr,
- temp_var_type=lp.auto)
+ temp_var_type=lp.Optional(None))
for i, expr in enumerate(exprs)]
def get_default_src_tgt_arguments(self):
diff --git a/sumpy/tools.py b/sumpy/tools.py
index c826ef2fe8ea4d3f96e3b8b816dd3500f5efaae8..540208f342e429326b7571168b754f02c0041d79 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -299,7 +299,7 @@ class KernelComputation(object):
lp.Assignment(id=None,
assignee="knl_%d_scaling" % i,
expression=sympy_conv(kernel.get_global_scaling_const()),
- temp_var_type=dtype)
+ temp_var_type=lp.Optional(dtype))
for i, (kernel, dtype) in enumerate(
zip(self.kernels, self.value_dtypes))]
@@ -559,10 +559,13 @@ class MatrixBlockIndexRanges(object):
# Author: Raymond Hettinger
# License: MIT
-import collections
+try:
+ from collections.abc import MutableSet
+except ImportError:
+ from collections import MutableSet
-class OrderedSet(collections.MutableSet):
+class OrderedSet(MutableSet):
def __init__(self, iterable=None):
self.end = end = []
diff --git a/test/test_fmm.py b/test/test_fmm.py
index 71e3f044432b4bd612337bdd55d56e421540bcc4..2bc6d69d6bcafedf338cc78aea0591d6d5ac3845 100644
--- a/test/test_fmm.py
+++ b/test/test_fmm.py
@@ -42,10 +42,6 @@ from sumpy.expansion.local import (
LaplaceConformingVolumeTaylorLocalExpansion,
HelmholtzConformingVolumeTaylorLocalExpansion)
-from sumpy.expansion.level_to_order import (
- h2d_level_to_order_lookup,
- l2d_level_to_order_lookup)
-
import pytest
import logging
@@ -60,35 +56,6 @@ 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):
- pytest.importorskip("pyfmmlib")
- 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(2), LaplaceConformingVolumeTaylorLocalExpansion,
diff --git a/test/test_misc.py b/test/test_misc.py
index 33e8535a8cbce064263958962421d9b43e77ee6c..a25b2e800ab6f9942395821cf0cb4345e80b6ccf 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
@@ -136,6 +138,28 @@ def test_order_finder(knl):
for level in range(30)]
print(orders)
+ # Order should not increase with level
+ assert (np.diff(orders) <= 0).all()
+
+
+@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)
+
+ # Order should not increase with level
+ assert (np.diff(orders) <= 0).all()
+
# {{{ expansion toys p2e2e2p test cases