From 749f37a4f34e67a91a89fc1dd5a39e3c89067229 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Fri, 17 Jul 2020 22:51:26 +0200 Subject: [PATCH] Revert "Merge branch 'test-matrix-refactor' into 'master'" This reverts merge request !207 --- pytential/linalg/proxy.py | 17 +- test/extra_int_eq_data.py | 23 -- test/extra_matrix_data.py | 112 ------- test/test_linalg_proxy.py | 141 +++------ test/test_matrix.py | 605 +++++++++++++++++++++++--------------- 5 files changed, 427 insertions(+), 471 deletions(-) delete mode 100644 test/extra_matrix_data.py diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py index 4183cd0a..0f79148e 100644 --- a/pytential/linalg/proxy.py +++ b/pytential/linalg/proxy.py @@ -39,7 +39,9 @@ Proxy Point Generation ~~~~~~~~~~~~~~~~~~~~~~ .. autoclass:: ProxyGenerator + .. autofunction:: partition_by_nodes +.. autofunction:: partition_from_coarse .. autofunction:: gather_block_neighbor_points .. autofunction:: gather_block_interaction_points @@ -48,14 +50,16 @@ Proxy Point Generation # {{{ point index partitioning -def partition_by_nodes(actx, discr, - tree_kind="adaptive-level-restricted", max_nodes_in_box=None): +def partition_by_nodes(actx, discr, use_tree=True, max_nodes_in_box=None): """Generate equally sized ranges of nodes. The partition is created at the lowest level of granularity, i.e. nodes. This results in balanced ranges of points, but will split elements across different ranges. :arg discr: a :class:`meshmode.discretization.Discretization`. - :arg tree_kind: if not *None*, it is passed to :class:`boxtree.TreeBuilder`. + :arg use_tree: if ``True``, node partitions are generated using a + :class:`boxtree.TreeBuilder`, which leads to geometrically close + points to belong to the same partition. If ``False``, a simple linear + partition is constructed. :arg max_nodes_in_box: passed to :class:`boxtree.TreeBuilder`. :return: a :class:`sumpy.tools.BlockIndexRanges`. @@ -65,7 +69,7 @@ def partition_by_nodes(actx, discr, # FIXME: this is just an arbitrary value max_nodes_in_box = 32 - if tree_kind is not None: + if use_tree: from boxtree import box_flags_enum from boxtree import TreeBuilder @@ -74,8 +78,7 @@ def partition_by_nodes(actx, discr, from meshmode.dof_array import flatten, thaw tree, _ = builder(actx.queue, flatten(thaw(actx, discr.nodes())), - max_particles_in_box=max_nodes_in_box, - kind=tree_kind) + max_particles_in_box=max_nodes_in_box) tree = tree.get(actx.queue) leaf_boxes, = (tree.box_flags @@ -96,7 +99,7 @@ def partition_by_nodes(actx, discr, ranges = actx.from_numpy(np.arange( 0, discr.ndofs + 1, - max_nodes_in_box, dtype=np.int)) + discr.ndofs // max_nodes_in_box, dtype=np.int)) assert ranges[-1] == discr.ndofs diff --git a/test/extra_int_eq_data.py b/test/extra_int_eq_data.py index dbc807a7..3942adba 100644 --- a/test/extra_int_eq_data.py +++ b/test/extra_int_eq_data.py @@ -357,29 +357,6 @@ class SphereTestCase(IntegralEquationTestCase): uniform_refinement_rounds=resolution) -class TorusTestCase(IntegralEquationTestCase): - name = "torus" - - # qbx - qbx_order = 4 - target_order = 7 - use_refinement = True - - # geometry - r_major = 10.0 - r_minor = 1.0 - - # test case - resolutions = [0, 1, 2] - - def get_mesh(self, resolution, mesh_order): - from meshmode.mesh.generation import generate_torus - mesh = generate_torus(self.r_major, self.r_minor, order=mesh_order) - - from meshmode.mesh.refinement import refine_uniformly - return refine_uniformly(mesh, resolution) - - class MergedCubesTestCase(Helmholtz3DTestCase): name = "merged_cubes" diff --git a/test/extra_matrix_data.py b/test/extra_matrix_data.py deleted file mode 100644 index 983c815e..00000000 --- a/test/extra_matrix_data.py +++ /dev/null @@ -1,112 +0,0 @@ -import numpy as np -import numpy.linalg as la - -from pytools.obj_array import make_obj_array -from sumpy.tools import BlockIndexRanges, MatrixBlockIndexRanges - -from pytential import sym - -import extra_int_eq_data as extra - - -# {{{ helpers - -def max_block_error(mat, blk, index_set, p=None): - error = -np.inf - for i in range(index_set.nblocks): - mat_i = index_set.take(mat, i) - blk_i = index_set.block_take(blk, i) - - error = max(error, la.norm(mat_i - blk_i, ord=p) / la.norm(mat_i, ord=p)) - - return error - -# }}} - - -# {{{ MatrixTestCase - -class MatrixTestCaseMixin: - # partitioning - approx_block_count = 10 - max_particles_in_box = None - tree_kind = "adaptive-level-restricted" - index_sparsity_factor = 1.0 - - # operators - op_type = "scalar" - - # disable fmm for matrix tests - fmm_backend = None - - def get_block_indices(self, actx, discr, matrix_indices=True): - max_particles_in_box = self.max_particles_in_box - if max_particles_in_box is None: - max_particles_in_box = discr.ndofs // self.approx_block_count - - from pytential.linalg.proxy import partition_by_nodes - indices = partition_by_nodes(actx, discr, - tree_kind=self.tree_kind, - max_nodes_in_box=max_particles_in_box) - - if abs(self.index_sparsity_factor - 1.0) < 1.0e-14: - if not matrix_indices: - return indices - return MatrixBlockIndexRanges(actx.context, indices, indices) - - # randomly pick a subset of points - indices = indices.get(actx.queue) - - subset = np.empty(indices.nblocks, dtype=np.object) - for i in range(indices.nblocks): - iidx = indices.block_indices(i) - isize = int(self.index_sparsity_factor * len(iidx)) - isize = max(1, min(isize, len(iidx))) - - subset[i] = np.sort(np.random.choice(iidx, size=isize, replace=False)) - - ranges = actx.from_numpy(np.cumsum([0] + [r.shape[0] for r in subset])) - indices = actx.from_numpy(np.hstack(subset)) - - indices = BlockIndexRanges(actx.context, - actx.freeze(indices), actx.freeze(ranges)) - - if not matrix_indices: - return indices - return MatrixBlockIndexRanges(actx.context, indices, indices) - - def get_operator(self, ambient_dim, qbx_forced_limit="avg"): - knl = self.knl_class(ambient_dim) - kwargs = self.knl_sym_kwargs.copy() - kwargs["qbx_forced_limit"] = qbx_forced_limit - - if self.op_type == "scalar": - sym_u = sym.var("u") - sym_op = sym.S(knl, sym_u, **kwargs) - elif self.op_type == "scalar_mixed": - sym_u = sym.var("u") - sym_op = sym.S(knl, 0.3 * sym_u, **kwargs) \ - + sym.D(knl, 0.5 * sym_u, **kwargs) - elif self.op_type == "vector": - sym_u = sym.make_sym_vector("u", ambient_dim) - - sym_op = make_obj_array([ - sym.Sp(knl, sym_u[0], **kwargs) - + sym.D(knl, sym_u[1], **kwargs), - sym.S(knl, 0.4 * sym_u[0], **kwargs) - + 0.3 * sym.D(knl, sym_u[0], **kwargs) - ]) - else: - raise ValueError(f"unknown operator type: '{self.op_type}'") - - return sym_u, sym_op - - -class CurveTestCase(MatrixTestCaseMixin, extra.CurveTestCase): - pass - - -class TorusTestCase(MatrixTestCaseMixin, extra.TorusTestCase): - pass - -# }}} diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py index 09a1367b..8e485d6e 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -22,46 +22,41 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from functools import partial - import numpy as np import numpy.linalg as la import pyopencl as cl +import pyopencl.array # noqa from pytential import bind, sym -from pytential import GeometryCollection from meshmode.array_context import PyOpenCLArrayContext -from meshmode.mesh.generation import ellipse +from meshmode.mesh.generation import ( # noqa + ellipse, NArmedStarfish, generate_torus, make_curve_mesh) import pytest from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) -import extra_matrix_data as extra -import logging -logger = logging.getLogger(__name__) -logging.basicConfig(level=logging.INFO) +from test_matrix import _build_geometry, _build_block_index -def plot_partition_indices(actx, discr, indices, **kwargs): - try: - import matplotlib.pyplot as pt - except ImportError: - return +def _plot_partition_indices(actx, discr, indices, **kwargs): + import matplotlib.pyplot as pt indices = indices.get(actx.queue) + args = [ - kwargs.get("tree_kind", "linear").replace("-", "_"), - kwargs.get("discr_stage", "stage1"), + kwargs.get("method", "unknown"), + "tree" if kwargs.get("use_tree", False) else "linear", + kwargs.get("pid", "stage1"), discr.ambient_dim ] pt.figure(figsize=(10, 8), dpi=300) pt.plot(np.diff(indices.ranges)) - pt.savefig("test_partition_{1}_{3}d_ranges_{2}.png".format(*args)) + pt.savefig("test_partition_{0}_{1}_{3}d_ranges_{2}.png".format(*args)) pt.clf() from pytential.utils import flatten_to_numpy @@ -69,18 +64,24 @@ def plot_partition_indices(actx, discr, indices, **kwargs): sources = flatten_to_numpy(actx, discr.nodes()) pt.figure(figsize=(10, 8), dpi=300) + if indices.indices.shape[0] != discr.ndofs: pt.plot(sources[0], sources[1], 'ko', alpha=0.5) - for i in range(indices.nblocks): isrc = indices.block_indices(i) pt.plot(sources[0][isrc], sources[1][isrc], 'o') pt.xlim([-1.5, 1.5]) pt.ylim([-1.5, 1.5]) - pt.savefig("test_partition_{1}_{3}d_{2}.png".format(*args)) + pt.savefig("test_partition_{0}_{1}_{3}d_{2}.png".format(*args)) pt.clf() elif discr.ambient_dim == 3: + from meshmode.discretization import NoninterpolatoryElementGroupError + try: + discr.groups[0].basis() + except NoninterpolatoryElementGroupError: + return + from meshmode.discretization.visualization import make_visualizer marker = -42.0 * np.ones(discr.ndofs) @@ -99,70 +100,38 @@ def plot_partition_indices(actx, discr, indices, **kwargs): ]) -PROXY_TEST_CASES = [ - extra.CurveTestCase( - name="ellipse", - target_order=7, - curve_fn=partial(ellipse, 3.0)), - extra.TorusTestCase( - target_order=2, - resolutions=[0]) - ] - - -@pytest.mark.skip(reason="only useful with visualize=True") -@pytest.mark.parametrize("tree_kind", ['adaptive', None]) -@pytest.mark.parametrize("case", PROXY_TEST_CASES) -def test_partition_points(ctx_factory, tree_kind, case, visualize=False): - """Tests that the points are correctly partitioned (by visualization).""" - +@pytest.mark.parametrize("use_tree", [True, False]) +@pytest.mark.parametrize("ambient_dim", [2, 3]) +def test_partition_points(ctx_factory, use_tree, ambient_dim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) - case = case.copy(tree_kind=tree_kind, index_sparsity_factor=0.6) - logger.info("\n%s", case) - - # {{{ - - qbx = case.get_layer_potential(actx, case.resolutions[-1], case.target_order) - places = GeometryCollection(qbx, auto_where=case.name) - - density_discr = places.get_discretization(case.name) - indices = case.get_block_indices(actx, density_discr) + places, dofdesc = _build_geometry(actx, ambient_dim=ambient_dim) + discr = places.get_discretization(dofdesc.geometry, dofdesc.discr_stage) + indices = _build_block_index(actx, discr, use_tree=use_tree, factor=0.6) if visualize: - plot_partition_indices(actx, density_discr, indices, tree_kind=tree_kind) + _plot_partition_indices(actx, discr, indices, use_tree=use_tree) - # }}} - - -@pytest.mark.parametrize("case", PROXY_TEST_CASES) -@pytest.mark.parametrize("index_sparsity_factor", [1.0, 0.6]) -def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=False): - """Tests that the proxies generated are all at the correct radius from the - points in the cluster, etc. - """ +@pytest.mark.parametrize("ambient_dim", [2, 3]) +@pytest.mark.parametrize("factor", [1.0, 0.6]) +def test_proxy_generator(ctx_factory, ambient_dim, factor, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) - case = case.copy(index_sparsity_factor=index_sparsity_factor) - logger.info("\n%s", case) - - # {{{ generate proxies + places, dofdesc = _build_geometry(actx, ambient_dim=ambient_dim) + dofdesc = dofdesc.to_stage1() - qbx = case.get_layer_potential(actx, case.resolutions[-1], case.target_order) - places = GeometryCollection(qbx, auto_where=case.name) - - density_discr = places.get_discretization(case.name) - srcindices = case.get_block_indices(actx, density_discr, matrix_indices=False) + density_discr = places.get_discretization(dofdesc.geometry, dofdesc.discr_stage) + srcindices = _build_block_index(actx, density_discr, factor=factor) from pytential.linalg.proxy import ProxyGenerator generator = ProxyGenerator(places) proxies, pxyranges, pxycenters, pxyradii = \ - generator(actx, places.auto_source, srcindices) + generator(actx, dofdesc, srcindices) proxies = np.vstack([actx.to_numpy(p) for p in proxies]) pxyranges = actx.to_numpy(pxyranges) @@ -175,13 +144,8 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal r = la.norm(proxies[:, ipxy] - pxycenters[:, i].reshape(-1, 1), axis=0) assert np.allclose(r - pxyradii[i], 0.0, atol=1.0e-14) - # }}} - - # {{{ visualization - srcindices = srcindices.get(queue) if visualize: - ambient_dim = places.ambient_dim if ambient_dim == 2: import matplotlib.pyplot as pt @@ -248,45 +212,32 @@ def test_proxy_generator(ctx_factory, case, index_sparsity_factor, visualize=Fal ambient_dim, i) vis.write_vtk_file(filename, []) - # }}} - - -@pytest.mark.parametrize("case", PROXY_TEST_CASES) -@pytest.mark.parametrize("index_sparsity_factor", [1.0, 0.6]) -def test_interaction_points(ctx_factory, - case, index_sparsity_factor, visualize=False): - """Test that neighboring points (inside the proxy balls, but outside the - current block/cluster) are actually inside. - """ +@pytest.mark.parametrize("ambient_dim", [2, 3]) +@pytest.mark.parametrize("factor", [1.0, 0.6]) +def test_interaction_points(ctx_factory, ambient_dim, factor, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) - case = case.copy(index_sparsity_factor=index_sparsity_factor) - logger.info("\n%s", case) - - # {{{ check neighboring points + places, dofdesc = _build_geometry(actx, ambient_dim=ambient_dim) + dofdesc = dofdesc.to_stage1() - qbx = case.get_layer_potential(actx, case.resolutions[-1], case.target_order) - places = GeometryCollection(qbx, auto_where=case.name) - - density_discr = places.get_discretization(case.name) - srcindices = case.get_block_indices(actx, density_discr, matrix_indices=False) + density_discr = places.get_discretization(dofdesc.geometry, dofdesc.discr_stage) + srcindices = _build_block_index(actx, density_discr, factor=factor) # generate proxy points from pytential.linalg.proxy import ProxyGenerator generator = ProxyGenerator(places) - _, _, pxycenters, pxyradii = generator(actx, places.auto_source, srcindices) + _, _, pxycenters, pxyradii = generator(actx, dofdesc, srcindices) - # get neighboring points from pytential.linalg.proxy import ( # noqa gather_block_neighbor_points, gather_block_interaction_points) nbrindices = gather_block_neighbor_points(actx, density_discr, srcindices, pxycenters, pxyradii) nodes, ranges = gather_block_interaction_points(actx, - places, places.auto_source, srcindices) + places, dofdesc, srcindices) srcindices = srcindices.get(queue) nbrindices = nbrindices.get(queue) @@ -297,13 +248,8 @@ def test_interaction_points(ctx_factory, assert not np.any(np.isin(inbr, isrc)) - # }}} - - # {{{ visualize - from pytential.utils import flatten_to_numpy if visualize: - ambient_dim = places.ambient_dim if ambient_dim == 2: import matplotlib.pyplot as pt density_nodes = flatten_to_numpy(actx, density_discr.nodes()) @@ -354,7 +300,6 @@ def test_interaction_points(ctx_factory, vis.write_vtk_file(filename, [ ("marker", marker_dev), ]) - # }}} if __name__ == "__main__": diff --git a/test/test_matrix.py b/test/test_matrix.py index 291cefc7..ce429642 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -31,40 +31,171 @@ import numpy as np import numpy.linalg as la import pyopencl as cl +import pyopencl.array +from pytools.obj_array import make_obj_array, is_obj_array + +from sumpy.tools import BlockIndexRanges, MatrixBlockIndexRanges from sumpy.symbolic import USE_SYMENGINE -from sumpy.tools import MatrixBlockIndexRanges -from pytools.obj_array import make_obj_array from pytential import bind, sym from pytential import GeometryCollection from meshmode.array_context import PyOpenCLArrayContext -from meshmode.mesh.generation import ellipse, NArmedStarfish +from meshmode.mesh.generation import ( # noqa + ellipse, NArmedStarfish, make_curve_mesh, generate_torus) import pytest from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) -import extra_matrix_data as extra -import logging -logger = logging.getLogger(__name__) -logging.basicConfig(level=logging.INFO) +try: + import matplotlib.pyplot as pt +except ImportError: + pass + + +def _build_geometry(actx, + ambient_dim=2, + nelements=30, + target_order=7, + qbx_order=4, + curve_f=None, + auto_where=None): + + if curve_f is None: + curve_f = NArmedStarfish(5, 0.25) + + if ambient_dim == 2: + mesh = make_curve_mesh(curve_f, + np.linspace(0, 1, nelements + 1), + target_order) + elif ambient_dim == 3: + mesh = generate_torus(10.0, 2.0, order=target_order) + else: + raise ValueError("unsupported ambient dimension") + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + from pytential.qbx import QBXLayerPotentialSource + density_discr = Discretization(actx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + qbx = QBXLayerPotentialSource(density_discr, + fine_order=4 * target_order, + qbx_order=qbx_order, + fmm_order=False) + places = GeometryCollection(qbx, auto_where=auto_where) + + return places, places.auto_source + + +def _build_block_index(actx, discr, + nblks=10, + factor=1.0, + use_tree=True): + max_particles_in_box = discr.ndofs // nblks + + # create index ranges + from pytential.linalg.proxy import partition_by_nodes + indices = partition_by_nodes(actx, discr, + use_tree=use_tree, + max_nodes_in_box=max_particles_in_box) + + if abs(factor - 1.0) < 1.0e-14: + return indices + + # randomly pick a subset of points + # FIXME: this needs porting in sumpy.tools.BlockIndexRanges + indices = indices.get(actx.queue) + + indices_ = np.empty(indices.nblocks, dtype=np.object) + for i in range(indices.nblocks): + iidx = indices.block_indices(i) + isize = int(factor * len(iidx)) + isize = max(1, min(isize, len(iidx))) + + indices_[i] = np.sort( + np.random.choice(iidx, size=isize, replace=False)) + + ranges_ = actx.from_numpy(np.cumsum([0] + [r.shape[0] for r in indices_])) + indices_ = actx.from_numpy(np.hstack(indices_)) + + indices = BlockIndexRanges(actx.context, + actx.freeze(indices_), actx.freeze(ranges_)) + + return indices + + +def _build_op(lpot_id, + k=0, + ambient_dim=2, + source=sym.DEFAULT_SOURCE, + target=sym.DEFAULT_TARGET, + qbx_forced_limit="avg"): + from sumpy.kernel import LaplaceKernel, HelmholtzKernel + + if k: + knl = HelmholtzKernel(ambient_dim) + knl_kwargs = {"k": k} + else: + knl = LaplaceKernel(ambient_dim) + knl_kwargs = {} + + lpot_kwargs = { + "qbx_forced_limit": qbx_forced_limit, + "source": source, + "target": target} + lpot_kwargs.update(knl_kwargs) + if lpot_id == 1: + # scalar single-layer potential + u_sym = sym.var("u") + op = sym.S(knl, u_sym, **lpot_kwargs) + elif lpot_id == 2: + # scalar combination of layer potentials + u_sym = sym.var("u") + op = sym.S(knl, 0.3 * u_sym, **lpot_kwargs) \ + + sym.D(knl, 0.5 * u_sym, **lpot_kwargs) + elif lpot_id == 3: + # vector potential + u_sym = sym.make_sym_vector("u", 2) + u0_sym, u1_sym = u_sym + + op = make_obj_array([ + sym.Sp(knl, u0_sym, **lpot_kwargs) + + sym.D(knl, u1_sym, **lpot_kwargs), + sym.S(knl, 0.4 * u0_sym, **lpot_kwargs) + + 0.3 * sym.D(knl, u0_sym, **lpot_kwargs) + ]) + else: + raise ValueError("Unknown lpot_id: {}".format(lpot_id)) + + op = 0.5 * u_sym + op + + return op, u_sym, knl_kwargs + + +def _max_block_error(mat, blk, index_set): + error = -np.inf + for i in range(index_set.nblocks): + mat_i = index_set.take(mat, i) + blk_i = index_set.block_take(blk, i) + + error = max(error, la.norm(mat_i - blk_i) / la.norm(mat_i)) + + return error @pytest.mark.skipif(USE_SYMENGINE, reason="https://gitlab.tiker.net/inducer/sumpy/issues/25") @pytest.mark.parametrize("k", [0, 42]) -@pytest.mark.parametrize("curve_fn", [ +@pytest.mark.parametrize("curve_f", [ partial(ellipse, 3), NArmedStarfish(5, 0.25)]) -@pytest.mark.parametrize("op_type", ["scalar_mixed", "vector"]) -def test_build_matrix(ctx_factory, k, curve_fn, op_type, visualize=False): - """Checks that the matrix built with `symbolic.execution.build_matrix` - gives the same (to tolerance) answer as a direct evaluation. - """ - +@pytest.mark.parametrize("lpot_id", [2, 3]) +def test_matrix_build(ctx_factory, k, curve_f, lpot_id, visualize=False): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) @@ -73,63 +204,48 @@ def test_build_matrix(ctx_factory, k, curve_fn, op_type, visualize=False): from sympy.core.cache import clear_cache clear_cache() - case = extra.CurveTestCase( - name="curve", - knl_class_or_helmholtz_k=k, - curve_fn=curve_fn, - op_type=op_type, - target_order=7, - qbx_order=4, - resolutions=[30]) - - logger.info("\n%s", case) - - # {{{ geometry - - qbx = case.get_layer_potential(actx, case.resolutions[-1], case.target_order) + target_order = 7 + qbx_order = 4 + nelements = 30 + mesh = make_curve_mesh(curve_f, + np.linspace(0, 1, nelements + 1), + target_order) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + pre_density_discr = Discretization(actx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + from pytential.qbx import QBXLayerPotentialSource + qbx = QBXLayerPotentialSource(pre_density_discr, + 4 * target_order, + qbx_order=qbx_order, + # Don't use FMM for now + fmm_order=False) from pytential.qbx.refinement import refine_geometry_collection - places = GeometryCollection(qbx, auto_where=case.name) + places = GeometryCollection(qbx) places = refine_geometry_collection(places, kernel_length_scale=(5 / k if k else None)) - dd = places.auto_source.to_stage1() - density_discr = places.get_discretization(dd.geometry) - - logger.info("nelements: %d", density_discr.mesh.nelements) - logger.info("ndofs: %d", density_discr.ndofs) - - # }}} - - # {{{ symbolic - - sym_u, sym_op = case.get_operator(places.ambient_dim) - bound_op = bind(places, sym_op) + source = places.auto_source.to_stage1() + density_discr = places.get_discretization(source.geometry) - # }}} - - # {{{ dense matrix + op, u_sym, knl_kwargs = _build_op(lpot_id, k=k, + source=places.auto_source, + target=places.auto_target) + bound_op = bind(places, op) from pytential.symbolic.execution import build_matrix - mat = actx.to_numpy( - build_matrix(actx, places, sym_op, sym_u, - context=case.knl_concrete_kwargs)) - - if visualize: - try: - import matplotlib.pyplot as pt - except ImportError: - visualize = False + mat = build_matrix(actx, places, op, u_sym).get() if visualize: from sumpy.tools import build_matrix as build_matrix_via_matvec - mat2 = bound_op.scipy_op(actx, "u", dtype=mat.dtype, - **case.knl_concrete_kwargs) + mat2 = bound_op.scipy_op(actx, "u", dtype=mat.dtype, **knl_kwargs) mat2 = build_matrix_via_matvec(mat2) - - logger.info("real %.5e imag %.5e", - la.norm((mat - mat2).real, "fro") / la.norm(mat2.real, "fro"), - la.norm((mat - mat2).imag, "fro") / la.norm(mat2.imag, "fro")) + print(la.norm((mat - mat2).real, "fro") / la.norm(mat2.real, "fro"), + la.norm((mat - mat2).imag, "fro") / la.norm(mat2.imag, "fro")) pt.subplot(121) pt.imshow(np.log10(np.abs(1.0e-20 + (mat - mat2).real))) @@ -138,7 +254,6 @@ def test_build_matrix(ctx_factory, k, curve_fn, op_type, visualize=False): pt.imshow(np.log10(np.abs(1.0e-20 + (mat - mat2).imag))) pt.colorbar() pt.show() - pt.clf() if visualize: pt.subplot(121) @@ -148,47 +263,36 @@ def test_build_matrix(ctx_factory, k, curve_fn, op_type, visualize=False): pt.imshow(mat.imag) pt.colorbar() pt.show() - pt.clf() - - # }}} - - # {{{ check from pytential.utils import unflatten_from_numpy, flatten_to_numpy - np.random.seed(12) for i in range(5): - if isinstance(sym_u, np.ndarray): + if is_obj_array(u_sym): u = make_obj_array([ np.random.randn(density_discr.ndofs) - for _ in range(len(sym_u)) + for _ in range(len(u_sym)) ]) else: u = np.random.randn(density_discr.ndofs) u_dev = unflatten_from_numpy(actx, density_discr, u) - res_matvec = np.hstack(flatten_to_numpy(actx, - bound_op(actx, u=u_dev, **case.knl_concrete_kwargs) - )) + res_matvec = np.hstack( + flatten_to_numpy(actx, bound_op(actx, u=u_dev)) + ) res_mat = mat.dot(np.hstack(u)) abs_err = la.norm(res_mat - res_matvec, np.inf) rel_err = abs_err / la.norm(res_matvec, np.inf) - logger.info("AbsErr {:.5e} RelErr {:.5e}".format(abs_err, rel_err)) - assert rel_err < 1.0e-13, 'iteration: {}'.format(i) - - # }}} + print("AbsErr {:.5e} RelErr {:.5e}".format(abs_err, rel_err)) + assert rel_err < 1e-13, 'iteration: {}'.format(i) +@pytest.mark.parametrize("factor", [1.0, 0.6]) @pytest.mark.parametrize("ambient_dim", [2, 3]) -@pytest.mark.parametrize("block_builder_type", ["qbx", "p2p"]) -@pytest.mark.parametrize("index_sparsity_factor", [1.0, 0.6]) -@pytest.mark.parametrize("op_type", ["scalar", "scalar_mixed"]) -def test_block_builder(ctx_factory, ambient_dim, - block_builder_type, index_sparsity_factor, op_type, visualize=False): - """Test that block builders and full matrix builders actually match.""" - +@pytest.mark.parametrize("lpot_id", [1, 2]) +def test_p2p_block_builder(ctx_factory, factor, ambient_dim, lpot_id, + visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) @@ -197,85 +301,53 @@ def test_block_builder(ctx_factory, ambient_dim, from sympy.core.cache import clear_cache clear_cache() - if ambient_dim == 2: - case = extra.CurveTestCase( - name="ellipse", - target_order=7, - index_sparsity_factor=index_sparsity_factor, - op_type=op_type, - resolutions=[32], - curve_fn=partial(ellipse, 3.0), - ) - elif ambient_dim == 3: - case = extra.TorusTestCase( - index_sparsity_factor=index_sparsity_factor, - op_type=op_type, - target_order=2, - resolutions=[0], - ) - else: - raise ValueError(f"unsupported dimension: {ambient_dim}") - - logger.info("\n%s", case) - - # {{{ geometry + place_ids = ( + sym.DOFDescriptor( + geometry=sym.DEFAULT_SOURCE, + discr_stage=sym.QBX_SOURCE_STAGE2), + sym.DOFDescriptor(geometry=sym.DEFAULT_TARGET) + ) + target_order = 2 if ambient_dim == 3 else 7 - dd = sym.DOFDescriptor(case.name, discr_stage=sym.QBX_SOURCE_STAGE2) - qbx = case.get_layer_potential(actx, case.resolutions[-1], case.target_order) + places, dofdesc = _build_geometry(actx, + target_order=target_order, + ambient_dim=ambient_dim, + auto_where=place_ids) + op, u_sym, _ = _build_op(lpot_id, + ambient_dim=ambient_dim, + source=places.auto_source, + target=places.auto_target) - places = GeometryCollection(qbx, auto_where=(dd, dd.to_stage1())) + dd = places.auto_source density_discr = places.get_discretization(dd.geometry, dd.discr_stage) - - logger.info("nelements: %d", density_discr.mesh.nelements) - logger.info("ndofs: %d", density_discr.ndofs) - - # }}} - - # {{{ symbolic - - sym_u, sym_op = case.get_operator(ambient_dim) + index_set = _build_block_index(actx, density_discr, factor=factor) + index_set = MatrixBlockIndexRanges(ctx, index_set, index_set) from pytential.symbolic.execution import _prepare_expr - sym_prep_op = _prepare_expr(places, sym_op) - - # }}} + expr = _prepare_expr(places, op) - # {{{ matrix - - index_set = case.get_block_indices(actx, density_discr) - kwargs = dict( - dep_expr=sym_u, + from pytential.symbolic.matrix import P2PMatrixBuilder + mbuilder = P2PMatrixBuilder(actx, + dep_expr=u_sym, other_dep_exprs=[], dep_source=places.get_geometry(dd.geometry), - dep_discr=density_discr, + dep_discr=places.get_discretization(dd.geometry, dd.discr_stage), places=places, - context=case.knl_concrete_kwargs - ) - - if block_builder_type == "qbx": - from pytential.symbolic.matrix import MatrixBuilder - from pytential.symbolic.matrix import \ - NearFieldBlockBuilder as BlockMatrixBuilder - elif block_builder_type == "p2p": - from pytential.symbolic.matrix import P2PMatrixBuilder as MatrixBuilder - from pytential.symbolic.matrix import \ - FarFieldBlockBuilder as BlockMatrixBuilder - kwargs["exclude_self"] = True - else: - raise ValueError(f"unknown block builder type: '{block_builder_type}'") + context={}, + exclude_self=True) + mat = mbuilder(expr) - mat = MatrixBuilder(actx, **kwargs)(sym_prep_op) - blk = BlockMatrixBuilder(actx, index_set=index_set, **kwargs)(sym_prep_op) - - # }}} - - # {{{ check - - if visualize and ambient_dim == 2: - try: - import matplotlib.pyplot as pt - except ImportError: - visualize = False + from pytential.symbolic.matrix import FarFieldBlockBuilder + mbuilder = FarFieldBlockBuilder(actx, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places.get_geometry(dd.geometry), + dep_discr=places.get_discretization(dd.geometry, dd.discr_stage), + places=places, + index_set=index_set, + context={}, + exclude_self=True) + blk = mbuilder(expr) index_set = index_set.get(actx.queue) if visualize and ambient_dim == 2: @@ -291,27 +363,19 @@ def test_block_builder(ctx_factory, ambient_dim, _, (ax1, ax2) = pt.subplots(1, 2, figsize=(10, 8), dpi=300, constrained_layout=True) ax1.imshow(blk_full) - ax1.set_title(type(BlockMatrixBuilder).__name__) + ax1.set_title('FarFieldBlockBuilder') ax2.imshow(mat_full) - ax2.set_title(type(MatrixBuilder).__name__) - - filename = f"matrix_block_{block_builder_type}_{ambient_dim}d" - pt.savefig(filename) + ax2.set_title('P2PMatrixBuilder') + pt.savefig("test_p2p_block_{}d_{:.1f}.png".format(ambient_dim, factor)) - assert extra.max_block_error(mat, blk, index_set) < 1.0e-14 + assert _max_block_error(mat, blk, index_set) < 1.0e-14 - # }}} - - -@pytest.mark.parametrize(('source_discr_stage', 'target_discr_stage'), [ - (sym.QBX_SOURCE_STAGE1, sym.QBX_SOURCE_STAGE1), - (sym.QBX_SOURCE_STAGE2, sym.QBX_SOURCE_STAGE2), - # (sym.QBX_SOURCE_STAGE2, sym.QBX_SOURCE_STAGE1), - ]) -def test_build_matrix_fixed_stage(ctx_factory, - source_discr_stage, target_discr_stage, visualize=False): - """Checks that the block builders match for difference stages.""" +@pytest.mark.parametrize("factor", [1.0, 0.6]) +@pytest.mark.parametrize("ambient_dim", [2, 3]) +@pytest.mark.parametrize("lpot_id", [1, 2]) +def test_qbx_block_builder(ctx_factory, factor, ambient_dim, lpot_id, + visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) @@ -320,91 +384,170 @@ def test_build_matrix_fixed_stage(ctx_factory, from sympy.core.cache import clear_cache clear_cache() - case = extra.CurveTestCase( - name="starfish", - curve_fn=NArmedStarfish(5, 0.25), - - target_order=4, - resolutions=[32], - - index_sparsity_factor=0.6, - op_type="scalar", - tree_kind=None, + place_ids = ( + sym.DOFDescriptor( + geometry=sym.DEFAULT_SOURCE, + discr_stage=sym.QBX_SOURCE_STAGE2), + sym.DOFDescriptor(geometry=sym.DEFAULT_TARGET) ) + target_order = 2 if ambient_dim == 3 else 7 + + places, _ = _build_geometry(actx, + target_order=target_order, + ambient_dim=ambient_dim, + auto_where=place_ids) + op, u_sym, _ = _build_op(lpot_id, + ambient_dim=ambient_dim, + source=places.auto_source, + target=places.auto_target, + qbx_forced_limit="avg") - logger.info("\n%s", case) - - # {{{ geometry - - dd = sym.DOFDescriptor(case.name) - qbx = case.get_layer_potential(actx, case.resolutions[-1], case.target_order) - - places = GeometryCollection({case.name: qbx}, - auto_where=( - dd.copy(discr_stage=source_discr_stage), - dd.copy(discr_stage=target_discr_stage))) + from pytential.symbolic.execution import _prepare_expr + expr = _prepare_expr(places, op) dd = places.auto_source density_discr = places.get_discretization(dd.geometry, dd.discr_stage) + index_set = _build_block_index(actx, density_discr, factor=factor) + index_set = MatrixBlockIndexRanges(ctx, index_set, index_set) - # }}} + from pytential.symbolic.matrix import NearFieldBlockBuilder + mbuilder = NearFieldBlockBuilder(actx, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places.get_geometry(dd.geometry), + dep_discr=places.get_discretization(dd.geometry, dd.discr_stage), + places=places, + index_set=index_set, + context={}) + blk = mbuilder(expr) - # {{{ symbolic + from pytential.symbolic.matrix import MatrixBuilder + mbuilder = MatrixBuilder(actx, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places.get_geometry(dd.geometry), + dep_discr=places.get_discretization(dd.geometry, dd.discr_stage), + places=places, + context={}) + mat = mbuilder(expr) - if source_discr_stage is target_discr_stage: - qbx_forced_limit = -1 - else: - qbx_forced_limit = None + index_set = index_set.get(queue) + if visualize: + blk_full = np.zeros_like(mat) + mat_full = np.zeros_like(mat) - sym_u, sym_op = case.get_operator(places.ambient_dim, qbx_forced_limit) + for i in range(index_set.nblocks): + itgt, isrc = index_set.block_indices(i) - from pytential.symbolic.execution import _prepare_expr - sym_prep_op = _prepare_expr(places, sym_op) + blk_full[np.ix_(itgt, isrc)] = index_set.block_take(blk, i) + mat_full[np.ix_(itgt, isrc)] = index_set.take(mat, i) - # }}} + _, (ax1, ax2) = pt.subplots(1, 2, + figsize=(10, 8), constrained_layout=True) + ax1.imshow(mat_full) + ax1.set_title('MatrixBuilder') + ax2.imshow(blk_full) + ax2.set_title('NearFieldBlockBuilder') + pt.savefig("test_qbx_block_builder.png", dpi=300) - # {{{ check + assert _max_block_error(mat, blk, index_set) < 1.0e-14 - source_discr = places.get_discretization(case.name, source_discr_stage) - target_discr = places.get_discretization(case.name, target_discr_stage) - logger.info("nelements: %d", density_discr.mesh.nelements) - logger.info("ndofs: %d", source_discr.ndofs) - logger.info("ndofs: %d", target_discr.ndofs) +@pytest.mark.parametrize(('source_discr_stage', 'target_discr_stage'), + [(sym.QBX_SOURCE_STAGE1, sym.QBX_SOURCE_STAGE1), + (sym.QBX_SOURCE_STAGE2, sym.QBX_SOURCE_STAGE2)]) +def test_build_matrix_places(ctx_factory, + source_discr_stage, target_discr_stage, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) - icols = case.get_block_indices(actx, source_discr, matrix_indices=False) - irows = case.get_block_indices(actx, target_discr, matrix_indices=False) - index_set = MatrixBlockIndexRanges(actx.context, icols, irows) + # prevent cache explosion + from sympy.core.cache import clear_cache + clear_cache() - kwargs = dict( - dep_expr=sym_u, - other_dep_exprs=[], - dep_source=places.get_geometry(case.name), - dep_discr=density_discr, - places=places, - context=case.knl_concrete_kwargs, + qbx_forced_limit = -1 + place_ids = ( + sym.DOFDescriptor( + geometry=sym.DEFAULT_SOURCE, + discr_stage=source_discr_stage), + sym.DOFDescriptor(geometry=sym.DEFAULT_TARGET) ) - # qbx - from pytential.symbolic import matrix - mat = matrix.MatrixBuilder( - actx, **kwargs)(sym_prep_op) - blk = matrix.NearFieldBlockBuilder( - actx, index_set=index_set, **kwargs)(sym_prep_op) + # build test operators + places, _ = _build_geometry(actx, + nelements=8, + target_order=2, + ambient_dim=2, + curve_f=partial(ellipse, 1.0), + auto_where=place_ids) + op, u_sym, _ = _build_op(lpot_id=1, + ambient_dim=2, + source=places.auto_source, + target=places.auto_target, + qbx_forced_limit=qbx_forced_limit) + + dd = places.auto_target + target_discr = places.get_discretization(dd.geometry, dd.discr_stage) + dd = places.auto_source + source_discr = places.get_discretization(dd.geometry, dd.discr_stage) - assert mat.shape == (target_discr.ndofs, source_discr.ndofs) - assert extra.max_block_error(mat, blk, index_set.get(queue)) < 1.0e-14 + index_set = _build_block_index(actx, source_discr, factor=0.6) + index_set = MatrixBlockIndexRanges(ctx, index_set, index_set) - # p2p - mat = matrix.P2PMatrixBuilder( - actx, exclude_self=True, **kwargs)(sym_prep_op) - blk = matrix.FarFieldBlockBuilder( - actx, index_set=index_set, exclude_self=True, **kwargs)(sym_prep_op) + from pytential.symbolic.execution import _prepare_expr + op = _prepare_expr(places, op) + + # build full QBX matrix + from pytential.symbolic.matrix import MatrixBuilder + mbuilder = MatrixBuilder(actx, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places.get_geometry(dd.geometry), + dep_discr=places.get_discretization(dd.geometry, dd.discr_stage), + places=places, + context={}) + qbx_mat = mbuilder(op) - assert mat.shape == (target_discr.ndofs, source_discr.ndofs) - assert extra.max_block_error(mat, blk, index_set.get(queue)) < 1.0e-14 + # build full p2p matrix + from pytential.symbolic.matrix import P2PMatrixBuilder + mbuilder = P2PMatrixBuilder(actx, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places.get_geometry(dd.geometry), + dep_discr=places.get_discretization(dd.geometry, dd.discr_stage), + places=places, + context={}) + p2p_mat = mbuilder(op) - # }}} + assert p2p_mat.shape == (target_discr.ndofs, source_discr.ndofs) + + # build block qbx and p2p matrices + from pytential.symbolic.matrix import NearFieldBlockBuilder + mbuilder = NearFieldBlockBuilder(actx, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places.get_geometry(dd.geometry), + dep_discr=places.get_discretization(dd.geometry, dd.discr_stage), + places=places, + index_set=index_set, + context={}) + mat = mbuilder(op) + if dd.discr_stage is not None: + assert _max_block_error(qbx_mat, mat, index_set.get(queue)) < 1.0e-14 + + from pytential.symbolic.matrix import FarFieldBlockBuilder + mbuilder = FarFieldBlockBuilder(actx, + dep_expr=u_sym, + other_dep_exprs=[], + dep_source=places.get_geometry(dd.geometry), + dep_discr=places.get_discretization(dd.geometry, dd.discr_stage), + places=places, + index_set=index_set, + context={}, + exclude_self=True) + mat = mbuilder(op) + assert _max_block_error(p2p_mat, mat, index_set.get(queue)) < 1.0e-14 if __name__ == "__main__": -- GitLab