From 1799dfa95b39892ef42c5b228c232094a83e35a5 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 7 Feb 2019 16:37:27 -0600 Subject: [PATCH 1/2] Add a test of the toy infrastructure. Tests the most important functionality (P2E, E2P, E2P). Adapted from the 3D paper experiments. --- test/test_misc.py | 115 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 115 insertions(+) diff --git a/test/test_misc.py b/test/test_misc.py index 7b35b2dd..35b9137a 100644 --- a/test/test_misc.py +++ b/test/test_misc.py @@ -32,6 +32,8 @@ import pyopencl as cl # noqa: F401 from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) +from pytools import Record + from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, BiharmonicKernel, YukawaKernel) @@ -135,6 +137,119 @@ def test_order_finder(knl): print(orders) +# {{{ expansion toys p2e2e2p test cases + +def approx_convergence_factor(orders, errors): + poly = np.polyfit(orders, np.log(errors), deg=1) + return np.exp(poly[0]) + + +class P2E2E2PTestCase(Record): + + @property + def dim(self): + return len(self.source) + + @staticmethod + def eval(expr, source, center1, center2, target): + from pymbolic import parse, evaluate + context = { + "s": source, + "c1": center1, + "c2": center2, + "t": target, + "norm": la.norm} + + return evaluate(parse(expr), context) + + def __init__(self, + source, center1, center2, target, expansion1, expansion2, conv_factor): + + if isinstance(conv_factor, str): + conv_factor = self.eval(conv_factor, source, center1, center2, target) + + Record.__init__(self, + source=source, + center1=center1, + center2=center2, + target=target, + expansion1=expansion1, + expansion2=expansion2, + conv_factor=conv_factor) + + +P2E2E2P_TEST_CASES = ( + # local to local, 3D + P2E2E2PTestCase( + source=np.array([3., 4., 5.]), + center1=np.array([1., 0., 0.]), + center2=np.array([1., 3., 0.]), + target=np.array([1., 1., 1.]), + expansion1=t.local_expand, + expansion2=t.local_expand, + conv_factor="norm(t-c1)/norm(s-c1)"), + # multipole to multipole, 3D + P2E2E2PTestCase( + source=np.array([1., 1., 1.]), + center1=np.array([1., 0., 0.]), + center2=np.array([1., 0., 3.]), + target=np.array([3., 4., 5.]), + expansion1=t.multipole_expand, + expansion2=t.multipole_expand, + conv_factor="norm(s-c2)/norm(t-c2)"), + # multipole to local, 3D + P2E2E2PTestCase( + source=np.array([-2., 2., 1.]), + center1=np.array([-2., 5., 3.]), + center2=np.array([0., 0., 0.]), + target=np.array([0., 0., -1]), + expansion1=t.multipole_expand, + expansion2=t.local_expand, + conv_factor="norm(t-c2)/(norm(c2-c1)-norm(c1-s))"), +) + +# }}} + + +ORDERS_P2E2E2P = (3, 4, 5) +RTOL_P2E2E2P = 1e-2 + + +@pytest.mark.parametrize("case", P2E2E2P_TEST_CASES) +def test_toy_p2e2e2p(ctx_getter, case): + dim = case.dim + + src = case.source.reshape(dim, -1) + tgt = case.target.reshape(dim, -1) + + if not 0 <= case.conv_factor <= 1: + raise ValueError( + "convergence factor not in valid range: %e" % case.conv_factor) + + from sumpy.expansion.local import VolumeTaylorLocalExpansion + from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion + + cl_ctx = ctx_getter() + ctx = t.ToyContext(cl_ctx, + LaplaceKernel(dim), + VolumeTaylorMultipoleExpansion, + VolumeTaylorLocalExpansion) + + errors = [] + + src_pot = t.PointSources(ctx, src, weights=np.array([1.])) + pot_actual = src_pot.eval(tgt).item() + + for order in ORDERS_P2E2E2P: + expn = case.expansion1(src_pot, case.center1, order=order) + expn2 = case.expansion2(expn, case.center2, order=order) + pot_p2e2e2p = expn2.eval(tgt).item() + errors.append(np.abs(pot_actual - pot_p2e2e2p)) + + conv_factor = approx_convergence_factor(ORDERS_P2E2E2P, errors) + assert conv_factor < min(1, case.conv_factor * (1 + RTOL_P2E2E2P)) + + # You can test individual routines by typing # $ python test_misc.py 'test_p2p(cl.create_some_context)' -- GitLab From 762c398c1d872da26910b54bd3cb08e28094a48d Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 11 Feb 2019 11:03:26 -0600 Subject: [PATCH 2/2] A few small changes --- test/test_misc.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/test_misc.py b/test/test_misc.py index 35b9137a..33e8535a 100644 --- a/test/test_misc.py +++ b/test/test_misc.py @@ -188,6 +188,7 @@ P2E2E2P_TEST_CASES = ( expansion1=t.local_expand, expansion2=t.local_expand, conv_factor="norm(t-c1)/norm(s-c1)"), + # multipole to multipole, 3D P2E2E2PTestCase( source=np.array([1., 1., 1.]), @@ -197,6 +198,7 @@ P2E2E2P_TEST_CASES = ( expansion1=t.multipole_expand, expansion2=t.multipole_expand, conv_factor="norm(s-c2)/norm(t-c2)"), + # multipole to local, 3D P2E2E2PTestCase( source=np.array([-2., 2., 1.]), @@ -246,8 +248,9 @@ def test_toy_p2e2e2p(ctx_getter, case): pot_p2e2e2p = expn2.eval(tgt).item() errors.append(np.abs(pot_actual - pot_p2e2e2p)) - conv_factor = approx_convergence_factor(ORDERS_P2E2E2P, errors) - assert conv_factor < min(1, case.conv_factor * (1 + RTOL_P2E2E2P)) + conv_factor = approx_convergence_factor(1 + np.array(ORDERS_P2E2E2P), errors) + assert conv_factor <= min(1, case.conv_factor * (1 + RTOL_P2E2E2P)), \ + (conv_factor, case.conv_factor * (1 + RTOL_P2E2E2P)) # You can test individual routines by typing -- GitLab