From 8d659e84f5cd45d4502e4a7e1bf87114960f09f7 Mon Sep 17 00:00:00 2001
From: Matt Wala <wala1@illinois.edu>
Date: Fri, 6 Feb 2015 21:00:41 -0600
Subject: [PATCH] Get the FMM working with Bessel expansions on the 2D
Helmholtz kernel.
- Add a H2DMultipoleExpansion class.
- Implement L2L, M2L, M2M operators.
- Teach loopy about casting to complex numbers (to handle a corner
case where some multipole coefficients may look real valued).
- Update the tests.
---
sumpy/codegen.py | 20 +++++++++++
sumpy/expansion/local.py | 68 +++++++++++++++++++++++++-----------
sumpy/expansion/multipole.py | 61 +++++++++++++++++++++++++++++++-
sumpy/kernel.py | 5 +--
test/test_fmm.py | 25 +++++++------
test/test_kernels.py | 20 ++++++-----
6 files changed, 156 insertions(+), 43 deletions(-)
diff --git a/sumpy/codegen.py b/sumpy/codegen.py
index 29a2543e..1c91e7bf 100644
--- a/sumpy/codegen.py
+++ b/sumpy/codegen.py
@@ -65,6 +65,26 @@ class SympyToPymbolicMapper(SympyToPymbolicMapperBase):
# }}}
+# {{{ complex handling
+
+def complex_mangler(identifier, arg_dtypes):
+ """A function "mangler" to make casting to complex values
+ digestible for :mod:`loopy`.
+
+ See argument *function_manglers* to :func:`loopy.make_kernel`.
+ """
+
+ if identifier == "complex" and arg_dtypes == (np.dtype(np.float64),):
+ return (np.dtype(np.complex128), "cdouble_fromreal",
+ (np.dtype(np.float64),))
+ elif identifier == "complex" and arg_dtypes == (np.dtype(np.complex128),):
+ return (np.dtype(np.complex128), "cdouble_cast",
+ (np.dtype(np.complex128),))
+ else:
+ return None
+
+# }}}
+
# {{{ bessel handling
diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py
index 699fa3f4..bff7acee 100644
--- a/sumpy/expansion/local.py
+++ b/sumpy/expansion/local.py
@@ -123,34 +123,60 @@ class H2DLocalExpansion(LocalExpansionBase):
return range(-self.order, self.order+1)
def coefficients_from_source(self, avec, bvec):
- hankel_1 = sp.Function("hankel_1")
- center_source_angle = sp.atan2(-avec[1], -avec[0])
-
from sumpy.symbolic import sympy_real_norm_2
- u = sympy_real_norm_2(avec)
-
- e_i_csangle = sp.exp(sp.I*center_source_angle)
- return [
- self.kernel.postprocess_at_source(
- hankel_1(i, sp.Symbol("k")*u)*e_i_csangle**i,
- avec)
- for i in self.get_coefficient_identifiers()]
+ hankel_1 = sp.Function("hankel_1")
+ # The coordinates are negated since avec points from source to center.
+ source_angle_rel_center = sp.atan2(-avec[1], -avec[0])
+ avec_len = sympy_real_norm_2(avec)
+ return [self.kernel.postprocess_at_source(
+ hankel_1(l, sp.Symbol("k") * avec_len)
+ * sp.exp(sp.I * l * source_angle_rel_center), avec)
+ for l in self.get_coefficient_identifiers()]
def evaluate(self, coeffs, bvec):
+ from sumpy.symbolic import sympy_real_norm_2
bessel_j = sp.Function("bessel_j")
+ bvec_len = sympy_real_norm_2(bvec)
+ target_angle_rel_center = sp.atan2(bvec[1], bvec[0])
+ return sum(coeffs[self.get_storage_index(l)]
+ * self.kernel.postprocess_at_target(
+ bessel_j(l, sp.Symbol("k") * bvec_len)
+ * sp.exp(sp.I * l * -target_angle_rel_center), bvec)
+ for l in self.get_coefficient_identifiers())
+ def translate_from(self, src_expansion, src_coeff_exprs, dvec):
from sumpy.symbolic import sympy_real_norm_2
- v = sympy_real_norm_2(bvec)
-
- center_target_angle = sp.atan2(bvec[1], bvec[0])
- e_i_ctangle = sp.exp(-sp.I*center_target_angle)
- return sum(
- coeffs[self.get_storage_index(i)]
- * self.kernel.postprocess_at_target(
- bessel_j(i, sp.Symbol("k")*v)
- * e_i_ctangle**i, bvec)
- for i in self.get_coefficient_identifiers())
+ if isinstance(src_expansion, H2DLocalExpansion):
+ dvec_len = sympy_real_norm_2(dvec)
+ bessel_j = sp.Function("bessel_j")
+ new_center_angle_rel_old_center = sp.atan2(dvec[1], dvec[0])
+ translated_coeffs = []
+ for l in self.get_coefficient_identifiers():
+ translated_coeffs.append(
+ sum(src_coeff_exprs[src_expansion.get_storage_index(m)]
+ * bessel_j(m - l, sp.Symbol("k") * dvec_len)
+ * sp.exp(sp.I * (m - l) * -new_center_angle_rel_old_center)
+ for m in src_expansion.get_coefficient_identifiers()))
+ return translated_coeffs
+
+ from sumpy.expansion.multipole import H2DMultipoleExpansion
+ if isinstance(src_expansion, H2DMultipoleExpansion):
+ dvec_len = sympy_real_norm_2(dvec)
+ hankel_1 = sp.Function("hankel_1")
+ new_center_angle_rel_old_center = sp.atan2(dvec[1], dvec[0])
+ translated_coeffs = []
+ for l in self.get_coefficient_identifiers():
+ translated_coeffs.append(
+ sum((-1) ** l * hankel_1(m + l, sp.Symbol("k") * dvec_len)
+ * sp.exp(sp.I * (m + l) * new_center_angle_rel_old_center)
+ * src_coeff_exprs[src_expansion.get_storage_index(m)]
+ for m in src_expansion.get_coefficient_identifiers()))
+ return translated_coeffs
+
+ raise RuntimeError("do not know how to translate %s to "
+ "local 2D Hankel expansion"
+ % type(src_expansion).__name__)
# }}}
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index c89e4c8b..c572e115 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -87,7 +87,7 @@ class VolumeTaylorMultipoleExpansion(
if not isinstance(src_expansion, type(self)):
raise RuntimeError("do not know how to translate %s to "
"Taylor multipole expansion"
- % type(src_expansion).__name)
+ % type(src_expansion).__name__)
logger.info("building translation operator: %s(%d) -> %s(%d): start"
% (type(src_expansion).__name__,
@@ -129,4 +129,63 @@ class VolumeTaylorMultipoleExpansion(
# }}}
+
+# {{{ 2D J-expansion
+
+class H2DMultipoleExpansion(MultipoleExpansionBase):
+ def __init__(self, kernel, order):
+ from sumpy.kernel import HelmholtzKernel
+ assert (isinstance(kernel.get_base_kernel(), HelmholtzKernel)
+ and kernel.dim == 2)
+
+ MultipoleExpansionBase.__init__(self, kernel, order)
+
+ def get_storage_index(self, k):
+ return self.order+k
+
+ def get_coefficient_identifiers(self):
+ return range(-self.order, self.order+1)
+
+ def coefficients_from_source(self, avec, bvec):
+ from sumpy.symbolic import sympy_real_norm_2
+ bessel_j = sp.Function("bessel_j")
+ avec_len = sympy_real_norm_2(avec)
+ # The coordinates are negated since avec points from source to center.
+ source_angle_rel_center = sp.atan2(-avec[1], -avec[0])
+ return [sp.Function("complex")(self.kernel.postprocess_at_source(
+ bessel_j(l, sp.Symbol("k") * avec_len) *
+ sp.exp(sp.I * l * -source_angle_rel_center), avec))
+ for l in self.get_coefficient_identifiers()]
+
+ def evaluate(self, coeffs, bvec):
+ from sumpy.symbolic import sympy_real_norm_2
+ hankel_1 = sp.Function("hankel_1")
+ bvec_len = sympy_real_norm_2(bvec)
+ target_angle_rel_center = sp.atan2(bvec[1], bvec[0])
+ return sum(coeffs[self.get_storage_index(l)]
+ * self.kernel.postprocess_at_target(
+ hankel_1(l, sp.Symbol("k") * bvec_len)
+ * sp.exp(sp.I * l * target_angle_rel_center), bvec)
+ for l in self.get_coefficient_identifiers())
+
+ def translate_from(self, src_expansion, src_coeff_exprs, dvec):
+ if not isinstance(src_expansion, type(self)):
+ raise RuntimeError("do not know how to translate %s to "
+ "multipole 2D Hankel expansion"
+ % type(src_expansion).__name__)
+ from sumpy.symbolic import sympy_real_norm_2
+ dvec_len = sympy_real_norm_2(dvec)
+ bessel_j = sp.Function("bessel_j")
+ new_center_angle_rel_old_center = sp.atan2(dvec[1], dvec[0])
+ translated_coeffs = []
+ for l in self.get_coefficient_identifiers():
+ translated_coeffs.append(
+ sum(src_coeff_exprs[src_expansion.get_storage_index(m)]
+ * bessel_j(m - l, sp.Symbol("k") * dvec_len)
+ * sp.exp(sp.I * (m - l) * new_center_angle_rel_old_center)
+ for m in src_expansion.get_coefficient_identifiers()))
+ return translated_coeffs
+
+# }}}
+
# vim: fdm=marker
diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index 371a3069..4fbbaedc 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -366,9 +366,10 @@ class HelmholtzKernel(ExpressionKernel):
return "HelmKnl(%s)" % (self.helmholtz_k_name)
def prepare_loopy_kernel(self, loopy_knl):
- from sumpy.codegen import bessel_preamble_generator, bessel_mangler
+ from sumpy.codegen import (bessel_preamble_generator, bessel_mangler,
+ complex_mangler)
loopy_knl = lp.register_function_manglers(loopy_knl,
- [bessel_mangler])
+ [bessel_mangler, complex_mangler])
loopy_knl = lp.register_preamble_generators(loopy_knl,
[bessel_preamble_generator])
diff --git a/test/test_fmm.py b/test/test_fmm.py
index 8ba2a34b..0726ca76 100644
--- a/test/test_fmm.py
+++ b/test/test_fmm.py
@@ -29,6 +29,11 @@ import pyopencl as cl
from pyopencl.tools import ( # noqa
pytest_generate_tests_for_pyopencl as pytest_generate_tests)
from sumpy.kernel import LaplaceKernel, HelmholtzKernel
+from sumpy.expansion.multipole import (
+ VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion)
+from sumpy.expansion.local import (
+ VolumeTaylorLocalExpansion, H2DLocalExpansion)
+
import pytest
import logging
@@ -43,13 +48,14 @@ else:
faulthandler.enable()
-@pytest.mark.parametrize("knl", [
- LaplaceKernel(2),
- LaplaceKernel(3),
- HelmholtzKernel(2),
- HelmholtzKernel(3),
+@pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class", [
+ (LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),
+ (LaplaceKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),
+ (HelmholtzKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),
+ (HelmholtzKernel(2), H2DLocalExpansion, H2DMultipoleExpansion),
+ (HelmholtzKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion)
])
-def test_sumpy_fmm(ctx_getter, knl):
+def test_sumpy_fmm(ctx_getter, knl, local_expn_class, mpole_expn_class):
logging.basicConfig(level=logging.INFO)
ctx = ctx_getter()
@@ -118,9 +124,6 @@ def test_sumpy_fmm(ctx_getter, knl):
logger.info("computing direct (reference) result")
- from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion
- from sumpy.expansion.local import VolumeTaylorLocalExpansion
-
from pytools.convergence import PConvergenceVerifier
pconv_verifier = PConvergenceVerifier()
@@ -136,8 +139,8 @@ def test_sumpy_fmm(ctx_getter, knl):
order_values = [1, 2]
for order in order_values:
- mpole_expn = VolumeTaylorMultipoleExpansion(knl, order)
- local_expn = VolumeTaylorLocalExpansion(knl, order)
+ mpole_expn = mpole_expn_class(knl, order)
+ local_expn = local_expn_class(knl, order)
out_kernels = [knl]
from sumpy.fmm import SumpyExpansionWranglerCodeContainer
diff --git a/test/test_kernels.py b/test/test_kernels.py
index ec1ab2d3..44861576 100644
--- a/test/test_kernels.py
+++ b/test/test_kernels.py
@@ -31,7 +31,8 @@ import pyopencl as cl
from pyopencl.tools import ( # noqa
pytest_generate_tests_for_pyopencl as pytest_generate_tests)
-from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion
+from sumpy.expansion.multipole import (
+ VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion)
from sumpy.expansion.local import VolumeTaylorLocalExpansion, H2DLocalExpansion
from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative,
DirectionalSourceDerivative)
@@ -95,6 +96,7 @@ def test_p2p(ctx_getter):
(HelmholtzKernel(2), VolumeTaylorMultipoleExpansion),
(HelmholtzKernel(2), VolumeTaylorLocalExpansion),
(HelmholtzKernel(2), H2DLocalExpansion),
+ (HelmholtzKernel(2), H2DMultipoleExpansion)
])
@pytest.mark.parametrize("with_source_derivative", [
False,
@@ -279,11 +281,12 @@ def test_p2e2p(ctx_getter, knl, expn_class, order, with_source_derivative):
assert eoc_rec_grad_x.order_estimate() > tgt_order_grad - grad_slack
-@pytest.mark.parametrize("knl", [
- LaplaceKernel(2),
- HelmholtzKernel(2)
+@pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class", [
+ (LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),
+ (HelmholtzKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),
+ (HelmholtzKernel(2), H2DLocalExpansion, H2DMultipoleExpansion)
])
-def test_translations(ctx_getter, knl):
+def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):
logging.basicConfig(level=logging.INFO)
from sympy.core.cache import clear_cache
@@ -336,7 +339,8 @@ def test_translations(ctx_getter, knl):
del eval_offset
- if isinstance(knl, HelmholtzKernel):
+ if isinstance(knl, HelmholtzKernel) and \
+ isinstance(local_expn_class, VolumeTaylorLocalExpansion):
# FIXME: Embarrassing--but we run out of memory for higher orders.
orders = [2, 3]
else:
@@ -366,8 +370,8 @@ def test_translations(ctx_getter, knl):
return pot
for order in orders:
- m_expn = VolumeTaylorMultipoleExpansion(knl, order=order)
- l_expn = VolumeTaylorLocalExpansion(knl, order=order)
+ m_expn = mpole_expn_class(knl, order=order)
+ l_expn = local_expn_class(knl, order=order)
from sumpy import P2EFromSingleBox, E2PFromSingleBox, P2P, E2EFromCSR
p2m = P2EFromSingleBox(ctx, m_expn)
--
GitLab