diff --git a/test/test_kernels.py b/test/test_kernels.py
index 2ddbd9c5f804269b3c91ae3cad323f6453fa2a65..d1768fd7c5040615a993f1ebb105b9bb1498c5b3 100644
--- a/test/test_kernels.py
+++ b/test/test_kernels.py
@@ -22,6 +22,7 @@ THE SOFTWARE.
import pytest
import sys
+from functools import partial
import numpy as np
import numpy.linalg as la
@@ -51,6 +52,8 @@ from sumpy.kernel import (
AxisTargetDerivative,
DirectionalSourceDerivative)
+import sumpy.toys as t
+
import logging
logger = logging.getLogger(__name__)
@@ -499,8 +502,6 @@ def test_translations(actx_factory, knl, local_expn_class, mpole_expn_class,
res = 20
nsources = 15
- target_kernels = [knl]
-
extra_kwargs = {}
if isinstance(knl, HelmholtzKernel):
extra_kwargs["k"] = 0.05
@@ -510,11 +511,10 @@ def test_translations(actx_factory, knl, local_expn_class, mpole_expn_class,
# Just to make sure things also work away from the origin
rng = np.random.default_rng(18)
origin = np.array([2, 1, 0][:knl.dim], np.float64)
- sources = actx.from_numpy(
+ sources = (
0.7 * (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64))
+ origin[:, np.newaxis])
- strengths = actx.from_numpy(
- np.ones(nsources, dtype=np.float64) * (1/nsources))
+ strengths = np.ones(nsources, dtype=np.float64) * (1/nsources)
pconv_verifier_p2m2p = PConvergenceVerifier()
pconv_verifier_p2m2m2p = PConvergenceVerifier()
@@ -525,8 +525,9 @@ def test_translations(actx_factory, knl, local_expn_class, mpole_expn_class,
eval_offset = np.array([5.5, 0.0, 0][:knl.dim])
fp = FieldPlotter(eval_offset + origin, extent=0.3, npoints=res)
+ targets = fp.points
- centers = actx.from_numpy((np.array(
+ centers = (np.array(
[
# box 0: particles, first mpole here
[0, 0, 0][:knl.dim],
@@ -540,7 +541,7 @@ def test_translations(actx_factory, knl, local_expn_class, mpole_expn_class,
# box 3: second local and eval here
eval_offset
],
- dtype=np.float64) + origin).T.copy())
+ dtype=np.float64) + origin).T.copy()
del eval_offset
@@ -549,192 +550,50 @@ def test_translations(actx_factory, knl, local_expn_class, mpole_expn_class,
else:
orders = [3, 4, 5]
- nboxes = centers.shape[-1]
-
- def eval_at(e2p, source_box_nr, rscale):
- e2p_target_boxes = actx.from_numpy(
- np.array([source_box_nr], dtype=np.int32))
-
- # These are indexed by global box numbers.
- e2p_box_target_starts = actx.from_numpy(
- np.array([0, 0, 0, 0], dtype=np.int32))
- e2p_box_target_counts_nonchild = actx.from_numpy(
- np.array([0, 0, 0, 0], dtype=np.int32))
- e2p_box_target_counts_nonchild[source_box_nr] = fp.points.shape[-1]
-
- evt, (pot,) = e2p(
- actx.queue,
-
- src_expansions=mpoles,
- src_base_ibox=0,
-
- target_boxes=e2p_target_boxes,
- box_target_starts=e2p_box_target_starts,
- box_target_counts_nonchild=e2p_box_target_counts_nonchild,
- centers=centers,
- targets=targets,
-
- rscale=rscale,
- **extra_kwargs
- )
- pot = actx.to_numpy(pot)
-
- return pot
-
m2l_factory = NonFFTM2LTranslationClassFactory()
m2l_translation = m2l_factory.get_m2l_translation_class(knl, local_expn_class)()
- for order in orders:
- m_expn = mpole_expn_class(knl, order=order)
- l_expn = local_expn_class(knl, order=order, m2l_translation=m2l_translation)
-
- from sumpy import P2EFromSingleBox, E2PFromSingleBox, P2P, E2EFromCSR
- p2m = P2EFromSingleBox(actx.context, m_expn)
- m2m = E2EFromCSR(actx.context, m_expn, m_expn)
- m2p = E2PFromSingleBox(actx.context, m_expn, target_kernels)
- m2l = E2EFromCSR(actx.context, m_expn, l_expn)
- l2l = E2EFromCSR(actx.context, l_expn, l_expn)
- l2p = E2PFromSingleBox(actx.context, l_expn, target_kernels)
- p2p = P2P(actx.context, target_kernels, exclude_self=False)
-
- targets = actx.from_numpy(make_obj_array(fp.points))
-
- # {{{ compute (direct) reference solution
-
- evt, (pot_direct,) = p2p(
- actx.queue,
- targets, sources, (strengths,),
- **extra_kwargs)
- pot_direct = actx.to_numpy(pot_direct)
-
- # }}}
-
- m1_rscale = 0.5
- m2_rscale = 0.25
- l1_rscale = 0.5
- l2_rscale = 0.25
-
- # {{{ apply P2M
-
- p2m_source_boxes = actx.from_numpy(np.array([0], dtype=np.int32))
+ toy_ctx = t.ToyContext(
+ actx.context,
+ kernel=knl,
+ local_expn_class=partial(local_expn_class,
+ m2l_translation=m2l_translation),
+ mpole_expn_class=mpole_expn_class,
+ extra_kernel_kwargs=extra_kwargs,
+ )
- # These are indexed by global box numbers.
- p2m_box_source_starts = actx.from_numpy(
- np.array([0, 0, 0, 0], dtype=np.int32))
- p2m_box_source_counts_nonchild = actx.from_numpy(
- np.array([nsources, 0, 0, 0], dtype=np.int32))
+ p = t.PointSources(toy_ctx, sources, weights=strengths)
+ p2p = p.eval(targets)
- evt, (mpoles,) = p2m(actx.queue,
- source_boxes=p2m_source_boxes,
- box_source_starts=p2m_box_source_starts,
- box_source_counts_nonchild=p2m_box_source_counts_nonchild,
- centers=centers,
- sources=sources,
- strengths=(strengths,),
- nboxes=nboxes,
- rscale=m1_rscale,
-
- tgt_base_ibox=0,
- **extra_kwargs)
-
- # }}}
-
- pot = eval_at(m2p, 0, m1_rscale)
-
- err = la.norm((pot - pot_direct) / res**2)
- err = err / (la.norm(pot_direct) / res**2)
+ m1_rscale = 0.5
+ m2_rscale = 0.25
+ l1_rscale = 0.5
+ l2_rscale = 0.25
+ for order in orders:
+ print(centers[:, 0].shape)
+ p2m = t.multipole_expand(p, centers[:, 0], order, m1_rscale)
+ p2m2p = p2m.eval(targets)
+ err = la.norm((p2m2p - p2p) / res**2)
+ err = err / (la.norm(p2p) / res**2)
pconv_verifier_p2m2p.add_data_point(order, err)
- # {{{ apply M2M
-
- m2m_target_boxes = actx.from_numpy(np.array([1], dtype=np.int32))
- m2m_src_box_starts = actx.from_numpy(np.array([0, 1], dtype=np.int32))
- m2m_src_box_lists = actx.from_numpy(np.array([0], dtype=np.int32))
-
- evt, (mpoles,) = m2m(actx.queue,
- src_expansions=mpoles,
- src_base_ibox=0,
- tgt_base_ibox=0,
- ntgt_level_boxes=mpoles.shape[0],
-
- target_boxes=m2m_target_boxes,
-
- src_box_starts=m2m_src_box_starts,
- src_box_lists=m2m_src_box_lists,
- centers=centers,
-
- src_rscale=m1_rscale,
- tgt_rscale=m2_rscale,
- **extra_kwargs)
-
- # }}}
-
- pot = eval_at(m2p, 1, m2_rscale)
-
- err = la.norm((pot - pot_direct)/res**2)
- err = err / (la.norm(pot_direct) / res**2)
-
+ p2m2m = t.multipole_expand(p2m, centers[:, 1], order, m2_rscale)
+ p2m2m2p = p2m2m.eval(targets)
+ err = la.norm((p2m2m2p - p2p)/res**2)
+ err = err / (la.norm(p2p) / res**2)
pconv_verifier_p2m2m2p.add_data_point(order, err)
- # {{{ apply M2L
-
- m2l_target_boxes = actx.from_numpy(np.array([2], dtype=np.int32))
- m2l_src_box_starts = actx.from_numpy(np.array([0, 1], dtype=np.int32))
- m2l_src_box_lists = actx.from_numpy(np.array([1], dtype=np.int32))
-
- evt, (mpoles,) = m2l(actx.queue,
- src_expansions=mpoles,
- src_base_ibox=0,
- tgt_base_ibox=0,
- ntgt_level_boxes=mpoles.shape[0],
-
- target_boxes=m2l_target_boxes,
- src_box_starts=m2l_src_box_starts,
- src_box_lists=m2l_src_box_lists,
- centers=centers,
-
- src_rscale=m2_rscale,
- tgt_rscale=l1_rscale,
- **extra_kwargs)
-
- # }}}
-
- pot = eval_at(l2p, 2, l1_rscale)
-
- err = la.norm((pot - pot_direct)/res**2)
- err = err / (la.norm(pot_direct) / res**2)
-
+ p2m2m2l = t.local_expand(p2m2m, centers[:, 2], order, l1_rscale)
+ p2m2m2l2p = p2m2m2l.eval(targets)
+ err = la.norm((p2m2m2l2p - p2p)/res**2)
+ err = err / (la.norm(p2p) / res**2)
pconv_verifier_p2m2m2l2p.add_data_point(order, err)
- # {{{ apply L2L
-
- l2l_target_boxes = actx.from_numpy(np.array([3], dtype=np.int32))
- l2l_src_box_starts = actx.from_numpy(np.array([0, 1], dtype=np.int32))
- l2l_src_box_lists = actx.from_numpy(np.array([2], dtype=np.int32))
-
- evt, (mpoles,) = l2l(actx.queue,
- src_expansions=mpoles,
- src_base_ibox=0,
- tgt_base_ibox=0,
- ntgt_level_boxes=mpoles.shape[0],
-
- target_boxes=l2l_target_boxes,
- src_box_starts=l2l_src_box_starts,
- src_box_lists=l2l_src_box_lists,
- centers=centers,
-
- src_rscale=l1_rscale,
- tgt_rscale=l2_rscale,
- **extra_kwargs)
-
- # }}}
-
- pot = eval_at(l2p, 3, l2_rscale)
-
- err = la.norm((pot - pot_direct)/res**2)
- err = err / (la.norm(pot_direct) / res**2)
-
+ p2m2m2l2l = t.local_expand(p2m2m2l, centers[:, 3], order, l2_rscale)
+ p2m2m2l2l2p = p2m2m2l2l.eval(targets)
+ err = la.norm((p2m2m2l2l2p - p2p)/res**2)
+ err = err / (la.norm(p2p) / res**2)
pconv_verifier_full.add_data_point(order, err)
for name, verifier in [