From c91c1498b692a202ddbd172f628ef4435b6f1775 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Sun, 21 Jul 2013 22:47:51 -0400
Subject: [PATCH] Build, test M2M
---
sumpy/expansion/multipole.py | 56 ++++++++++++++++++++++++++++++++----
test/test_kernels.py | 22 +++++++++++++-
2 files changed, 72 insertions(+), 6 deletions(-)
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index 924055dc..425980e9 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -26,6 +26,9 @@ import sympy as sp # noqa
from sumpy.expansion import ExpansionBase, VolumeTaylorExpansionBase
+import logging
+logger = logging.getLogger(__name__)
+
class MultipoleExpansionBase(ExpansionBase):
pass
@@ -61,26 +64,69 @@ class VolumeTaylorMultipoleExpansion(
for iaxis, mi_i in enumerate(mi))
result[i] += (
- - 1/mi_factorial(derivative_mi)
- * mi_power(avec, derivative_mi)
+ - mi_power(avec, derivative_mi)
* dir_vec[idim])
return result
else:
return [
- mi_power(avec, mi) / mi_factorial(mi)
+ mi_power(avec, mi)
for mi in self.get_coefficient_identifiers()]
def evaluate(self, coeffs, bvec):
ppkernel = self.kernel.postprocess_at_target(
self.kernel.get_expression(bvec), bvec)
- from sumpy.tools import mi_derivative
+ from sumpy.tools import mi_derivative, mi_factorial
return sum(
- coeff
+ coeff / mi_factorial(mi)
* mi_derivative(ppkernel, bvec, mi)
for coeff, mi in zip(coeffs, 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 "
+ "Taylor multipole expansion"
+ % type(src_expansion).__name)
+
+ logger.info("building translation operator: %s(%d) -> %s(%d): start"
+ % (type(src_expansion).__name__,
+ src_expansion.order,
+ type(self).__name__,
+ self.order))
+
+ src_mi_to_index = dict((mi, i) for i, mi in enumerate(
+ src_expansion.get_coefficient_identifiers()))
+
+ result = [0] * len(self.get_coefficient_identifiers())
+ from pytools import generate_nonnegative_integer_tuples_below as gnitb
+
+ for i, tgt_mi in enumerate(
+ self.get_coefficient_identifiers()):
+
+ tgt_mi_plus_one = tuple(mi_i + 1 for mi_i in tgt_mi)
+
+ for src_mi in gnitb(tgt_mi_plus_one):
+ try:
+ src_index = src_mi_to_index[src_mi]
+ except KeyError:
+ # tgt higher-order than source--dumb, but not life-threatening
+ continue
+
+ contrib = src_coeff_exprs[src_index]
+
+ for idim in xrange(self.dim):
+ n = tgt_mi[idim]
+ k = src_mi[idim]
+ assert n >= k
+ contrib *= (sp.binomial(n, k)
+ * dvec[idim]**(n-k))
+
+ result[i] += contrib
+
+ logger.info("building translation operator: done")
+ return result
+
# }}}
# vim: fdm=marker
diff --git a/test/test_kernels.py b/test/test_kernels.py
index d74c45ba..c4f08b1d 100644
--- a/test/test_kernels.py
+++ b/test/test_kernels.py
@@ -365,6 +365,7 @@ def test_translations(ctx_getter, knl):
from sumpy import P2E, E2P, P2P, E2E
p2m = P2E(ctx, m_expn, out_kernels)
+ m2m = E2E(ctx, m_expn, m_expn)
m2l = E2E(ctx, m_expn, l_expn)
l2l = E2E(ctx, l_expn, l_expn)
l2p = E2P(ctx, l_expn, out_kernels)
@@ -391,11 +392,28 @@ def test_translations(ctx_getter, knl):
# }}}
+ # {{{ apply M2M
+
+ m2m_target_boxes = np.array([1], dtype=np.int32)
+ m2m_src_box_starts = np.array([0, 1], dtype=np.int32)
+ m2m_src_box_lists = np.array([0], dtype=np.int32)
+
+ evt, (mpoles,) = m2m(queue,
+ src_expansions=mpoles,
+ target_boxes=m2m_target_boxes,
+ src_box_starts=m2m_src_box_starts,
+ src_box_lists=m2m_src_box_lists,
+ centers=centers,
+ #flags="print_hl_cl",
+ out_host=True, **extra_kwargs)
+
+ # }}}
+
# {{{ apply M2L
m2l_target_boxes = np.array([2], dtype=np.int32)
m2l_src_box_starts = np.array([0, 1], dtype=np.int32)
- m2l_src_box_lists = np.array([0], dtype=np.int32)
+ m2l_src_box_lists = np.array([1], dtype=np.int32)
evt, (mpoles,) = m2l(queue,
src_expansions=mpoles,
@@ -456,6 +474,8 @@ def test_translations(ctx_getter, knl):
err = la.norm((pot - pot_direct)/res**2)
+ err = err / (la.norm(pot_direct) / res**2)
+
# }}}
pconv_verifier.add_data_point(order, err)
--
GitLab