From 210e4931d1425dc03c14254c6432a27502c73913 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 20 Aug 2022 20:01:30 +0300 Subject: [PATCH 1/8] add array context skeleton --- requirements.txt | 1 + setup.py | 1 + sumpy/array_context.py | 77 +++++++++++++ sumpy/expansion/multipole.py | 2 +- sumpy/qbx.py | 2 +- test/test_codegen.py | 15 ++- test/test_cse.py | 213 +++++++++++++++++++++++++---------- 7 files changed, 245 insertions(+), 66 deletions(-) create mode 100644 sumpy/array_context.py diff --git a/requirements.txt b/requirements.txt index a26c33e3b..ff9ca871b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -8,4 +8,5 @@ git+https://github.com/inducer/islpy.git#egg=islpy git+https://github.com/inducer/pyopencl.git#egg=pyopencl git+https://github.com/inducer/boxtree.git#egg=boxtree git+https://github.com/inducer/loopy.git#egg=loopy +git+https://github.com/inducer/arraycontext.git#egg=arraycontext git+https://github.com/inducer/pyfmmlib.git#egg=pyfmmlib diff --git a/setup.py b/setup.py index f9f0be9be..19ca094fc 100644 --- a/setup.py +++ b/setup.py @@ -101,6 +101,7 @@ def write_git_revision(package_name): "pytools>=2021.1.1", "loopy>=2021.1", "boxtree>=2018.1", + "arraycontext", "pytest>=2.3", "pyrsistent>=0.16.0", "dataclasses>=0.7;python_version<='3.6'", diff --git a/sumpy/array_context.py b/sumpy/array_context.py new file mode 100644 index 000000000..760f7d5d8 --- /dev/null +++ b/sumpy/array_context.py @@ -0,0 +1,77 @@ +__copyright__ = "Copyright (C) 2022 Alexandru Fikl" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from boxtree.array_context import PyOpenCLArrayContext as PyOpenCLArrayContextBase +from arraycontext.pytest import ( + _PytestPyOpenCLArrayContextFactoryWithClass, + register_pytest_array_context_factory) + +__doc__ = """ +.. autoclass:: PyOpenCLArrayContext +""" + + +# {{{ PyOpenCLArrayContext + +class PyOpenCLArrayContext(PyOpenCLArrayContextBase): + def transform_loopy_program(self, t_unit): + default_ep = t_unit.default_entrypoint + options = default_ep.options + + if not (options.return_dict and options.no_numpy): + raise ValueError("Loopy kernel passed to call_loopy must " + "have return_dict and no_numpy options set. " + "Did you use arraycontext.make_loopy_program " + "to create this kernel?") + + return super().transform_loopy_program(t_unit) + +# }}} + + +# {{{ pytest + +def _acf(): + import pyopencl as cl + ctx = cl.create_some_context() + queue = cl.CommandQueue(ctx) + + return PyOpenCLArrayContext(queue, force_device_scalars=True) + + +class PytestPyOpenCLArrayContextFactory( + _PytestPyOpenCLArrayContextFactoryWithClass): + actx_class = PyOpenCLArrayContext + + def __call__(self): + # NOTE: prevent any cache explosions during testing! + from sympy.core.cache import clear_cache + clear_cache() + + return super().__call__() + + +register_pytest_array_context_factory( + "sumpy.pyopencl", + PytestPyOpenCLArrayContextFactory) + +# }}} diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index fd5a3b1d1..1924e6c7f 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -203,7 +203,7 @@ def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, # └───⬏ ↑ # └─────┘ # - # For the second hyperplane, data is propogated rightwards first + # For the second hyperplane, data is propagated rightwards first # and then upwards second which is opposite to that of the first # hyperplane. # diff --git a/sumpy/qbx.py b/sumpy/qbx.py index b479e1db3..de2598881 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -104,7 +104,7 @@ def _evaluate(self, sac, avec, bvec, rscale, expansion_nr, coefficients): # In LineTaylorLocalExpansion.evaluate, we can't run # postprocess_at_target because the coefficients are assigned # symbols and postprocess with a derivative will make them zero. - # Instead run postprocess here before the coeffients are assigned. + # Instead run postprocess here before the coefficients are assigned. coefficients = [tgt_knl.postprocess_at_target(coeff, bvec) for coeff in coefficients] diff --git a/test/test_codegen.py b/test/test_codegen.py index 815b92dd9..9c092e58e 100644 --- a/test/test_codegen.py +++ b/test/test_codegen.py @@ -20,13 +20,15 @@ THE SOFTWARE. """ - +import pytest import sys import logging logger = logging.getLogger(__name__) +# {{{ test_symbolic_assignment_name_uniqueness + def test_symbolic_assignment_name_uniqueness(): # https://gitlab.tiker.net/inducer/sumpy/issues/13 from sumpy.assignment_collection import SymbolicAssignmentCollection @@ -43,6 +45,10 @@ def test_symbolic_assignment_name_uniqueness(): assert len(sac.assignments) == 3 +# }}} + + +# {{{ test_line_taylor_coeff_growth def test_line_taylor_coeff_growth(): # Regression test for LineTaylorLocalExpansion. @@ -70,15 +76,16 @@ def test_line_taylor_coeff_growth(): max_order = 2 assert np.polyfit(np.log(indices), np.log(counts), deg=1)[0] < max_order +# }}} + # You can test individual routines by typing -# $ python test_fmm.py "test_sumpy_fmm(cl.create_some_context)" +# $ python test_codegen.py 'test_line_taylor_coeff_growth()' if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from pytest import main - main([__file__]) + pytest.main([__file__]) # vim: fdm=marker diff --git a/test/test_cse.py b/test/test_cse.py index ba9986bed..68352c905 100644 --- a/test/test_cse.py +++ b/test/test_cse.py @@ -67,26 +67,22 @@ import pytest import sys -from sumpy.symbolic import ( - Add, Pow, exp, sqrt, symbols, sympify, cos, sin, Function, USE_SYMENGINE) +import sumpy.symbolic as sym +from sumpy.cse import cse, preprocess_for_cse, postprocess_for_cse -if not USE_SYMENGINE: +if not sym.USE_SYMENGINE: from sympy.simplify.cse_opts import sub_pre, sub_post from sympy.functions.special.hyper import meijerg from sympy.simplify import cse_opts -from sumpy.cse import ( - cse, preprocess_for_cse, postprocess_for_cse) +import logging +logger = logging.getLogger(__name__) - -w, x, y, z = symbols("w,x,y,z") -x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12 = symbols("x:13") +w, x, y, z = sym.symbols("w,x,y,z") +x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12 = sym.symbols("x:13") sympyonly = ( - pytest.mark.skipif(USE_SYMENGINE, reason="uses a sympy-only feature")) - - -# Dummy "optimization" functions for testing. + pytest.mark.skipif(sym.USE_SYMENGINE, reason="uses a sympy-only feature")) def opt1(expr): @@ -97,6 +93,8 @@ def opt2(expr): return expr*z +# {{{ test_preprocess_for_cse + def test_preprocess_for_cse(): assert preprocess_for_cse(x, [(opt1, None)]) == x + y assert preprocess_for_cse(x, [(None, opt1)]) == x @@ -105,6 +103,10 @@ def test_preprocess_for_cse(): assert preprocess_for_cse( x, [(opt1, None), (opt2, None)]) == (x + y)*z +# }}} + + +# {{{ test_postprocess_for_cse def test_postprocess_for_cse(): assert postprocess_for_cse(x, [(opt1, None)]) == x @@ -115,19 +117,27 @@ def test_postprocess_for_cse(): assert postprocess_for_cse( x, [(None, opt1), (None, opt2)]) == x*z + y +# }}} + + +# {{{ test_cse_single def test_cse_single(): # Simple substitution. - e = Add(Pow(x + y, 2), sqrt(x + y)) + e = sym.Add(sym.Pow(x + y, 2), sym.sqrt(x + y)) substs, reduced = cse([e]) assert substs == [(x0, x + y)] - assert reduced == [sqrt(x0) + x0**2] + assert reduced == [sym.sqrt(x0) + x0**2] + +# }}} + +# {{{ @sympyonly def test_cse_not_possible(): # No substitution possible. - e = Add(x, y) + e = sym.Add(x, y) substs, reduced = cse([e]) assert substs == [] assert reduced == [x + y] @@ -136,34 +146,46 @@ def test_cse_not_possible(): + meijerg((1, 3), (y, 4), (5,), [], x)) assert cse(eq) == ([], [eq]) +# }}} + + +# {{{ test_nested_substitution def test_nested_substitution(): # Substitution within a substitution. - e = Add(Pow(w*x + y, 2), sqrt(w*x + y)) + e = sym.Add(sym.Pow(w*x + y, 2), sym.sqrt(w*x + y)) substs, reduced = cse([e]) assert substs == [(x0, w*x + y)] - assert reduced == [sqrt(x0) + x0**2] + assert reduced == [sym.sqrt(x0) + x0**2] +# }}} + + +# {{{ test_subtraction_opt @sympyonly def test_subtraction_opt(): # Make sure subtraction is optimized. - e = (x - y)*(z - y) + exp((x - y)*(z - y)) + e = (x - y)*(z - y) + sym.exp((x - y)*(z - y)) substs, reduced = cse( [e], optimizations=[(cse_opts.sub_pre, cse_opts.sub_post)]) assert substs == [(x0, (x - y)*(y - z))] - assert reduced == [-x0 + exp(-x0)] - e = -(x - y)*(z - y) + exp(-(x - y)*(z - y)) + assert reduced == [-x0 + sym.exp(-x0)] + e = -(x - y)*(z - y) + sym.exp(-(x - y)*(z - y)) substs, reduced = cse( [e], optimizations=[(cse_opts.sub_pre, cse_opts.sub_post)]) assert substs == [(x0, (x - y)*(y - z))] - assert reduced == [x0 + exp(x0)] + assert reduced == [x0 + sym.exp(x0)] # issue 4077 n = -1 + 1/x e = n/x/(-n)**2 - 1/n/x assert cse(e, optimizations=[(cse_opts.sub_pre, cse_opts.sub_post)]) == \ ([], [0]) +# }}} + + +# {{{ test_multiple_expressions def test_multiple_expressions(): e1 = (x + y)*z @@ -181,7 +203,7 @@ def test_multiple_expressions(): rsubsts, _ = cse(reversed(l_)) assert substs == rsubsts assert reduced == [x1, x1 + z, x0] - f = Function("f") + f = sym.Function("f") l_ = [f(x - z, y - z), x - z, y - z] substs, reduced = cse(l_) rsubsts, _ = cse(reversed(l_)) @@ -195,31 +217,46 @@ def test_multiple_expressions(): assert cse([x*y, z + x*y, x*y*z + 3]) == \ ([(x0, x*y)], [x0, z + x0, 3 + x0*z]) +# }}} + + +# {{{ test_issue_4203 def test_issue_4203(): - assert cse(sin(x**x)/x**x) == ([(x0, x**x)], [sin(x0)/x0]) + assert cse(sym.sin(x**x)/x**x) == ([(x0, x**x)], [sym.sin(x0)/x0]) + +# }}} +# {{{ test_dont_cse_subs + def test_dont_cse_subs(): - f = Function("f") - g = Function("g") + f = sym.Function("f") + g = sym.Function("g") name_val, (expr,) = cse(f(x+y).diff(x) + g(x+y).diff(x)) assert name_val == [] +# }}} + + +# {{{ test_dont_cse_derivative def test_dont_cse_derivative(): - from sumpy.symbolic import Derivative - f = Function("f") + f = sym.Function("f") - deriv = Derivative(f(x+y), x) + deriv = sym.Derivative(f(x+y), x) name_val, (expr,) = cse(x + y + deriv) assert name_val == [] assert expr == x + y + deriv +# }}} + + +# {{{ test_pow_invpow def test_pow_invpow(): assert cse(1/x**2 + x**2) == \ @@ -228,39 +265,64 @@ def test_pow_invpow(): ([(x0, x**2), (x1, 1/x0)], [x0 + x1*(x1 + 1)]) assert cse(1/x**2 + (1 + 1/x**2)*x**2) == \ ([(x0, x**2), (x1, 1/x0)], [x0*(x1 + 1) + x1]) - assert cse(cos(1/x**2) + sin(1/x**2)) == \ - ([(x0, x**(-2))], [sin(x0) + cos(x0)]) - assert cse(cos(x**2) + sin(x**2)) == \ - ([(x0, x**2)], [sin(x0) + cos(x0)]) + assert cse(sym.cos(1/x**2) + sym.sin(1/x**2)) == \ + ([(x0, x**(-2))], [sym.sin(x0) + sym.cos(x0)]) + assert cse(sym.cos(x**2) + sym.sin(x**2)) == \ + ([(x0, x**2)], [sym.sin(x0) + sym.cos(x0)]) assert cse(y/(2 + x**2) + z/x**2/y) == \ ([(x0, x**2)], [y/(x0 + 2) + z/(x0*y)]) - assert cse(exp(x**2) + x**2*cos(1/x**2)) == \ - ([(x0, x**2)], [x0*cos(1/x0) + exp(x0)]) + assert cse(sym.exp(x**2) + x**2*sym.cos(1/x**2)) == \ + ([(x0, x**2)], [x0*sym.cos(1/x0) + sym.exp(x0)]) assert cse((1 + 1/x**2)/x**2) == \ ([(x0, x**(-2))], [x0*(x0 + 1)]) assert cse(x**(2*y) + x**(-2*y)) == \ ([(x0, x**(2*y))], [x0 + 1/x0]) +# }}} + + +# {{{ test_issue_4499 @sympyonly def test_issue_4499(): # previously, this gave 16 constants from sympy.abc import a, b from sympy import Tuple, S - B = Function("B") # noqa - G = Function("G") # noqa - t = Tuple( - *(a, a + S(1)/2, 2*a, b, 2*a - b + 1, (sqrt(z)/2)**(-2*a + 1)*B(2*a - - b, sqrt(z))*B(b - 1, sqrt(z))*G(b)*G(2*a - b + 1), # noqa - sqrt(z)*(sqrt(z)/2)**(-2*a + 1)*B(b, sqrt(z))*B(2*a - b, # noqa - sqrt(z))*G(b)*G(2*a - b + 1), sqrt(z)*(sqrt(z)/2)**(-2*a + 1)*B(b - 1, - sqrt(z))*B(2*a - b + 1, sqrt(z))*G(b)*G(2*a - b + 1), - (sqrt(z)/2)**(-2*a + 1)*B(b, sqrt(z))*B(2*a - b + 1, # noqa - sqrt(z))*G(b)*G(2*a - b + 1), 1, 0, S(1)/2, z/2, -b + 1, -2*a + b, # noqa - -2*a)) # noqa + + B = sym.Function("B") # noqa: N806 + G = sym.Function("G") # noqa: N806 + t = Tuple(*( + a, + a + S(1)/2, + 2*a, + b, + 2*a - b + 1, + (sym.sqrt(z)/2)**(-2*a + 1) + * B(2*a-b, sym.sqrt(z)) + * B(b - 1, sym.sqrt(z))*G(b)*G(2*a - b + 1), + sym.sqrt(z)*(sym.sqrt(z)/2)**(-2*a + 1) + * B(b, sym.sqrt(z)) + * B(2*a - b, sym.sqrt(z))*G(b)*G(2*a - b + 1), + sym.sqrt(z)*(sym.sqrt(z)/2)**(-2*a + 1) + * B(b - 1, sym.sqrt(z)) + * B(2*a - b + 1, sym.sqrt(z))*G(b)*G(2*a - b + 1), + (sym.sqrt(z)/2)**(-2*a + 1) + * B(b, sym.sqrt(z)) + * B(2*a - b + 1, sym.sqrt(z))*G(b)*G(2*a - b + 1), + 1, + 0, + S(1)/2, + z/2, + -b + 1, + -2*a + b, + -2*a)) c = cse(t) assert len(c[0]) == 11 +# }}} + + +# {{{ test_issue_6169 @sympyonly def test_issue_6169(): @@ -271,14 +333,17 @@ def test_issue_6169(): # mechanism assert sub_post(sub_pre((-x - y)*z - x - y)) == -z*(x + y) - x - y +# }}} + + +# {{{ test_cse_Indexed @sympyonly -def test_cse_Indexed(): # noqa +def test_cse_indexed(): from sympy import IndexedBase, Idx len_y = 5 y = IndexedBase("y", shape=(len_y,)) x = IndexedBase("x", shape=(len_y,)) - Dy = IndexedBase("Dy", shape=(len_y-1,)) # noqa i = Idx("i", len_y-1) expr1 = (y[i+1]-y[i])/(x[i+1]-x[i]) @@ -286,9 +351,13 @@ def test_cse_Indexed(): # noqa replacements, reduced_exprs = cse([expr1, expr2]) assert len(replacements) > 0 +# }}} + + +# {{{ test_Piecewise @sympyonly -def test_Piecewise(): # noqa +def test_piecewise(): from sympy import Piecewise, Eq f = Piecewise((-z + x*y, Eq(y, 0)), (-z - x*y, True)) ans = cse(f) @@ -296,41 +365,57 @@ def test_Piecewise(): # noqa [Piecewise((x0+x1, Eq(y, 0)), (x0 - x1, True))]) assert ans == actual_ans +# }}} + + +# {{{ test_name_conflict def test_name_conflict(): z1 = x0 + y z2 = x2 + x3 - l_ = [cos(z1) + z1, cos(z2) + z2, x0 + x2] + l_ = [sym.cos(z1) + z1, sym.cos(z2) + z2, x0 + x2] substs, reduced = cse(l_) assert [e.subs(dict(substs)) for e in reduced] == l_ +# }}} + + +# {{{ test_name_conflict_cust_symbols def test_name_conflict_cust_symbols(): z1 = x0 + y z2 = x2 + x3 - l_ = [cos(z1) + z1, cos(z2) + z2, x0 + x2] - substs, reduced = cse(l_, symbols("x:10")) + l_ = [sym.cos(z1) + z1, sym.cos(z2) + z2, x0 + x2] + substs, reduced = cse(l_, sym.symbols("x:10")) assert [e.subs(dict(substs)) for e in reduced] == l_ +# }}} + + +# {{{ test_symbols_exhausted_error def test_symbols_exhausted_error(): - l_ = cos(x+y)+x+y+cos(w+y)+sin(w+y) - sym = [x, y, z] + l_ = sym.cos(x+y)+x+y+sym.cos(w+y)+sym.sin(w+y) + s = [x, y, z] with pytest.raises(ValueError): - print(cse(l_, symbols=sym)) + logger.info("cse:\n%s", cse(l_, symbols=s)) + +# }}} + +# {{{ test_issue_7840 @sympyonly def test_issue_7840(): # daveknippers' example - C393 = sympify( # noqa + C393 = sym.sympify( # noqa: N806 "Piecewise((C391 - 1.65, C390 < 0.5), (Piecewise((C391 - 1.65, \ C391 > 2.35), (C392, True)), True))" ) - C391 = sympify( # noqa + C391 = sym.sympify( # noqa: N806 "Piecewise((2.05*C390**(-1.03), C390 < 0.5), (2.5*C390**(-0.625), True))" ) - C393 = C393.subs("C391",C391) # noqa + C393 = C393.subs("C391", C391) # noqa: N806 # simple substitution sub = {} sub["C390"] = 0.703451854 @@ -345,7 +430,7 @@ def test_issue_7840(): assert ss_answer == cse_answer # GitRay's example - expr = sympify( + expr = sym.sympify( "Piecewise((Symbol('ON'), Equality(Symbol('mode'), Symbol('ON'))), \ (Piecewise((Piecewise((Symbol('OFF'), StrictLessThan(Symbol('x'), \ Symbol('threshold'))), (Symbol('ON'), S.true)), Equality(Symbol('mode'), \ @@ -357,6 +442,10 @@ def test_issue_7840(): # there should not be any replacements assert len(substitutions) < 1 +# }}} + + +# {{{ test_recursive_matching def test_recursive_matching(): assert cse([x+y, 2+x+y, x+y+z, 3+x+y+z]) == \ @@ -372,12 +461,16 @@ def test_recursive_matching(): assert cse([2*x*x, x*x*y, x*x*y*w, x*x*y*w*x0, x*x*y*w*x2]) == \ ([(x1, x**2), (x3, x1*y), (x4, w*x3)], [2*x1, x3, x4, x0*x4, x2*x4]) +# }}} + + +# You can test individual routines by typing +# $ python test_cse.py 'test_recursive_matching()' if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from pytest import main - main([__file__]) + pytest.main([__file__]) # vim: fdm=marker From 3f4df12c9573ebeb843e88719c3a5d7368de3329 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 20 Aug 2022 20:01:42 +0300 Subject: [PATCH 2/8] port test_fmm to arraycontext --- test/test_fmm.py | 369 ++++++++++++++++++++++++----------------------- 1 file changed, 185 insertions(+), 184 deletions(-) diff --git a/test/test_fmm.py b/test/test_fmm.py index 58ec8b3ad..f447999cf 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -20,42 +20,47 @@ THE SOFTWARE. """ - +import pytest import sys import os -from unittest.mock import patch +from functools import partial + import numpy as np import numpy.linalg as la -import pyopencl as cl -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) -from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, YukawaKernel, + +from arraycontext import pytest_generate_tests_for_array_contexts +from sumpy.array_context import ( # noqa: F401 + PytestPyOpenCLArrayContextFactory, _acf) + +from sumpy.kernel import ( + LaplaceKernel, + HelmholtzKernel, + YukawaKernel, BiharmonicKernel) from sumpy.expansion.multipole import ( VolumeTaylorMultipoleExpansion, - H2DMultipoleExpansion, Y2DMultipoleExpansion, + H2DMultipoleExpansion, + Y2DMultipoleExpansion, LinearPDEConformingVolumeTaylorMultipoleExpansion) from sumpy.expansion.local import ( VolumeTaylorLocalExpansion, - H2DLocalExpansion, Y2DLocalExpansion, + H2DLocalExpansion, + Y2DLocalExpansion, LinearPDEConformingVolumeTaylorLocalExpansion) from sumpy.fmm import ( - SumpyTreeIndependentDataForWrangler, - SumpyExpansionWrangler) + SumpyTreeIndependentDataForWrangler, + SumpyExpansionWrangler) -import pytest import logging logger = logging.getLogger(__name__) +pytest_generate_tests = pytest_generate_tests_for_array_contexts([ + PytestPyOpenCLArrayContextFactory, + ]) -try: - import faulthandler -except ImportError: - pass -else: - faulthandler.enable() +# {{{ test_sumpy_fmm @pytest.mark.parametrize("use_translation_classes, use_fft, fft_backend", [(False, False, None), (True, False, None), (True, True, "loopy"), @@ -84,10 +89,11 @@ (YukawaKernel(2), Y2DLocalExpansion, Y2DMultipoleExpansion, False), ]) -def test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class, +def test_sumpy_fmm(actx_factory, knl, local_expn_class, mpole_expn_class, order_varies_with_level, use_translation_classes, use_fft, - fft_backend): - logging.basicConfig(level=logging.INFO) + fft_backend, visualize=False): + if visualize: + logging.basicConfig(level=logging.INFO) if local_expn_class == VolumeTaylorLocalExpansion and use_fft: pytest.skip("VolumeTaylorExpansion with FFT takes a lot of resources.") @@ -96,68 +102,66 @@ def test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class, pytest.skip("Fourier/Bessel based expansions with FFT is not supported yet.") if use_fft: + from unittest.mock import patch + with patch.dict(os.environ, {"SUMPY_FFT_BACKEND": fft_backend}): - _test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class, + _test_sumpy_fmm(actx_factory, knl, local_expn_class, mpole_expn_class, order_varies_with_level, use_translation_classes, use_fft, fft_backend) else: - _test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class, + _test_sumpy_fmm(actx_factory, knl, local_expn_class, mpole_expn_class, order_varies_with_level, use_translation_classes, use_fft, fft_backend) -def _test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class, +def _test_sumpy_fmm(actx_factory, knl, local_expn_class, mpole_expn_class, order_varies_with_level, use_translation_classes, use_fft, fft_backend): - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) + actx = actx_factory() nsources = 1000 ntargets = 300 dtype = np.float64 - from boxtree.tools import ( - make_normal_particle_array as p_normal) + from boxtree.tools import make_normal_particle_array as p_normal + sources = p_normal(actx.queue, nsources, knl.dim, dtype, seed=15) - sources = p_normal(queue, nsources, knl.dim, dtype, seed=15) if 1: offset = np.zeros(knl.dim) offset[0] = 0.1 - - targets = ( - p_normal(queue, ntargets, knl.dim, dtype, seed=18) - + offset) + targets = offset + p_normal(actx.queue, ntargets, knl.dim, dtype, seed=18) del offset else: from sumpy.visualization import FieldPlotter fp = FieldPlotter(np.array([0.5, 0]), extent=3, npoints=200) + from pytools.obj_array import make_obj_array - targets = make_obj_array( - [fp.points[i] for i in range(knl.dim)]) + targets = make_obj_array([fp.points[i] for i in range(knl.dim)]) from boxtree import TreeBuilder - tb = TreeBuilder(ctx) + tb = TreeBuilder(actx.context) - tree, _ = tb(queue, sources, targets=targets, + tree, _ = tb(actx.queue, sources, targets=targets, max_particles_in_box=30, debug=True) from boxtree.traversal import FMMTraversalBuilder - tbuild = FMMTraversalBuilder(ctx) - trav, _ = tbuild(queue, tree, debug=True) + tbuild = FMMTraversalBuilder(actx.context) + trav, _ = tbuild(actx.queue, tree, debug=True) # {{{ plot tree if 0: - host_tree = tree.get(queue) - host_trav = trav.get(queue) + host_tree = tree.get(actx.queue) + host_trav = trav.get(actx.queue) if 0: - print("src_box", host_tree.find_box_nr_for_source(403)) - print("tgt_box", host_tree.find_box_nr_for_target(28)) - print(list(host_trav.target_or_target_parent_boxes).index(37)) - print(host_trav.get_box_list("sep_bigger", 22)) + logger.info("src_box: %s", host_tree.find_box_nr_for_source(403)) + logger.info("tgt_box: %s", host_tree.find_box_nr_for_target(28)) + logger.info("%s", + list(host_trav.target_or_target_parent_boxes).index(37)) + logger.info("%s", host_trav.get_box_list("sep_bigger", 22)) from boxtree.visualization import TreePlotter plotter = TreePlotter(host_tree) @@ -170,14 +174,11 @@ def _test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class, # }}} - from pyopencl.clrandom import PhiloxGenerator - rng = PhiloxGenerator(ctx, seed=44) - weights = rng.uniform(queue, nsources, dtype=np.float64) - + rng = np.random.default_rng(44) + weights = actx.from_numpy(rng.random(nsources, dtype=np.float64)) logger.info("computing direct (reference) result") from pytools.convergence import PConvergenceVerifier - pconv_verifier = PConvergenceVerifier() extra_kwargs = {} @@ -201,7 +202,6 @@ def _test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class, elif knl.dim == 2 and issubclass(local_expn_class, Y2DLocalExpansion): order_values = [10, 12] - from functools import partial for order in order_values: target_kernels = [knl] @@ -216,7 +216,7 @@ def _test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class, knl, local_expn_class)() tree_indep = SumpyTreeIndependentDataForWrangler( - ctx, + actx.context, partial(mpole_expn_class, knl), partial(local_expn_class, knl, m2l_translation=m2l_translation), target_kernels) @@ -238,60 +238,60 @@ def fmm_level_to_order(kernel, kernel_args, tree, lev): pot, = drive_fmm(wrangler, (weights,)) from sumpy import P2P - p2p = P2P(ctx, target_kernels, exclude_self=False) - evt, (ref_pot,) = p2p(queue, targets, sources, (weights,), + p2p = P2P(actx.context, target_kernels, exclude_self=False) + evt, (ref_pot,) = p2p(actx.queue, targets, sources, (weights,), **extra_kwargs) - pot = pot.get() - ref_pot = ref_pot.get() + pot = actx.to_numpy(pot) + ref_pot = actx.to_numpy(ref_pot) rel_err = la.norm(pot - ref_pot, np.inf) / la.norm(ref_pot, np.inf) logger.info("order %d -> relative l2 error: %g", order, rel_err) pconv_verifier.add_data_point(order, rel_err) - print(pconv_verifier) + logger.info("\n%s", pconv_verifier) pconv_verifier() +# }}} + + +# {{{ test_coeff_magnitude_rscale @pytest.mark.parametrize("knl", [LaplaceKernel(2), BiharmonicKernel(2)]) -def test_coeff_magnitude_rscale(ctx_factory, knl): +def test_coeff_magnitude_rscale(actx_factory, knl): """Checks that the rscale used keeps the coefficient magnitude difference small """ local_expn_class = LinearPDEConformingVolumeTaylorLocalExpansion mpole_expn_class = LinearPDEConformingVolumeTaylorMultipoleExpansion - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) + actx = actx_factory() nsources = 1000 ntargets = 300 dtype = np.float64 - from boxtree.tools import ( - make_normal_particle_array as p_normal) + from boxtree.tools import make_normal_particle_array as p_normal - sources = p_normal(queue, nsources, knl.dim, dtype, seed=15) + sources = p_normal(actx.queue, nsources, knl.dim, dtype, seed=15) offset = np.zeros(knl.dim) offset[0] = 0.1 - targets = ( - p_normal(queue, ntargets, knl.dim, dtype, seed=18) + offset) + targets = offset + p_normal(actx.queue, ntargets, knl.dim, dtype, seed=18) from boxtree import TreeBuilder - tb = TreeBuilder(ctx) + tb = TreeBuilder(actx.context) - tree, _ = tb(queue, sources, targets=targets, + tree, _ = tb(actx.queue, sources, targets=targets, max_particles_in_box=30, debug=True) from boxtree.traversal import FMMTraversalBuilder - tbuild = FMMTraversalBuilder(ctx) - trav, _ = tbuild(queue, tree, debug=True) + tbuild = FMMTraversalBuilder(actx.context) + trav, _ = tbuild(actx.queue, tree, debug=True) - from pyopencl.clrandom import PhiloxGenerator - rng = PhiloxGenerator(ctx, seed=44) - weights = rng.uniform(queue, nsources, dtype=np.float64) + rng = np.random.default_rng(31) + weights = actx.from_numpy(rng.random(nsources, dtype=np.float64)) extra_kwargs = {} dtype = np.float64 @@ -304,11 +304,10 @@ def test_coeff_magnitude_rscale(ctx_factory, knl): extra_kwargs["lam"] = 2 dtype = np.complex128 - from functools import partial target_kernels = [knl] tree_indep = SumpyTreeIndependentDataForWrangler( - ctx, + actx.context, partial(mpole_expn_class, knl), partial(local_expn_class, knl), target_kernels) @@ -330,21 +329,28 @@ def fmm_level_to_order(kernel, kernel_args, tree, lev): trav.from_sep_bigger_lists, (weights,)) - result = np.abs(wrangler.local_expansions_view(local_result, 5)[1][0]) + result = actx.to_numpy( + actx.np.abs(wrangler.local_expansions_view(local_result, 5)[1][0]) + ) + + result_ratio = np.max(result) / np.min(result) + assert result_ratio < 10**6, result_ratio + +# }}} - assert np.max(result) / np.min(result) < 10**6 +# {{{ test_unified_single_and_double -def test_unified_single_and_double(ctx_factory): +def test_unified_single_and_double(actx_factory, visualize=False): """ Test that running one FMM for single layer + double layer gives the same result as running one FMM for each and adding the results together at the end """ - logging.basicConfig(level=logging.INFO) + if visualize: + logging.basicConfig(level=logging.INFO) - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) + actx = actx_factory() knl = LaplaceKernel(2) local_expn_class = LinearPDEConformingVolumeTaylorLocalExpansion @@ -354,42 +360,36 @@ def test_unified_single_and_double(ctx_factory): ntargets = 300 dtype = np.float64 - from boxtree.tools import ( - make_normal_particle_array as p_normal) + from boxtree.tools import make_normal_particle_array as p_normal - sources = p_normal(queue, nsources, knl.dim, dtype, seed=15) + sources = p_normal(actx.queue, nsources, knl.dim, dtype, seed=15) offset = np.zeros(knl.dim) offset[0] = 0.1 - targets = ( - p_normal(queue, ntargets, knl.dim, dtype, seed=18) - + offset) - + targets = offset + p_normal(actx.queue, ntargets, knl.dim, dtype, seed=18) del offset from boxtree import TreeBuilder - tb = TreeBuilder(ctx) + tb = TreeBuilder(actx.context) - tree, _ = tb(queue, sources, targets=targets, + tree, _ = tb(actx.queue, sources, targets=targets, max_particles_in_box=30, debug=True) from boxtree.traversal import FMMTraversalBuilder - tbuild = FMMTraversalBuilder(ctx) - trav, _ = tbuild(queue, tree, debug=True) + tbuild = FMMTraversalBuilder(actx.context) + trav, _ = tbuild(actx.queue, tree, debug=True) - from pyopencl.clrandom import PhiloxGenerator - rng = PhiloxGenerator(ctx, seed=44) + rng = np.random.default_rng(44) weights = ( - rng.uniform(queue, nsources, dtype=np.float64), - rng.uniform(queue, nsources, dtype=np.float64), - ) + actx.from_numpy(rng.random(nsources, dtype=np.float64)), + actx.from_numpy(rng.random(nsources, dtype=np.float64)) + ) logger.info("computing direct (reference) result") dtype = np.float64 order = 3 - from functools import partial from sumpy.kernel import DirectionalSourceDerivative, AxisTargetDerivative deriv_knl = DirectionalSourceDerivative(knl, "dir_vec") @@ -407,7 +407,7 @@ def test_unified_single_and_double(ctx_factory): if deriv_knl in source_kernels: source_extra_kwargs["dir_vec"] = dir_vec tree_indep = SumpyTreeIndependentDataForWrangler( - ctx, + actx.context, partial(mpole_expn_class, knl), partial(local_expn_class, knl), target_kernels=target_kernels, source_kernels=source_kernels, @@ -419,7 +419,7 @@ def test_unified_single_and_double(ctx_factory): from boxtree.fmm import drive_fmm pot = drive_fmm(wrangler, weights) - results.append(np.array([pot[0].get(), pot[1].get()])) + results.append(np.array([actx.to_numpy(pot[0]), actx.to_numpy(pot[1])])) ref_pot = results[0] + results[1] pot = results[2] @@ -427,47 +427,50 @@ def test_unified_single_and_double(ctx_factory): assert rel_err < 1e-12 +# }}} + + +# {{{ test_sumpy_fmm_timing_data_collection @pytest.mark.parametrize("use_fft", [True, False]) -def test_sumpy_fmm_timing_data_collection(ctx_factory, use_fft): - logging.basicConfig(level=logging.INFO) +def test_sumpy_fmm_timing_data_collection(ctx_factory, use_fft, visualize=False): + if visualize: + logging.basicConfig(level=logging.INFO) + + import pyopencl as cl + from sumpy.array_context import PyOpenCLArrayContext ctx = ctx_factory() - queue = cl.CommandQueue( - ctx, - properties=cl.command_queue_properties.PROFILING_ENABLE) + queue = cl.CommandQueue(ctx, + properties=cl.command_queue_properties.PROFILING_ENABLE) + actx = PyOpenCLArrayContext(queue, force_device_scalars=True) nsources = 500 dtype = np.float64 - from boxtree.tools import ( - make_normal_particle_array as p_normal) + from boxtree.tools import make_normal_particle_array as p_normal knl = LaplaceKernel(2) local_expn_class = VolumeTaylorLocalExpansion mpole_expn_class = VolumeTaylorMultipoleExpansion order = 1 - sources = p_normal(queue, nsources, knl.dim, dtype, seed=15) + sources = p_normal(actx.queue, nsources, knl.dim, dtype, seed=15) from boxtree import TreeBuilder - tb = TreeBuilder(ctx) + tb = TreeBuilder(actx.context) - tree, _ = tb(queue, sources, - max_particles_in_box=30, debug=True) + tree, _ = tb(actx.queue, sources, max_particles_in_box=30, debug=True) from boxtree.traversal import FMMTraversalBuilder - tbuild = FMMTraversalBuilder(ctx) - trav, _ = tbuild(queue, tree, debug=True) + tbuild = FMMTraversalBuilder(actx.context) + trav, _ = tbuild(actx.queue, tree, debug=True) - from pyopencl.clrandom import PhiloxGenerator - rng = PhiloxGenerator(ctx) - weights = rng.uniform(queue, nsources, dtype=np.float64) + rng = np.random.default_rng(44) + weights = actx.from_numpy(rng.random(nsources, dtype=np.float64)) target_kernels = [knl] - from functools import partial - if use_fft: from sumpy.expansion.m2l import FFTM2LTranslationClassFactory m2l_translation_factory = FFTM2LTranslationClassFactory() @@ -479,7 +482,7 @@ def test_sumpy_fmm_timing_data_collection(ctx_factory, use_fft): knl, local_expn_class)() tree_indep = SumpyTreeIndependentDataForWrangler( - ctx, + actx.context, partial(mpole_expn_class, knl), partial(local_expn_class, knl, m2l_translation=m2l_translation), target_kernels) @@ -490,52 +493,48 @@ def test_sumpy_fmm_timing_data_collection(ctx_factory, use_fft): timing_data = {} pot, = drive_fmm(wrangler, (weights,), timing_data=timing_data) - print(timing_data) + logger.info("timing_data:\n%s", timing_data) + assert timing_data -def test_sumpy_fmm_exclude_self(ctx_factory): - logging.basicConfig(level=logging.INFO) +def test_sumpy_fmm_exclude_self(actx_factory, visualize=False): + if visualize: + logging.basicConfig(level=logging.INFO) - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) + actx = actx_factory() nsources = 500 dtype = np.float64 - from boxtree.tools import ( - make_normal_particle_array as p_normal) + from boxtree.tools import make_normal_particle_array as p_normal knl = LaplaceKernel(2) local_expn_class = VolumeTaylorLocalExpansion mpole_expn_class = VolumeTaylorMultipoleExpansion order = 10 - sources = p_normal(queue, nsources, knl.dim, dtype, seed=15) + sources = p_normal(actx.queue, nsources, knl.dim, dtype, seed=15) from boxtree import TreeBuilder - tb = TreeBuilder(ctx) + tb = TreeBuilder(actx.context) - tree, _ = tb(queue, sources, - max_particles_in_box=30, debug=True) + tree, _ = tb(actx.queue, sources, max_particles_in_box=30, debug=True) from boxtree.traversal import FMMTraversalBuilder - tbuild = FMMTraversalBuilder(ctx) - trav, _ = tbuild(queue, tree, debug=True) + tbuild = FMMTraversalBuilder(actx.context) + trav, _ = tbuild(actx.queue, tree, debug=True) - from pyopencl.clrandom import PhiloxGenerator - rng = PhiloxGenerator(ctx) - weights = rng.uniform(queue, nsources, dtype=np.float64) + rng = np.random.default_rng(44) + weights = actx.from_numpy(rng.random(nsources, dtype=np.float64)) target_to_source = np.arange(tree.ntargets, dtype=np.int32) self_extra_kwargs = {"target_to_source": target_to_source} target_kernels = [knl] - from functools import partial - tree_indep = SumpyTreeIndependentDataForWrangler( - ctx, + actx.context, partial(mpole_expn_class, knl), partial(local_expn_class, knl), target_kernels, @@ -550,65 +549,65 @@ def test_sumpy_fmm_exclude_self(ctx_factory): pot, = drive_fmm(wrangler, (weights,)) from sumpy import P2P - p2p = P2P(ctx, target_kernels, exclude_self=True) - evt, (ref_pot,) = p2p(queue, sources, sources, (weights,), + p2p = P2P(actx.context, target_kernels, exclude_self=True) + evt, (ref_pot,) = p2p(actx.queue, sources, sources, (weights,), **self_extra_kwargs) - pot = pot.get() - ref_pot = ref_pot.get() + pot = actx.to_numpy(pot) + ref_pot = actx.to_numpy(ref_pot) rel_err = la.norm(pot - ref_pot) / la.norm(ref_pot) logger.info("order %d -> relative l2 error: %g", order, rel_err) assert np.isclose(rel_err, 0, atol=1e-7) +# }}} -def test_sumpy_axis_source_derivative(ctx_factory): - logging.basicConfig(level=logging.INFO) - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) +# {{{ test_sumpy_axis_source_derivative + +def test_sumpy_axis_source_derivative(actx_factory, visualize=False): + if visualize: + logging.basicConfig(level=logging.INFO) + + actx = actx_factory() nsources = 500 dtype = np.float64 - from boxtree.tools import ( - make_normal_particle_array as p_normal) + from boxtree.tools import make_normal_particle_array as p_normal knl = LaplaceKernel(2) local_expn_class = VolumeTaylorLocalExpansion mpole_expn_class = VolumeTaylorMultipoleExpansion order = 10 - sources = p_normal(queue, nsources, knl.dim, dtype, seed=15) + sources = p_normal(actx.queue, nsources, knl.dim, dtype, seed=15) from boxtree import TreeBuilder - tb = TreeBuilder(ctx) + tb = TreeBuilder(actx.context) - tree, _ = tb(queue, sources, + tree, _ = tb(actx.queue, sources, max_particles_in_box=30, debug=True) from boxtree.traversal import FMMTraversalBuilder - tbuild = FMMTraversalBuilder(ctx) - trav, _ = tbuild(queue, tree, debug=True) + tbuild = FMMTraversalBuilder(actx.context) + trav, _ = tbuild(actx.queue, tree, debug=True) - from pyopencl.clrandom import PhiloxGenerator - rng = PhiloxGenerator(ctx, seed=12) - weights = rng.uniform(queue, nsources, dtype=np.float64) + rng = np.random.default_rng(12) + weights = actx.from_numpy(rng.random(nsources, dtype=np.float64)) target_to_source = np.arange(tree.ntargets, dtype=np.int32) self_extra_kwargs = {"target_to_source": target_to_source} - from functools import partial - from sumpy.kernel import AxisTargetDerivative, AxisSourceDerivative pots = [] - for tgt_knl, src_knl in [(AxisTargetDerivative(0, knl), knl), + for tgt_knl, src_knl in [ + (AxisTargetDerivative(0, knl), knl), (knl, AxisSourceDerivative(0, knl))]: - tree_indep = SumpyTreeIndependentDataForWrangler( - ctx, + actx.context, partial(mpole_expn_class, knl), partial(local_expn_class, knl), target_kernels=[tgt_knl], @@ -622,53 +621,53 @@ def test_sumpy_axis_source_derivative(ctx_factory): from boxtree.fmm import drive_fmm pot, = drive_fmm(wrangler, (weights,)) - pots.append(pot.get()) + pots.append(actx.to_numpy(pot)) rel_err = la.norm(pots[0] + pots[1]) / la.norm(pots[0]) logger.info("order %d -> relative l2 error: %g", order, rel_err) assert np.isclose(rel_err, 0, atol=1e-5) +# }}} + + +# {{{ test_sumpy_target_point_multiplier @pytest.mark.parametrize("deriv_axes", [(), (0,), (1,)]) -def test_sumpy_target_point_multiplier(ctx_factory, deriv_axes): - logging.basicConfig(level=logging.INFO) +def test_sumpy_target_point_multiplier(actx_factory, deriv_axes, visualize=False): + if visualize: + logging.basicConfig(level=logging.INFO) - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) + actx = actx_factory() nsources = 500 dtype = np.float64 - from boxtree.tools import ( - make_normal_particle_array as p_normal) + from boxtree.tools import make_normal_particle_array as p_normal knl = LaplaceKernel(2) local_expn_class = VolumeTaylorLocalExpansion mpole_expn_class = VolumeTaylorMultipoleExpansion order = 5 - sources = p_normal(queue, nsources, knl.dim, dtype, seed=15) + sources = p_normal(actx.queue, nsources, knl.dim, dtype, seed=15) from boxtree import TreeBuilder - tb = TreeBuilder(ctx) + tb = TreeBuilder(actx.context) - tree, _ = tb(queue, sources, + tree, _ = tb(actx.queue, sources, max_particles_in_box=30, debug=True) from boxtree.traversal import FMMTraversalBuilder - tbuild = FMMTraversalBuilder(ctx) - trav, _ = tbuild(queue, tree, debug=True) + tbuild = FMMTraversalBuilder(actx.context) + trav, _ = tbuild(actx.queue, tree, debug=True) - from pyopencl.clrandom import PhiloxGenerator - rng = PhiloxGenerator(ctx, seed=12) - weights = rng.uniform(queue, nsources, dtype=np.float64) + rng = np.random.default_rng(12) + weights = actx.from_numpy(rng.random(nsources, dtype=np.float64)) target_to_source = np.arange(tree.ntargets, dtype=np.int32) self_extra_kwargs = {"target_to_source": target_to_source} - from functools import partial - from sumpy.kernel import TargetPointMultiplier, AxisTargetDerivative tgt_knls = [TargetPointMultiplier(0, knl), knl, knl] @@ -677,7 +676,7 @@ def test_sumpy_target_point_multiplier(ctx_factory, deriv_axes): tgt_knls[1] = AxisTargetDerivative(axis, tgt_knls[1]) tree_indep = SumpyTreeIndependentDataForWrangler( - ctx, + actx.context, partial(mpole_expn_class, knl), partial(local_expn_class, knl), target_kernels=tgt_knls, @@ -691,27 +690,29 @@ def test_sumpy_target_point_multiplier(ctx_factory, deriv_axes): from boxtree.fmm import drive_fmm pot0, pot1, pot2 = drive_fmm(wrangler, (weights,)) - pot0, pot1, pot2 = pot0.get(), pot1.get(), pot2.get() + pot0, pot1, pot2 = actx.to_numpy(pot0), actx.to_numpy(pot1), actx.to_numpy(pot2) if deriv_axes == (0,): - ref_pot = pot1 * sources[0].get() + pot2 + ref_pot = pot1 * actx.to_numpy(sources[0]) + pot2 else: - ref_pot = pot1 * sources[0].get() + ref_pot = pot1 * actx.to_numpy(sources[0]) rel_err = la.norm(pot0 - ref_pot) / la.norm(ref_pot) logger.info("order %d -> relative l2 error: %g", order, rel_err) assert np.isclose(rel_err, 0, atol=1e-5) +# }}} + # You can test individual routines by typing -# $ python test_fmm.py 'test_sumpy_fmm(cl.create_some_context, LaplaceKernel(2), -# VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion, False, False)' +# $ python test_fmm.py 'test_sumpy_fmm(_acf, LaplaceKernel(2), +# VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion, False, False, +# visualize=True)' if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from pytest import main - main([__file__]) + pytest.main([__file__]) # vim: fdm=marker From c198ef075700b937a60a42a5d1a167098e6e2a34 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 20 Aug 2022 20:23:50 +0300 Subject: [PATCH 3/8] port test_kernels to arraycontext --- sumpy/symbolic.py | 6 + test/test_kernels.py | 287 ++++++++++++++++++++++--------------------- 2 files changed, 150 insertions(+), 143 deletions(-) diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py index 4e1234ae2..ab4c25dcf 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -137,9 +137,15 @@ def UnevaluatedExpr(x): # noqa if USE_SYMENGINE: + def doit(expr): + return expr + def unevaluated_pow(a, b): return sym.Pow(a, b) else: + def doit(expr): + return expr.doit() + def unevaluated_pow(a, b): return sym.Pow(a, b, evaluate=False) diff --git a/test/test_kernels.py b/test/test_kernels.py index f2d89b860..bf03f00d4 100644 --- a/test/test_kernels.py +++ b/test/test_kernels.py @@ -20,55 +20,62 @@ THE SOFTWARE. """ +import pytest +import sys + import numpy as np import numpy.linalg as la -import sys -import pytest -import pyopencl as cl -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) +from pytools.convergence import PConvergenceVerifier +from arraycontext import pytest_generate_tests_for_array_contexts +from sumpy.array_context import ( # noqa: F401 + PytestPyOpenCLArrayContextFactory, _acf) +import sumpy.symbolic as sym from sumpy.expansion.multipole import ( - VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion, - VolumeTaylorMultipoleExpansionBase, - LinearPDEConformingVolumeTaylorMultipoleExpansion) + VolumeTaylorMultipoleExpansion, + H2DMultipoleExpansion, + VolumeTaylorMultipoleExpansionBase, + LinearPDEConformingVolumeTaylorMultipoleExpansion) from sumpy.expansion.local import ( - VolumeTaylorLocalExpansion, H2DLocalExpansion, - LinearPDEConformingVolumeTaylorLocalExpansion) + VolumeTaylorLocalExpansion, + H2DLocalExpansion, + LinearPDEConformingVolumeTaylorLocalExpansion) from sumpy.expansion.m2l import NonFFTM2LTranslationClassFactory -from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative, - DirectionalSourceDerivative, BiharmonicKernel, StokesletKernel) -import sumpy.symbolic as sym -from pytools.convergence import PConvergenceVerifier +from sumpy.kernel import ( + LaplaceKernel, + HelmholtzKernel, + BiharmonicKernel, + StokesletKernel, + AxisTargetDerivative, + DirectionalSourceDerivative) import logging logger = logging.getLogger(__name__) -try: - import faulthandler -except ImportError: - pass -else: - faulthandler.enable() +pytest_generate_tests = pytest_generate_tests_for_array_contexts([ + PytestPyOpenCLArrayContextFactory, + ]) + +# {{{ test_p2p @pytest.mark.parametrize("exclude_self", (True, False)) -def test_p2p(ctx_factory, exclude_self): - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) +def test_p2p(actx_factory, exclude_self): + actx = actx_factory() dimensions = 3 n = 5000 from sumpy.p2p import P2P lknl = LaplaceKernel(dimensions) - knl = P2P(ctx, + knl = P2P(actx.context, [lknl, AxisTargetDerivative(0, lknl)], exclude_self=exclude_self) - targets = np.random.rand(dimensions, n) - sources = targets if exclude_self else np.random.rand(dimensions, n) + rng = np.random.default_rng(42) + targets = rng.random(size=(dimensions, n)) + sources = targets if exclude_self else rng.random(size=(dimensions, n)) strengths = np.ones(n, dtype=np.float64) @@ -78,7 +85,7 @@ def test_p2p(ctx_factory, exclude_self): extra_kwargs["target_to_source"] = np.arange(n, dtype=np.int32) evt, (potential, x_derivative) = knl( - queue, targets, sources, [strengths], + actx.queue, targets, sources, [strengths], out_host=True, **extra_kwargs) potential_ref = np.empty_like(potential) @@ -99,24 +106,22 @@ def test_p2p(ctx_factory, exclude_self): potential_ref *= 1/(4*np.pi) rel_err = la.norm(potential - potential_ref)/la.norm(potential_ref) - print(rel_err) + logger.info("error: %.12e", rel_err) + assert rel_err < 1e-3 +# }}} + + +# {{{ test_p2e_multiple @pytest.mark.parametrize(("base_knl", "expn_class"), [ (LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion), (LaplaceKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion), ]) -def test_p2e_multiple(ctx_factory, base_knl, expn_class): - - from sympy.core.cache import clear_cache - clear_cache() - +def test_p2e_multiple(actx_factory, base_knl, expn_class): order = 4 - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) - - np.random.seed(17) + actx = actx_factory() nsources = 100 @@ -139,9 +144,11 @@ def test_p2e_multiple(ctx_factory, base_knl, expn_class): from sumpy import P2EFromSingleBox + rng = np.random.default_rng(14) center = np.array([2, 1, 0][:knl.dim], np.float64) - sources = (0.7*(-0.5+np.random.rand(knl.dim, nsources).astype(np.float64)) - + center[:, np.newaxis]) + sources = ( + 0.7 * (-0.5 + rng.random(size=(knl.dim, nsources), dtype=np.float64)) + + center[:, np.newaxis]) strengths = [ np.ones(nsources, dtype=np.float64) * (1/nsources), @@ -167,8 +174,11 @@ def test_p2e_multiple(ctx_factory, base_knl, expn_class): rscale = 0.5 # pick something non-1 # apply p2e at the same time - p2e = P2EFromSingleBox(ctx, expn, kernels=source_kernels, strength_usage=[0, 1]) - evt, (mpoles,) = p2e(queue, + p2e = P2EFromSingleBox(actx.context, expn, + kernels=source_kernels, + strength_usage=[0, 1]) + + evt, (mpoles,) = p2e(actx.queue, source_boxes=source_boxes, box_source_starts=box_source_starts, box_source_counts_nonchild=box_source_counts_nonchild, @@ -179,7 +189,6 @@ def test_p2e_multiple(ctx_factory, base_knl, expn_class): tgt_base_ibox=0, rscale=rscale, - #flags="print_hl_cl", out_host=True, dir_vec=dir_vec, **extra_kwargs) @@ -192,9 +201,11 @@ def test_p2e_multiple(ctx_factory, base_knl, expn_class): extra_source_kwargs = extra_kwargs.copy() if isinstance(source_kernel, DirectionalSourceDerivative): extra_source_kwargs["dir_vec"] = dir_vec - p2e = P2EFromSingleBox(ctx, expn, + + p2e = P2EFromSingleBox(actx.context, expn, kernels=[source_kernel], strength_usage=[i]) - evt, (mpoles,) = p2e(queue, + + evt, (mpoles,) = p2e(actx.queue, source_boxes=source_boxes, box_source_starts=box_source_starts, box_source_counts_nonchild=box_source_counts_nonchild, @@ -205,13 +216,16 @@ def test_p2e_multiple(ctx_factory, base_knl, expn_class): tgt_base_ibox=0, rscale=rscale, - #flags="print_hl_cl", out_host=True, **extra_source_kwargs) expected_result += mpoles norm = la.norm(actual_result - expected_result)/la.norm(expected_result) assert norm < 1e-12 +# }}} + + +# {{{ test_p2e2p @pytest.mark.parametrize("order", [4]) @pytest.mark.parametrize(("base_knl", "expn_class"), [ @@ -245,17 +259,8 @@ def test_p2e_multiple(ctx_factory, base_knl, expn_class): False, True ]) -# Sample: test_p2e2p(cl._csc, LaplaceKernel(2), VolumeTaylorLocalExpansion, 4, False) -def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative): - #logging.basicConfig(level=logging.INFO) - - from sympy.core.cache import clear_cache - clear_cache() - - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) - - np.random.seed(17) +def test_p2e2p(actx_factory, base_knl, expn_class, order, with_source_derivative): + actx = actx_factory() res = 100 nsources = 100 @@ -281,9 +286,9 @@ def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative) expn = expn_class(knl, order=order) from sumpy import P2EFromSingleBox, E2PFromSingleBox, P2P - p2e = P2EFromSingleBox(ctx, expn, kernels=[knl]) - e2p = E2PFromSingleBox(ctx, expn, kernels=target_kernels) - p2p = P2P(ctx, target_kernels, exclude_self=False) + p2e = P2EFromSingleBox(actx.context, expn, kernels=[knl]) + e2p = E2PFromSingleBox(actx.context, expn, kernels=target_kernels) + p2p = P2P(actx.context, target_kernels, exclude_self=False) from pytools.convergence import EOCRecorder eoc_rec_pot = EOCRecorder() @@ -295,11 +300,13 @@ def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative) else: h_values = [1/2, 1/3, 1/5] + rng = np.random.default_rng(19) center = np.array([2, 1, 0][:knl.dim], np.float64) - sources = (0.7*(-0.5+np.random.rand(knl.dim, nsources).astype(np.float64)) - + center[:, np.newaxis]) + sources = ( + 0.7 * (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64)) + + center[:, np.newaxis]) - strengths = np.ones(nsources, dtype=np.float64) * (1/nsources) + strengths = np.ones(nsources, dtype=np.float64) / nsources source_boxes = np.array([0], dtype=np.int32) box_source_starts = np.array([0], dtype=np.int32) @@ -331,7 +338,7 @@ def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative) # {{{ apply p2e - evt, (mpoles,) = p2e(queue, + evt, (mpoles,) = p2e(actx.queue, source_boxes=source_boxes, box_source_starts=box_source_starts, box_source_counts_nonchild=box_source_counts_nonchild, @@ -342,7 +349,6 @@ def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative) tgt_base_ibox=0, rscale=rscale, - #flags="print_hl_cl", out_host=True, **extra_source_kwargs) # }}} @@ -355,7 +361,7 @@ def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative) box_target_counts_nonchild = np.array([ntargets], dtype=np.int32) evt, (pot, grad_x, ) = e2p( - queue, + actx.queue, src_expansions=mpoles, src_base_ibox=0, target_boxes=source_boxes, @@ -365,7 +371,6 @@ def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative) targets=targets, rscale=rscale, - #flags="print_hl_cl", out_host=True, **extra_kwargs) # }}} @@ -373,7 +378,7 @@ def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative) # {{{ compute (direct) reference solution evt, (pot_direct, grad_x_direct, ) = p2p( - queue, + actx.queue, targets, sources, (strengths,), out_host=True, **extra_source_kwargs) @@ -410,11 +415,11 @@ def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative) eoc_rec_pot.add_data_point(h, err_pot) eoc_rec_grad_x.add_data_point(h, err_grad_x) - print(expn_class, knl, order) - print("POTENTIAL:") - print(eoc_rec_pot) - print("X TARGET DERIVATIVE:") - print(eoc_rec_grad_x) + logger.info("expn_cls %s knl %s order %d", expn_class, knl, order) + logger.info("POTENTIAL:") + logger.info("%s", eoc_rec_pot) + logger.info("X TARGET DERIVATIVE:") + logger.info("%s", eoc_rec_grad_x) tgt_order = order + 1 if issubclass(expn_class, LocalExpansionBase): @@ -451,6 +456,10 @@ def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative) assert eoc_rec_pot.order_estimate() > tgt_order - slack assert eoc_rec_grad_x.order_estimate() > tgt_order_grad - grad_slack +# }}} + + +# {{{ test_translations @pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class", [ (LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), @@ -471,16 +480,12 @@ def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative) (StokesletKernel(2, 0, 0), LinearPDEConformingVolumeTaylorLocalExpansion, LinearPDEConformingVolumeTaylorMultipoleExpansion), ]) -def test_translations(ctx_factory, knl, local_expn_class, mpole_expn_class): - logging.basicConfig(level=logging.INFO) - - from sympy.core.cache import clear_cache - clear_cache() +def test_translations(actx_factory, knl, local_expn_class, mpole_expn_class, + visualize=False): + if visualize: + logging.basicConfig(level=logging.INFO) - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) - - np.random.seed(17) + actx = actx_factory() res = 20 nsources = 15 @@ -494,9 +499,11 @@ def test_translations(ctx_factory, knl, local_expn_class, mpole_expn_class): extra_kwargs["mu"] = 0.05 # 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 = (0.7*(-0.5+np.random.rand(knl.dim, nsources).astype(np.float64)) - + origin[:, np.newaxis]) + sources = ( + 0.7 * (-0.5 + rng.random((knl.dim, nsources), dtype=np.float64)) + + origin[:, np.newaxis]) strengths = np.ones(nsources, dtype=np.float64) * (1/nsources) pconv_verifier_p2m2p = PConvergenceVerifier() @@ -543,7 +550,7 @@ def eval_at(e2p, source_box_nr, rscale): e2p_box_target_counts_nonchild[source_box_nr] = ntargets evt, (pot,) = e2p( - queue, + actx.queue, src_expansions=mpoles, src_base_ibox=0, @@ -569,13 +576,13 @@ def eval_at(e2p, source_box_nr, rscale): l_expn = local_expn_class(knl, order=order, m2l_translation=m2l_translation) from sumpy import P2EFromSingleBox, E2PFromSingleBox, P2P, E2EFromCSR - p2m = P2EFromSingleBox(ctx, m_expn) - m2m = E2EFromCSR(ctx, m_expn, m_expn) - m2p = E2PFromSingleBox(ctx, m_expn, target_kernels) - m2l = E2EFromCSR(ctx, m_expn, l_expn) - l2l = E2EFromCSR(ctx, l_expn, l_expn) - l2p = E2PFromSingleBox(ctx, l_expn, target_kernels) - p2p = P2P(ctx, target_kernels, exclude_self=False) + 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) fp = FieldPlotter(centers[:, -1], extent=0.3, npoints=res) targets = fp.points @@ -583,7 +590,7 @@ def eval_at(e2p, source_box_nr, rscale): # {{{ compute (direct) reference solution evt, (pot_direct,) = p2p( - queue, + actx.queue, targets, sources, (strengths,), out_host=True, **extra_kwargs) @@ -603,7 +610,7 @@ def eval_at(e2p, source_box_nr, rscale): p2m_box_source_counts_nonchild = np.array([nsources, 0, 0, 0], dtype=np.int32) - evt, (mpoles,) = p2m(queue, + 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, @@ -615,7 +622,6 @@ def eval_at(e2p, source_box_nr, rscale): tgt_base_ibox=0, - #flags="print_hl_wrapper", out_host=True, **extra_kwargs) # }}} @@ -635,7 +641,7 @@ def eval_at(e2p, source_box_nr, rscale): 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, + evt, (mpoles,) = m2m(actx.queue, src_expansions=mpoles, src_base_ibox=0, tgt_base_ibox=0, @@ -650,7 +656,6 @@ def eval_at(e2p, source_box_nr, rscale): src_rscale=m1_rscale, tgt_rscale=m2_rscale, - #flags="print_hl_cl", out_host=True, **extra_kwargs) # }}} @@ -668,7 +673,7 @@ def eval_at(e2p, source_box_nr, rscale): m2l_src_box_starts = np.array([0, 1], dtype=np.int32) m2l_src_box_lists = np.array([1], dtype=np.int32) - evt, (mpoles,) = m2l(queue, + evt, (mpoles,) = m2l(actx.queue, src_expansions=mpoles, src_base_ibox=0, tgt_base_ibox=0, @@ -682,7 +687,6 @@ def eval_at(e2p, source_box_nr, rscale): src_rscale=m2_rscale, tgt_rscale=l1_rscale, - #flags="print_hl_cl", out_host=True, **extra_kwargs) # }}} @@ -700,7 +704,7 @@ def eval_at(e2p, source_box_nr, rscale): l2l_src_box_starts = np.array([0, 1], dtype=np.int32) l2l_src_box_lists = np.array([2], dtype=np.int32) - evt, (mpoles,) = l2l(queue, + evt, (mpoles,) = l2l(actx.queue, src_expansions=mpoles, src_base_ibox=0, tgt_base_ibox=0, @@ -714,7 +718,6 @@ def eval_at(e2p, source_box_nr, rscale): src_rscale=l1_rscale, tgt_rscale=l2_rscale, - #flags="print_hl_wrapper", out_host=True, **extra_kwargs) # }}} @@ -732,13 +735,18 @@ def eval_at(e2p, source_box_nr, rscale): ("p2m2m2l2p", pconv_verifier_p2m2m2l2p), ("full", pconv_verifier_full), ]: - print(30*"-") - print(name) - print(30*"-") - print(verifier) - print(30*"-") + logger.info(30*"-") + logger.info("name: %s", name) + logger.info(30*"-") + logger.info("result: %s", verifier) + logger.info(30*"-") + verifier() +# }}} + + +# {{{ test_m2m_and_l2l_exprs_simpler @pytest.mark.parametrize("order", [4]) @pytest.mark.parametrize(("base_knl", "local_expn_class", "mpole_expn_class"), [ @@ -750,12 +758,6 @@ def eval_at(e2p, source_box_nr, rscale): ]) def test_m2m_and_l2l_exprs_simpler(base_knl, local_expn_class, mpole_expn_class, order, with_source_derivative): - - from sympy.core.cache import clear_cache - clear_cache() - - np.random.seed(17) - extra_kwargs = {} if isinstance(base_knl, HelmholtzKernel): if base_knl.allow_evanescent: @@ -773,9 +775,8 @@ def test_m2m_and_l2l_exprs_simpler(base_knl, local_expn_class, mpole_expn_class, mpole_expn = mpole_expn_class(knl, order=order) local_expn = local_expn_class(knl, order=order) - from sumpy.symbolic import make_sym_vector, Symbol, USE_SYMENGINE - dvec = make_sym_vector("d", knl.dim) - src_coeff_exprs = [Symbol(f"src_coeff{i}") for i in range(len(mpole_expn))] + dvec = sym.make_sym_vector("d", knl.dim) + src_coeff_exprs = [sym.Symbol(f"src_coeff{i}") for i in range(len(mpole_expn))] src_rscale = 3 tgt_rscale = 2 @@ -785,30 +786,23 @@ def test_m2m_and_l2l_exprs_simpler(base_knl, local_expn_class, mpole_expn_class, slower_m2m = mpole_expn.translate_from(mpole_expn, src_coeff_exprs, src_rscale, dvec, tgt_rscale, _fast_version=False) - def _check_equal(expr1, expr2): - if USE_SYMENGINE: - return float((expr1 - expr2).expand()) == 0.0 - else: - # with sympy we are using UnevaluatedExpr and expand doesn't expand it - # Running doit replaces UnevaluatedExpr with evaluated exprs - return float((expr1 - expr2).doit().expand()) == 0.0 - for expr1, expr2 in zip(faster_m2m, slower_m2m): - assert _check_equal(expr1, expr2) + assert float(sym.doit(expr1 - expr2).expand()) == 0.0 faster_l2l = local_expn.translate_from(local_expn, src_coeff_exprs, src_rscale, dvec, tgt_rscale) slower_l2l = local_expn.translate_from(local_expn, src_coeff_exprs, src_rscale, dvec, tgt_rscale, _fast_version=False) for expr1, expr2 in zip(faster_l2l, slower_l2l): - assert _check_equal(expr1, expr2) + assert float(sym.doit(expr1 - expr2).expand()) == 0.0 + +# }}} # {{{ test toeplitz def _m2l_translate_simple(tgt_expansion, src_expansion, src_coeff_exprs, src_rscale, dvec, tgt_rscale): - if not tgt_expansion.use_rscale: src_rscale = 1 tgt_rscale = 1 @@ -839,12 +833,14 @@ def _m2l_translate_simple(tgt_expansion, src_expansion, src_coeff_exprs, src_rsc local_result.append( coeff * kernel_deriv * tgt_rscale**sum(deriv)) result.append(sym.Add(*local_result)) + return result def test_m2l_toeplitz(): dim = 3 knl = LaplaceKernel(dim) + local_expn_class = LinearPDEConformingVolumeTaylorLocalExpansion mpole_expn_class = LinearPDEConformingVolumeTaylorMultipoleExpansion m2l_factory = NonFFTM2LTranslationClassFactory() @@ -854,19 +850,24 @@ def test_m2l_toeplitz(): mpole_expn = mpole_expn_class(knl, order=5) dvec = sym.make_sym_vector("d", dim) - src_coeff_exprs = list(1 + np.random.randn(len(mpole_expn))) + + rng = np.random.default_rng(44) + src_coeff_exprs = list(1 + rng.standard_normal(len(mpole_expn))) src_rscale = 2.0 tgt_rscale = 1.0 - expected_output = _m2l_translate_simple(local_expn, mpole_expn, src_coeff_exprs, - src_rscale, dvec, tgt_rscale) - actual_output = local_expn.translate_from(mpole_expn, src_coeff_exprs, - src_rscale, dvec, tgt_rscale, sac=None) + expected_output = _m2l_translate_simple( + local_expn, mpole_expn, src_coeff_exprs, + src_rscale, dvec, tgt_rscale) + actual_output = local_expn.translate_from( + mpole_expn, src_coeff_exprs, + src_rscale, dvec, tgt_rscale, sac=None) - replace_dict = {d: np.random.rand(1)[0] for d in dvec} + replace_dict = {d: rng.random() for d in dvec} for sym_a, sym_b in zip(expected_output, actual_output): num_a = sym_a.xreplace(replace_dict) num_b = sym_b.xreplace(replace_dict) + assert abs(num_a - num_b)/abs(num_a) < 1e-10 # }}} @@ -876,10 +877,11 @@ def test_m2l_toeplitz(): @pytest.mark.parametrize("dim", [2, 3]) @pytest.mark.parametrize("order", [2, 4, 6]) -def test_m2m_compressed_error_helmholtz(ctx_factory, dim, order): - import sumpy.toys as t +def test_m2m_compressed_error_helmholtz(actx_factory, dim, order): + from sumpy import toys + + actx = actx_factory() - ctx = ctx_factory() knl = HelmholtzKernel(dim) extra_kernel_kwargs = {"k": 5} @@ -914,24 +916,24 @@ def test_m2m_compressed_error_helmholtz(ctx_factory, dim, order): m2m_vals = [0, 0] for i, (mpole_expn_class, local_expn_class) in \ enumerate(zip(mpole_expn_classes, local_expn_classes)): - tctx = t.ToyContext( - ctx, + tctx = toys.ToyContext( + actx.context, knl, extra_kernel_kwargs=extra_kernel_kwargs, local_expn_class=local_expn_class, mpole_expn_class=mpole_expn_class, ) - pt_src = t.PointSources( + pt_src = toys.PointSources( tctx, sources, np.ones(sources.shape[-1]) ) - mexp = t.multipole_expand(pt_src, + mexp = toys.multipole_expand(pt_src, center=mpole_center.reshape(dim), order=order, rscale=h) - mexp2 = t.multipole_expand(mexp, + mexp2 = toys.multipole_expand(mexp, center=second_center.reshape(dim), order=order, rscale=h) @@ -941,20 +943,19 @@ def test_m2m_compressed_error_helmholtz(ctx_factory, dim, order): / np.linalg.norm(m2m_vals[1]) eoc_rec.add_data_point(furthest_source, err) - print(eoc_rec) + logger.info("\n%s", eoc_rec) assert eoc_rec.order_estimate() >= order + 1 # }}} # You can test individual routines by typing -# $ python test_kernels.py 'test_p2p(cl.create_some_context)' +# $ python test_kernels.py 'test_p2p(_acf, True)' if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from pytest import main - main([__file__]) + pytest.main([__file__]) # vim: fdm=marker From 0a97945a2c9ba5640d44ebe1307b5098a3052422 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 20 Aug 2022 20:36:32 +0300 Subject: [PATCH 4/8] port test_matrixgen to arraycontext --- test/test_matrixgen.py | 151 ++++++++++++++++++++++------------------- 1 file changed, 81 insertions(+), 70 deletions(-) diff --git a/test/test_matrixgen.py b/test/test_matrixgen.py index ed90ee1e5..4baf30312 100644 --- a/test/test_matrixgen.py +++ b/test/test_matrixgen.py @@ -20,28 +20,25 @@ THE SOFTWARE. """ +import pytest import sys + import numpy as np import numpy.linalg as la -import pyopencl as cl -import pyopencl.array # noqa - -from sumpy.tools import vector_to_device - -import pytest -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) - +from arraycontext import pytest_generate_tests_for_array_contexts +from sumpy.array_context import ( # noqa: F401 + PytestPyOpenCLArrayContextFactory, _acf) import logging logger = logging.getLogger(__name__) -import faulthandler -faulthandler.enable() +pytest_generate_tests = pytest_generate_tests_for_array_contexts([ + PytestPyOpenCLArrayContextFactory, + ]) -def _build_geometry(queue, ntargets, nsources, mode, target_radius=1.0): +def _build_geometry(actx, ntargets, nsources, mode, target_radius=1.0): # source points t = np.linspace(0.0, 2.0 * np.pi, nsources, endpoint=False) sources = np.array([np.cos(t), np.sin(t)]) @@ -59,14 +56,14 @@ def _build_geometry(queue, ntargets, nsources, mode, target_radius=1.0): centers = (1.0 - radius) * targets expansion_radii = np.full(ntargets, radius) - return (cl.array.to_device(queue, targets), - cl.array.to_device(queue, sources), - vector_to_device(queue, centers), - cl.array.to_device(queue, expansion_radii), - cl.array.to_device(queue, sigma)) + return (actx.from_numpy(targets), + actx.from_numpy(sources), + actx.from_numpy(centers), + actx.from_numpy(expansion_radii), + actx.from_numpy(sigma)) -def _build_subset_indices(queue, ntargets, nsources, factor): +def _build_subset_indices(actx, ntargets, nsources, factor): tgtindices = np.arange(0, ntargets) srcindices = np.arange(0, nsources) @@ -82,17 +79,19 @@ def _build_subset_indices(queue, ntargets, nsources, factor): tgtindices, srcindices = np.meshgrid(tgtindices, srcindices) return ( - cl.array.to_device(queue, tgtindices.ravel()).with_queue(None), - cl.array.to_device(queue, srcindices.ravel()).with_queue(None)) + actx.freeze(actx.from_numpy(tgtindices.ravel())), + actx.freeze(actx.from_numpy(srcindices.ravel()))) +# {{{ test_qbx_direct + @pytest.mark.parametrize("factor", [1.0, 0.6]) @pytest.mark.parametrize("lpot_id", [1, 2]) -def test_qbx_direct(ctx_factory, factor, lpot_id): - logging.basicConfig(level=logging.INFO) +def test_qbx_direct(actx_factory, factor, lpot_id, visualize=False): + if visualize: + logging.basicConfig(level=logging.INFO) - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) + actx = actx_factory() ndim = 2 order = 12 @@ -106,79 +105,88 @@ def test_qbx_direct(ctx_factory, factor, lpot_id): base_knl = LaplaceKernel(ndim) knl = DirectionalSourceDerivative(base_knl, dir_vec_name="dsource_vec") else: - raise ValueError("unknow lpot_id") + raise ValueError(f"unknown lpot_id: {lpot_id}") from sumpy.expansion.local import LineTaylorLocalExpansion expn = LineTaylorLocalExpansion(knl, order) from sumpy.qbx import LayerPotential - lpot = LayerPotential(ctx, expansion=expn, source_kernels=(knl,), + lpot = LayerPotential(actx.context, expansion=expn, source_kernels=(knl,), target_kernels=(base_knl,)) from sumpy.qbx import LayerPotentialMatrixGenerator - mat_gen = LayerPotentialMatrixGenerator(ctx, expansion=expn, - source_kernels=(knl,), target_kernels=(base_knl,)) + mat_gen = LayerPotentialMatrixGenerator(actx.context, + expansion=expn, + source_kernels=(knl,), + target_kernels=(base_knl,)) from sumpy.qbx import LayerPotentialMatrixSubsetGenerator - blk_gen = LayerPotentialMatrixSubsetGenerator(ctx, expansion=expn, - source_kernels=(knl,), target_kernels=(base_knl,)) + blk_gen = LayerPotentialMatrixSubsetGenerator(actx.context, + expansion=expn, + source_kernels=(knl,), + target_kernels=(base_knl,)) for n in [200, 300, 400]: targets, sources, centers, expansion_radii, sigma = \ - _build_geometry(queue, n, n, mode_nr, target_radius=1.2) + _build_geometry(actx, n, n, mode_nr, target_radius=1.2) h = 2 * np.pi / n strengths = (sigma * h,) - tgtindices, srcindices = _build_subset_indices(queue, + tgtindices, srcindices = _build_subset_indices(actx, ntargets=n, nsources=n, factor=factor) extra_kwargs = {} if lpot_id == 2: from pytools.obj_array import make_obj_array - extra_kwargs["dsource_vec"] = \ - vector_to_device(queue, make_obj_array(np.ones((ndim, n)))) + extra_kwargs["dsource_vec"] = ( + actx.from_numpy(make_obj_array(np.ones((ndim, n)))) + ) - _, (result_lpot,) = lpot(queue, + _, (result_lpot,) = lpot(actx.queue, targets=targets, sources=sources, centers=centers, expansion_radii=expansion_radii, strengths=strengths, **extra_kwargs) - result_lpot = result_lpot.get() + result_lpot = actx.to_numpy(result_lpot) - _, (mat,) = mat_gen(queue, + _, (mat,) = mat_gen(actx.queue, targets=targets, sources=sources, centers=centers, expansion_radii=expansion_radii, **extra_kwargs) - mat = mat.get() - result_mat = mat.dot(strengths[0].get()) + mat = actx.to_numpy(mat) + result_mat = mat @ actx.to_numpy(strengths[0]) - _, (blk,) = blk_gen(queue, + _, (blk,) = blk_gen(actx.queue, targets=targets, sources=sources, centers=centers, expansion_radii=expansion_radii, tgtindices=tgtindices, srcindices=srcindices, **extra_kwargs) - blk = blk.get() + blk = actx.to_numpy(blk) - tgtindices = tgtindices.get(queue) - srcindices = srcindices.get(queue) + tgtindices = actx.to_numpy(tgtindices) + srcindices = actx.to_numpy(srcindices) eps = 1.0e-10 * la.norm(result_lpot) assert la.norm(result_mat - result_lpot) < eps assert la.norm(blk - mat[tgtindices, srcindices]) < eps +# }}} + + +# {{{ test_p2p_direct @pytest.mark.parametrize("exclude_self", [True, False]) @pytest.mark.parametrize("factor", [1.0, 0.6]) @pytest.mark.parametrize("lpot_id", [1, 2]) -def test_p2p_direct(ctx_factory, exclude_self, factor, lpot_id): - logging.basicConfig(level=logging.INFO) +def test_p2p_direct(actx_factory, exclude_self, factor, lpot_id, visualize=False): + if visualize: + logging.basicConfig(level=logging.INFO) - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) + actx = actx_factory() ndim = 2 mode_nr = 25 @@ -190,70 +198,73 @@ def test_p2p_direct(ctx_factory, exclude_self, factor, lpot_id): lknl = LaplaceKernel(ndim) lknl = DirectionalSourceDerivative(lknl, dir_vec_name="dsource_vec") else: - raise ValueError("unknow lpot_id") + raise ValueError(f"unknown lpot_id: '{lpot_id}'") from sumpy.p2p import P2P - lpot = P2P(ctx, [lknl], exclude_self=exclude_self) + lpot = P2P(actx.context, [lknl], exclude_self=exclude_self) from sumpy.p2p import P2PMatrixGenerator - mat_gen = P2PMatrixGenerator(ctx, [lknl], exclude_self=exclude_self) + mat_gen = P2PMatrixGenerator(actx.context, [lknl], exclude_self=exclude_self) from sumpy.p2p import P2PMatrixSubsetGenerator - blk_gen = P2PMatrixSubsetGenerator(ctx, [lknl], exclude_self=exclude_self) + blk_gen = P2PMatrixSubsetGenerator( + actx.context, [lknl], exclude_self=exclude_self) for n in [200, 300, 400]: - targets, sources, _, _, sigma = \ - _build_geometry(queue, n, n, mode_nr, target_radius=1.2) + targets, sources, _, _, sigma = ( + _build_geometry(actx, n, n, mode_nr, target_radius=1.2)) h = 2 * np.pi / n strengths = (sigma * h,) - tgtindices, srcindices = _build_subset_indices(queue, + tgtindices, srcindices = _build_subset_indices(actx, ntargets=n, nsources=n, factor=factor) extra_kwargs = {} if exclude_self: - extra_kwargs["target_to_source"] = \ - cl.array.arange(queue, 0, n, dtype=np.int32) + extra_kwargs["target_to_source"] = ( + actx.from_numpy(np.arange(n, dtype=np.int32)) + ) if lpot_id == 2: from pytools.obj_array import make_obj_array - extra_kwargs["dsource_vec"] = \ - vector_to_device(queue, make_obj_array(np.ones((ndim, n)))) + extra_kwargs["dsource_vec"] = ( + actx.from_numpy(make_obj_array(np.ones((ndim, n))))) - _, (result_lpot,) = lpot(queue, + _, (result_lpot,) = lpot(actx.queue, targets=targets, sources=sources, strength=strengths, **extra_kwargs) - result_lpot = result_lpot.get() + result_lpot = actx.to_numpy(result_lpot) - _, (mat,) = mat_gen(queue, + _, (mat,) = mat_gen(actx.queue, targets=targets, sources=sources, **extra_kwargs) - mat = mat.get() - result_mat = mat.dot(strengths[0].get()) + mat = actx.to_numpy(mat) + result_mat = mat @ actx.to_numpy(strengths[0]) - _, (blk,) = blk_gen(queue, + _, (blk,) = blk_gen(actx.queue, targets=targets, sources=sources, tgtindices=tgtindices, srcindices=srcindices, **extra_kwargs) - blk = blk.get() + blk = actx.to_numpy(blk) - tgtindices = tgtindices.get(queue) - srcindices = srcindices.get(queue) + tgtindices = actx.to_numpy(tgtindices) + srcindices = actx.to_numpy(srcindices) eps = 1.0e-10 * la.norm(result_lpot) assert la.norm(result_mat - result_lpot) < eps assert la.norm(blk - mat[tgtindices, srcindices]) < eps +# }}} + # You can test individual routines by typing -# $ python test_kernels.py "test_p2p(cl.create_some_context)" +# $ python test_matrixgen.py 'test_p2p_direct(_acf, True, 1.0, 1, visualize=True)' if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from pytest import main - main([__file__]) + pytest.main([__file__]) # vim: fdm=marker From a3053d2f316a196a3ba27547dcaf088ebb466869 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 20 Aug 2022 20:47:35 +0300 Subject: [PATCH 5/8] port test_misc to arraycontext --- test/test_misc.py | 124 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 91 insertions(+), 33 deletions(-) diff --git a/test/test_misc.py b/test/test_misc.py index 9aeb7f124..5032e2fc9 100644 --- a/test/test_misc.py +++ b/test/test_misc.py @@ -20,6 +20,7 @@ THE SOFTWARE. """ +import pytest import sys from dataclasses import dataclass from typing import Any, Callable @@ -27,19 +28,33 @@ import numpy as np import numpy.linalg as la +from arraycontext import pytest_generate_tests_for_array_contexts +from sumpy.array_context import ( # noqa: F401 + PytestPyOpenCLArrayContextFactory, _acf) + import sumpy.toys as t import sumpy.symbolic as sym -import pytest -import pyopencl as cl # noqa: F401 -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) - -from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, - BiharmonicKernel, YukawaKernel, StokesletKernel, StressletKernel, - ElasticityKernel, LineOfCompressionKernel, ExpressionKernel) -from sumpy.expansion.diff_op import (make_identity_diff_op, gradient, - divergence, laplacian, concat, as_scalar_pde, curl, diff) +from sumpy.kernel import ( + LaplaceKernel, + HelmholtzKernel, + BiharmonicKernel, + YukawaKernel, + StokesletKernel, + StressletKernel, + ElasticityKernel, + LineOfCompressionKernel, + ExpressionKernel) +from sumpy.expansion.diff_op import ( + make_identity_diff_op, concat, as_scalar_pde, diff, + gradient, divergence, laplacian, curl) + +import logging +logger = logging.getLogger(__name__) + +pytest_generate_tests = pytest_generate_tests_for_array_contexts([ + PytestPyOpenCLArrayContextFactory, + ]) # {{{ pde check for kernels @@ -93,15 +108,17 @@ def nderivs(self): KernelInfo(LineOfCompressionKernel(3, 0), mu=5, nu=0.2), KernelInfo(LineOfCompressionKernel(3, 1), mu=5, nu=0.2), ]) -def test_pde_check_kernels(ctx_factory, knl_info, order=5): +def test_pde_check_kernels(actx_factory, knl_info, order=5): + actx = actx_factory() + dim = knl_info.kernel.dim - tctx = t.ToyContext(ctx_factory(), knl_info.kernel, + tctx = t.ToyContext(actx.context, knl_info.kernel, extra_source_kwargs=knl_info.extra_kwargs) - np.random.seed(17) + rng = np.random.default_rng(42) pt_src = t.PointSources( tctx, - np.random.rand(dim, 50) - 0.5, + rng.random(size=(dim, 50)) - 0.5, np.ones(50)) from pytools.convergence import EOCRecorder @@ -117,12 +134,14 @@ def test_pde_check_kernels(ctx_factory, knl_info, order=5): err = la.norm(pde) eoc_rec.add_data_point(h, err) - print(eoc_rec) + logger.info("eoc:\n%s", eoc_rec) assert eoc_rec.order_estimate() > order - knl_info.nderivs + 1 - 0.1 # }}} +# {{{ test_pde_check + @pytest.mark.parametrize("dim", [1, 2, 3]) def test_pde_check(dim, order=4): from sumpy.point_calculus import CalculusPatch @@ -138,15 +157,19 @@ def test_pde_check(dim, order=4): err = la.norm(df_num-df_true) eoc_rec.add_data_point(h, err) - print(eoc_rec) + logger.info("eoc:\n%s", eoc_rec) assert eoc_rec.order_estimate() > order-2-0.1 +# }}} + +# {{{ test_order_finder + +@dataclass(frozen=True) class FakeTree: - def __init__(self, dimensions, root_extent, stick_out_factor): - self.dimensions = dimensions - self.root_extent = root_extent - self.stick_out_factor = stick_out_factor + dimensions: int + root_extent: float + stick_out_factor: float @pytest.mark.parametrize("knl", [ @@ -161,7 +184,7 @@ def test_order_finder(knl): orders = [ ofind(knl, frozenset([("k", 5)]), tree, level) for level in range(30)] - print(orders) + logger.info("orders: %s", orders) # Order should not increase with level assert (np.diff(orders) <= 0).all() @@ -180,11 +203,13 @@ def test_fmmlib_order_finder(knl): orders = [ ofind(knl, frozenset([("k", 5)]), tree, level) for level in range(30)] - print(orders) + logger.info("orders: %s", orders) # Order should not increase with level assert (np.diff(orders) <= 0).all() +# }}} + # {{{ expansion toys p2e2e2p test cases @@ -193,7 +218,7 @@ def approx_convergence_factor(orders, errors): return np.exp(poly[0]) -@dataclass +@dataclass(frozen=True) class P2E2E2PTestCase: source: np.ndarray target: np.ndarray @@ -243,12 +268,14 @@ def dim(self): # }}} +# {{{ test_toy_p2e2e2p + ORDERS_P2E2E2P = (3, 4, 5) RTOL_P2E2E2P = 1e-2 @pytest.mark.parametrize("case", P2E2E2P_TEST_CASES) -def test_toy_p2e2e2p(ctx_factory, case): +def test_toy_p2e2e2p(actx_factory, case): dim = case.dim src = case.source.reshape(dim, -1) @@ -269,8 +296,8 @@ def test_toy_p2e2e2p(ctx_factory, case): from sumpy.expansion import VolumeTaylorExpansionFactory - cl_ctx = ctx_factory() - ctx = t.ToyContext(cl_ctx, + actx = actx_factory() + ctx = t.ToyContext(actx.context, LaplaceKernel(dim), expansion_factory=VolumeTaylorExpansionFactory()) @@ -289,6 +316,10 @@ def test_toy_p2e2e2p(ctx_factory, case): assert conv_factor <= min(1, case_conv_factor * (1 + RTOL_P2E2E2P)), \ (conv_factor, case_conv_factor * (1 + RTOL_P2E2E2P)) +# }}} + + +# {{{ test_cse_matvec def test_cse_matvec(): from sumpy.expansion import CSEMatVecOperator @@ -309,16 +340,21 @@ def test_cse_matvec(): op = CSEMatVecOperator(input_coeffs, output_coeffs, shape=(4, 2)) m = np.array([[2, 0], [6, 0], [0, 1], [30, 16]]) - vec = np.random.random(2) + rng = np.random.default_rng(42) + vec = rng.random(2) expected_result = m @ vec actual_result = op.matvec(vec) assert np.allclose(expected_result, actual_result) - vec = np.random.random(4) + vec = rng.random(4) expected_result = m.T @ vec actual_result = op.transpose_matvec(vec) assert np.allclose(expected_result, actual_result) +# }}} + + +# {{{ test_diff_op_stokes def test_diff_op_stokes(): from sumpy.symbolic import symbols, Function @@ -341,6 +377,10 @@ def test_diff_op_stokes(): assert expected_output == actual_output +# }}} + + +# {{{ test_as_scalar_pde_stokes def test_as_scalar_pde_stokes(): diff_op = make_identity_diff_op(3, 4) @@ -350,13 +390,17 @@ def test_as_scalar_pde_stokes(): # velocity components in Stokes should satisfy Biharmonic for i in range(3): - print(as_scalar_pde(pde, i)) - print(laplacian(laplacian(u[i]))) + logger.info("pde\n%s", as_scalar_pde(pde, i)) + logger.info("\n%s", laplacian(laplacian(u[i]))) assert as_scalar_pde(pde, i) == laplacian(laplacian(u[0])) # pressure should satisfy Laplace assert as_scalar_pde(pde, 3) == laplacian(u[0]) +# }}} + + +# {{{ test_as_scalar_pde_maxwell def test_as_scalar_pde_maxwell(): from sumpy.symbolic import symbols @@ -374,6 +418,10 @@ def test_as_scalar_pde_maxwell(): assert as_scalar_pde(pde, i) == \ laplacian(op[0]) - mu*epsilon*diff(diff(op[0], t), t) +# }}} + + +# {{{ test_as_scalar_pde_elasticity def test_as_scalar_pde_elasticity(): @@ -408,6 +456,10 @@ def test_as_scalar_pde_elasticity(): assert scalar_pde == laplacian(laplacian(diff_op[0])) assert scalar_pde.order == 4 +# }}} + + +# {{{ test_elasticity_new def test_elasticity_new(): from pickle import dumps, loads @@ -424,6 +476,10 @@ def test_elasticity_new(): assert not isinstance(knl, StokesletKernel) assert loads(dumps(knl)) == knl +# }}} + + +# {{{ test_weird_kernel w = make_identity_diff_op(2) @@ -457,15 +513,17 @@ def get_pde_as_diff_op(self): assert fft_size == order +# }}} + # You can test individual routines by typing -# $ python test_misc.py 'test_p2p(cl.create_some_context)' +# $ python test_misc.py 'test_pde_check_kernels(_acf, +# KernelInfo(HelmholtzKernel(2), k=5), order=5)' if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from pytest import main - main([__file__]) + pytest.main([__file__]) # vim: fdm=marker From da2104f5cd3ed81ca7ab26b7ed6d114f5922b1a6 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 20 Aug 2022 20:54:36 +0300 Subject: [PATCH 6/8] port test_qbx to arraycontext --- test/test_qbx.py | 78 +++++++++++++++++++++++++++++++----------------- 1 file changed, 50 insertions(+), 28 deletions(-) diff --git a/test/test_qbx.py b/test/test_qbx.py index 379ef0e22..c2b61c3b5 100644 --- a/test/test_qbx.py +++ b/test/test_qbx.py @@ -20,33 +20,41 @@ THE SOFTWARE. """ -import numpy as np +import pytest import sys -import pyopencl as cl -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) +import numpy as np + +from arraycontext import pytest_generate_tests_for_array_contexts +from sumpy.array_context import ( # noqa: F401 + PytestPyOpenCLArrayContextFactory, _acf) + +from sumpy.expansion.local import ( + LineTaylorLocalExpansion, + VolumeTaylorLocalExpansion) import logging logger = logging.getLogger(__name__) -from sumpy.expansion.local import ( - LineTaylorLocalExpansion, VolumeTaylorLocalExpansion) -import pytest +pytest_generate_tests = pytest_generate_tests_for_array_contexts([ + PytestPyOpenCLArrayContextFactory, + ]) + + +# {{{ test_direct_qbx_vs_eigval @pytest.mark.parametrize("expn_class", [ LineTaylorLocalExpansion, VolumeTaylorLocalExpansion, ]) -def test_direct_qbx_vs_eigval(ctx_factory, expn_class): +def test_direct_qbx_vs_eigval(actx_factory, expn_class, visualize=False): """This evaluates a single layer potential on a circle using a known eigenvalue/eigenvector combination. """ + if visualize: + logging.basicConfig(level=logging.INFO) - logging.basicConfig(level=logging.INFO) - - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) + actx = actx_factory() from sumpy.kernel import LaplaceKernel lknl = LaplaceKernel(2) @@ -55,8 +63,10 @@ def test_direct_qbx_vs_eigval(ctx_factory, expn_class): from sumpy.qbx import LayerPotential - lpot = LayerPotential(ctx, expansion=expn_class(lknl, order), - target_kernels=(lknl,), source_kernels=(lknl,)) + lpot = LayerPotential(actx.context, + expansion=expn_class(lknl, order), + target_kernels=(lknl,), + source_kernels=(lknl,)) mode_nr = 25 @@ -85,30 +95,36 @@ def test_direct_qbx_vs_eigval(ctx_factory, expn_class): expansion_radii = np.ones(n) * radius strengths = (sigma * h,) - evt, (result_qbx,) = lpot(queue, targets, sources, centers, strengths, + evt, (result_qbx,) = lpot( + actx.queue, + targets, sources, centers, strengths, expansion_radii=expansion_radii) eocrec.add_data_point(h, np.max(np.abs(result_ref - result_qbx))) - print(eocrec) + logger.info("eoc:\n%s", eocrec) slack = 1.5 assert eocrec.order_estimate() > order - slack +# }}} + + +# {{{ test_direct_qbx_vs_eigval_with_tgt_deriv @pytest.mark.parametrize("expn_class", [ LineTaylorLocalExpansion, VolumeTaylorLocalExpansion, ]) -def test_direct_qbx_vs_eigval_with_tgt_deriv(ctx_factory, expn_class): +def test_direct_qbx_vs_eigval_with_tgt_deriv( + actx_factory, expn_class, visualize=False): """This evaluates a single layer potential on a circle using a known eigenvalue/eigenvector combination. """ + if visualize: + logging.basicConfig(level=logging.INFO) - logging.basicConfig(level=logging.INFO) - - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) + actx = actx_factory() from sumpy.kernel import LaplaceKernel, AxisTargetDerivative lknl = LaplaceKernel(2) @@ -117,9 +133,9 @@ def test_direct_qbx_vs_eigval_with_tgt_deriv(ctx_factory, expn_class): from sumpy.qbx import LayerPotential - lpot_dx = LayerPotential(ctx, expansion=expn_class(lknl, order), + lpot_dx = LayerPotential(actx.context, expansion=expn_class(lknl, order), target_kernels=(AxisTargetDerivative(0, lknl),), source_kernels=(lknl,)) - lpot_dy = LayerPotential(ctx, expansion=expn_class(lknl, order), + lpot_dy = LayerPotential(actx.context, expansion=expn_class(lknl, order), target_kernels=(AxisTargetDerivative(1, lknl),), source_kernels=(lknl,)) mode_nr = 15 @@ -151,9 +167,11 @@ def test_direct_qbx_vs_eigval_with_tgt_deriv(ctx_factory, expn_class): strengths = (sigma * h,) - evt, (result_qbx_dx,) = lpot_dx(queue, targets, sources, centers, strengths, + evt, (result_qbx_dx,) = lpot_dx(actx.queue, + targets, sources, centers, strengths, expansion_radii=expansion_radii) - evt, (result_qbx_dy,) = lpot_dy(queue, targets, sources, centers, strengths, + evt, (result_qbx_dy,) = lpot_dy(actx.queue, + targets, sources, centers, strengths, expansion_radii=expansion_radii) normals = unit_circle @@ -162,17 +180,21 @@ def test_direct_qbx_vs_eigval_with_tgt_deriv(ctx_factory, expn_class): eocrec.add_data_point(h, np.max(np.abs(result_ref - result_qbx))) if expn_class is not LineTaylorLocalExpansion: - print(eocrec) + logger.info("eoc:\n%s", eocrec) slack = 1.5 assert eocrec.order_estimate() > order - slack +# }}} + + +# You can test individual routines by typing +# $ python test_qbx.py 'test_direct_qbx_vs_eigval(_acf, LineTaylorLocalExpansion)' if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from pytest import main - main([__file__]) + pytest.main([__file__]) # vim: fdm=marker From 43e13a771d49a3b759261a65dbd6eecf66498492 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 20 Aug 2022 20:54:42 +0300 Subject: [PATCH 7/8] port test_tools to arraycontext --- test/test_tools.py | 71 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 53 insertions(+), 18 deletions(-) diff --git a/test/test_tools.py b/test/test_tools.py index 08a90689a..09c88b78e 100644 --- a/test/test_tools.py +++ b/test/test_tools.py @@ -20,26 +20,38 @@ THE SOFTWARE. """ -import logging -logger = logging.getLogger(__name__) +import pytest +import sys -import sumpy.symbolic as sym -from sumpy.tools import (fft_toeplitz_upper_triangular, - matvec_toeplitz_upper_triangular, loopy_fft, fft) import numpy as np -import pyopencl as cl -import pyopencl.array as cla -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests) +from arraycontext import pytest_generate_tests_for_array_contexts +from sumpy.array_context import ( # noqa: F401 + PytestPyOpenCLArrayContextFactory, _acf) -import pytest +import sumpy.symbolic as sym +from sumpy.tools import ( + fft_toeplitz_upper_triangular, + matvec_toeplitz_upper_triangular, + loopy_fft, + fft) +import logging +logger = logging.getLogger(__name__) + +pytest_generate_tests = pytest_generate_tests_for_array_contexts([ + PytestPyOpenCLArrayContextFactory, + ]) + + +# {{{ test_matvec_fft def test_matvec_fft(): k = 5 - v = np.random.rand(k) - x = np.random.rand(k) + + rng = np.random.default_rng(42) + v = rng.random(k) + x = rng.random(k) fft = fft_toeplitz_upper_triangular(v, x) matvec = matvec_toeplitz_upper_triangular(v, x) @@ -47,6 +59,10 @@ def test_matvec_fft(): for i in range(k): assert abs(fft[i] - matvec[i]) < 1e-14 +# }}} + + +# {{{ test_matvec_fft_small_floats def test_matvec_fft_small_floats(): k = 5 @@ -60,15 +76,34 @@ def test_matvec_fft_small_floats(): continue assert abs(f) > 1e-10 +# }}} + + +# {{{ test_fft @pytest.mark.parametrize("size", [1, 2, 7, 10, 30, 210]) -def test_fft(ctx_factory, size): - ctx = ctx_factory() - queue = cl.CommandQueue(ctx) +def test_fft(actx_factory, size): + actx = actx_factory() + inp = np.arange(size, dtype=np.complex64) - inp_dev = cla.to_device(queue, inp) + inp_dev = actx.from_numpy(inp) out = fft(inp) fft_func = loopy_fft(inp.shape, inverse=False, complex_dtype=inp.dtype.type) - evt, (out_dev,) = fft_func(queue, y=inp_dev) - assert np.allclose(out_dev.get(), out) + evt, (out_dev,) = fft_func(actx.queue, y=inp_dev) + + assert np.allclose(actx.to_numpy(out_dev), out) + +# }}} + + +# You can test individual routines by typing +# $ python test_tools.py 'test_fft(_acf, 30)' + +if __name__ == "__main__": + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + pytest.main([__file__]) + +# vim: fdm=marker From fa276fd0dbebb0b74987b27de89912e866f05456 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 20 Aug 2022 21:22:05 +0300 Subject: [PATCH 8/8] port examples to arraycontext --- examples/curve-pot.py | 82 ++++++++++++++++++++-------------- examples/expansion-toys.py | 35 ++++++++++----- examples/sym-exp-complexity.py | 19 +++++--- 3 files changed, 85 insertions(+), 51 deletions(-) diff --git a/examples/curve-pot.py b/examples/curve-pot.py index 810990ad1..8d1a91a21 100644 --- a/examples/curve-pot.py +++ b/examples/curve-pot.py @@ -1,16 +1,22 @@ -import pyopencl as cl import numpy as np import numpy.linalg as la +import pyopencl as cl + try: import matplotlib.pyplot as plt -except ModuleNotFoundError: - plt = None + USE_MATPLOTLIB = True +except ImportError: + USE_MATPLOTLIB = False try: from mayavi import mlab -except ModuleNotFoundError: - mlab = None + USE_MAYAVI = True +except ImportError: + USE_MAYAVI = False + +import logging +logging.basicConfig(level=logging.INFO) def process_kernel(knl, what_operator): @@ -45,17 +51,16 @@ def draw_pot_figure(aspect_ratio, ovsmp_center_exp=0.66, force_center_side=None): - import logging - logging.basicConfig(level=logging.INFO) - if novsmp is None: novsmp = 4*nsrc if what_operator_lpot is None: what_operator_lpot = what_operator + from sumpy.array_context import PyOpenCLArrayContext ctx = cl.create_some_context() queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue, force_device_scalars=True) # {{{ make plot targets @@ -86,7 +91,7 @@ def draw_pot_figure(aspect_ratio, knl_kwargs = {} vol_source_knl, vol_target_knl = process_kernel(knl, what_operator) - p2p = P2P(ctx, source_kernels=(vol_source_knl,), + p2p = P2P(actx.context, source_kernels=(vol_source_knl,), target_kernels=(vol_target_knl,), exclude_self=False, value_dtypes=np.complex128) @@ -94,8 +99,10 @@ def draw_pot_figure(aspect_ratio, lpot_source_knl, lpot_target_knl = process_kernel(knl, what_operator_lpot) from sumpy.qbx import LayerPotential - lpot = LayerPotential(ctx, expansion=expn_class(knl, order=order), - source_kernels=(lpot_source_knl,), target_kernels=(lpot_target_knl,), + lpot = LayerPotential(actx.context, + expansion=expn_class(knl, order=order), + source_kernels=(lpot_source_knl,), + target_kernels=(lpot_target_knl,), value_dtypes=np.complex128) # }}} @@ -142,8 +149,9 @@ def map_to_curve(t): + center_side[:, np.newaxis] * center_dist*native_curve.normal) - #native_curve.plot() - #plt.show() + if 0: + native_curve.plot() + plt.show() volpot_kwargs = knl_kwargs.copy() lpot_kwargs = knl_kwargs.copy() @@ -169,7 +177,9 @@ def map_to_curve(t): def apply_lpot(x): xovsmp = np.dot(fim, x) - evt, (y,) = lpot(queue, native_curve.pos, ovsmp_curve.pos, + evt, (y,) = lpot(actx.queue, + native_curve.pos, + ovsmp_curve.pos, centers, [xovsmp * ovsmp_curve.speed * ovsmp_weights], expansion_radii=np.ones(centers.shape[1]), @@ -191,10 +201,14 @@ def apply_lpot(x): mode_nr = 0 density = np.cos(mode_nr*2*np.pi*native_t).astype(np.complex128) ovsmp_density = np.cos(mode_nr*2*np.pi*ovsmp_t).astype(np.complex128) - evt, (vol_pot,) = p2p(queue, fp.points, native_curve.pos, + evt, (vol_pot,) = p2p(actx.queue, + fp.points, + native_curve.pos, [native_curve.speed*native_weights*density], **volpot_kwargs) - evt, (curve_pot,) = lpot(queue, native_curve.pos, ovsmp_curve.pos, + evt, (curve_pot,) = lpot(actx.queue, + native_curve.pos, + ovsmp_curve.pos, centers, [ovsmp_density * ovsmp_curve.speed * ovsmp_weights], expansion_radii=np.ones(centers.shape[1]), @@ -202,7 +216,7 @@ def apply_lpot(x): # }}} - if 0: + if USE_MATPLOTLIB: # {{{ plot on-surface potential in 2D plt.plot(curve_pot, label="pot") @@ -216,7 +230,7 @@ def apply_lpot(x): ("potential", vol_pot.real) ]) - if 0: + if USE_MATPLOTLIB: # {{{ 2D false-color plot plt.clf() @@ -230,12 +244,8 @@ def apply_lpot(x): # close the curve plt.plot(src[-1::-len(src)+1, 0], src[-1::-len(src)+1, 1], "o-k") - #plt.gca().set_aspect("equal", "datalim") cb = plt.colorbar(shrink=0.9) cb.set_label(r"$\log_{10}(\mathdefault{Error})$") - #from matplotlib.ticker import NullFormatter - #plt.gca().xaxis.set_major_formatter(NullFormatter()) - #plt.gca().yaxis.set_major_formatter(NullFormatter()) fp.set_matplotlib_limits() # }}} @@ -261,7 +271,7 @@ def apply_lpot(x): plotval_vol[outlier_flag] = sum( nb[outlier_flag] for nb in neighbors)/len(neighbors) - if mlab is not None: + if USE_MAYAVI: fp.show_scalar_in_mayavi(scale*plotval_vol, max_val=1) mlab.colorbar() if 1: @@ -275,17 +285,23 @@ def apply_lpot(x): if __name__ == "__main__": - draw_pot_figure(aspect_ratio=1, nsrc=100, novsmp=100, helmholtz_k=(35+4j)*0.3, + draw_pot_figure( + aspect_ratio=1, nsrc=100, novsmp=100, helmholtz_k=(35+4j)*0.3, what_operator="D", what_operator_lpot="D", force_center_side=1) + if USE_MATPLOTLIB: + plt.savefig("eigvals-ext-nsrc100-novsmp100.pdf") + plt.clf() -# plt.savefig("eigvals-ext-nsrc100-novsmp100.pdf") - #plt.clf() - #draw_pot_figure(aspect_ratio=1, nsrc=100, novsmp=100, helmholtz_k=0, - # what_operator="D", what_operator_lpot="D", force_center_side=-1) - #plt.savefig("eigvals-int-nsrc100-novsmp100.pdf") - #plt.clf() - #draw_pot_figure(aspect_ratio=1, nsrc=100, novsmp=200, helmholtz_k=0, - # what_operator="D", what_operator_lpot="D", force_center_side=-1) - #plt.savefig("eigvals-int-nsrc100-novsmp200.pdf") + # draw_pot_figure( + # aspect_ratio=1, nsrc=100, novsmp=100, helmholtz_k=0, + # what_operator="D", what_operator_lpot="D", force_center_side=-1) + # plt.savefig("eigvals-int-nsrc100-novsmp100.pdf") + # plt.clf() + + # draw_pot_figure( + # aspect_ratio=1, nsrc=100, novsmp=200, helmholtz_k=0, + # what_operator="D", what_operator_lpot="D", force_center_side=-1) + # plt.savefig("eigvals-int-nsrc100-novsmp200.pdf") + # plt.clf() # vim: fdm=marker diff --git a/examples/expansion-toys.py b/examples/expansion-toys.py index e774b17ae..48647d0ee 100644 --- a/examples/expansion-toys.py +++ b/examples/expansion-toys.py @@ -1,21 +1,32 @@ +import numpy as np + import pyopencl as cl + import sumpy.toys as t -import numpy as np from sumpy.visualization import FieldPlotter +from sumpy.kernel import ( # noqa: F401 + YukawaKernel, + HelmholtzKernel, + LaplaceKernel) + try: import matplotlib.pyplot as plt -except ModuleNotFoundError: - plt = None + USE_MATPLOTLIB = True +except ImportError: + USE_MATPLOTLIB = False def main(): - from sumpy.kernel import ( # noqa: F401 - YukawaKernel, HelmholtzKernel, LaplaceKernel) + from sumpy.array_context import PyOpenCLArrayContext + ctx = cl.create_some_context() + queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue, force_device_scalars=True) + tctx = t.ToyContext( - cl.create_some_context(), - #LaplaceKernel(2), + actx.context, + # LaplaceKernel(2), YukawaKernel(2), extra_kernel_kwargs={"lam": 5}, - #HelmholtzKernel(2), extra_kernel_kwargs={"k": 0.3}, + # HelmholtzKernel(2), extra_kernel_kwargs={"k": 0.3}, ) pt_src = t.PointSources( @@ -25,7 +36,7 @@ def main(): fp = FieldPlotter([3, 0], extent=8) - if 0 and plt is not None: + if USE_MATPLOTLIB: t.logplot(fp, pt_src, cmap="jet") plt.colorbar() plt.show() @@ -35,12 +46,12 @@ def main(): lexp = t.local_expand(mexp, [3, 0]) lexp2 = t.local_expand(lexp, [3, 1], 3) - #diff = mexp - pt_src - #diff = mexp2 - pt_src + # diff = mexp - pt_src + # diff = mexp2 - pt_src diff = lexp2 - pt_src print(t.l_inf(diff, 1.2, center=lexp2.center)) - if 1 and plt is not None: + if USE_MATPLOTLIB: t.logplot(fp, diff, cmap="jet", vmin=-3, vmax=0) plt.colorbar() plt.show() diff --git a/examples/sym-exp-complexity.py b/examples/sym-exp-complexity.py index ae91a932c..779bfc360 100644 --- a/examples/sym-exp-complexity.py +++ b/examples/sym-exp-complexity.py @@ -1,6 +1,8 @@ import numpy as np -import pyopencl as cl import loopy as lp + +import pyopencl as cl + from sumpy.kernel import LaplaceKernel, HelmholtzKernel from sumpy.expansion.local import ( LinearPDEConformingVolumeTaylorLocalExpansion, @@ -9,14 +11,19 @@ LinearPDEConformingVolumeTaylorMultipoleExpansion, ) from sumpy.e2e import E2EFromCSR + try: import matplotlib.pyplot as plt -except ModuleNotFoundError: - plt = None + USE_MATPLOTLIB = True +except ImportError: + USE_MATPLOTLIB = False def find_flops(): + from sumpy.array_context import PyOpenCLArrayContext ctx = cl.create_some_context() + queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue, force_device_scalars=True) if 0: knl = LaplaceKernel(2) @@ -35,7 +42,7 @@ def find_flops(): print(order) m_expn = m_expn_cls(knl, order) l_expn = l_expn_cls(knl, order) - m2l = E2EFromCSR(ctx, m_expn, l_expn) + m2l = E2EFromCSR(actx.context, m_expn, l_expn) loopy_knl = m2l.get_kernel() loopy_knl = lp.add_and_infer_dtypes( @@ -74,7 +81,7 @@ def plot_flops(): flops = [45, 194, 474, 931, 1650, 2632, 3925, 5591, 7706, 10272] filename = "helmholtz-m2l-complexity-2d.pdf" - if plt is not None: + if USE_MATPLOTLIB: plt.rc("font", size=16) plt.title(case) plt.ylabel("Flop count") @@ -86,5 +93,5 @@ def plot_flops(): if __name__ == "__main__": - #find_flops() + # find_flops() plot_flops()