diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py index 0f79148e8961433de65fc960487c5e4cd6202e94..4183cd0a268fbba2807d9f1baf4d5fdf3f67fa0a 100644 --- a/pytential/linalg/proxy.py +++ b/pytential/linalg/proxy.py @@ -39,9 +39,7 @@ Proxy Point Generation ~~~~~~~~~~~~~~~~~~~~~~ .. autoclass:: ProxyGenerator - .. autofunction:: partition_by_nodes -.. autofunction:: partition_from_coarse .. autofunction:: gather_block_neighbor_points .. autofunction:: gather_block_interaction_points @@ -50,16 +48,14 @@ Proxy Point Generation # {{{ point index partitioning -def partition_by_nodes(actx, discr, use_tree=True, max_nodes_in_box=None): +def partition_by_nodes(actx, discr, + tree_kind="adaptive-level-restricted", 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 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 tree_kind: if not *None*, it is passed to :class:`boxtree.TreeBuilder`. :arg max_nodes_in_box: passed to :class:`boxtree.TreeBuilder`. :return: a :class:`sumpy.tools.BlockIndexRanges`. @@ -69,7 +65,7 @@ def partition_by_nodes(actx, discr, use_tree=True, max_nodes_in_box=None): # FIXME: this is just an arbitrary value max_nodes_in_box = 32 - if use_tree: + if tree_kind is not None: from boxtree import box_flags_enum from boxtree import TreeBuilder @@ -78,7 +74,8 @@ def partition_by_nodes(actx, discr, use_tree=True, max_nodes_in_box=None): from meshmode.dof_array import flatten, thaw tree, _ = builder(actx.queue, flatten(thaw(actx, discr.nodes())), - max_particles_in_box=max_nodes_in_box) + max_particles_in_box=max_nodes_in_box, + kind=tree_kind) tree = tree.get(actx.queue) leaf_boxes, = (tree.box_flags @@ -99,7 +96,7 @@ def partition_by_nodes(actx, discr, use_tree=True, max_nodes_in_box=None): ranges = actx.from_numpy(np.arange( 0, discr.ndofs + 1, - discr.ndofs // max_nodes_in_box, dtype=np.int)) + 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 3942adba7d81e06372aa7e54cf8c1cc50a49cf76..dbc807a760212432bb879fc3b39c4b9122a75aa0 100644 --- a/test/extra_int_eq_data.py +++ b/test/extra_int_eq_data.py @@ -357,6 +357,29 @@ 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 new file mode 100644 index 0000000000000000000000000000000000000000..983c815ec503eb735941f50875d7785280b30c73 --- /dev/null +++ b/test/extra_matrix_data.py @@ -0,0 +1,112 @@ +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 8e485d6e9fb3e3cb8fc6f541971796bc8fcf4f63..09a1367bc01d1c970dc62f39149ae968d740c86d 100644 --- a/test/test_linalg_proxy.py +++ b/test/test_linalg_proxy.py @@ -22,41 +22,46 @@ 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 ( # noqa - ellipse, NArmedStarfish, generate_torus, make_curve_mesh) +from meshmode.mesh.generation import ellipse 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("method", "unknown"), - "tree" if kwargs.get("use_tree", False) else "linear", - kwargs.get("pid", "stage1"), + kwargs.get("tree_kind", "linear").replace("-", "_"), + kwargs.get("discr_stage", "stage1"), discr.ambient_dim ] pt.figure(figsize=(10, 8), dpi=300) pt.plot(np.diff(indices.ranges)) - pt.savefig("test_partition_{0}_{1}_{3}d_ranges_{2}.png".format(*args)) + pt.savefig("test_partition_{1}_{3}d_ranges_{2}.png".format(*args)) pt.clf() from pytential.utils import flatten_to_numpy @@ -64,24 +69,18 @@ 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_{0}_{1}_{3}d_{2}.png".format(*args)) + pt.savefig("test_partition_{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) @@ -100,38 +99,70 @@ def _plot_partition_indices(actx, discr, indices, **kwargs): ]) -@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): +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).""" + ctx = ctx_factory() queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) - 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) + 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) if visualize: - _plot_partition_indices(actx, discr, indices, use_tree=use_tree) + plot_partition_indices(actx, density_discr, indices, tree_kind=tree_kind) + # }}} + + +@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) - places, dofdesc = _build_geometry(actx, ambient_dim=ambient_dim) - dofdesc = dofdesc.to_stage1() + case = case.copy(index_sparsity_factor=index_sparsity_factor) + logger.info("\n%s", case) + + # {{{ generate proxies - density_discr = places.get_discretization(dofdesc.geometry, dofdesc.discr_stage) - srcindices = _build_block_index(actx, density_discr, factor=factor) + 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) from pytential.linalg.proxy import ProxyGenerator generator = ProxyGenerator(places) proxies, pxyranges, pxycenters, pxyradii = \ - generator(actx, dofdesc, srcindices) + generator(actx, places.auto_source, srcindices) proxies = np.vstack([actx.to_numpy(p) for p in proxies]) pxyranges = actx.to_numpy(pxyranges) @@ -144,8 +175,13 @@ def test_proxy_generator(ctx_factory, ambient_dim, factor, visualize=False): 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 @@ -212,32 +248,45 @@ def test_proxy_generator(ctx_factory, ambient_dim, factor, visualize=False): 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) - places, dofdesc = _build_geometry(actx, ambient_dim=ambient_dim) - dofdesc = dofdesc.to_stage1() + case = case.copy(index_sparsity_factor=index_sparsity_factor) + logger.info("\n%s", case) + + # {{{ check neighboring points - density_discr = places.get_discretization(dofdesc.geometry, dofdesc.discr_stage) - srcindices = _build_block_index(actx, density_discr, factor=factor) + 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) # generate proxy points from pytential.linalg.proxy import ProxyGenerator generator = ProxyGenerator(places) - _, _, pxycenters, pxyradii = generator(actx, dofdesc, srcindices) + _, _, pxycenters, pxyradii = generator(actx, places.auto_source, 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, dofdesc, srcindices) + places, places.auto_source, srcindices) srcindices = srcindices.get(queue) nbrindices = nbrindices.get(queue) @@ -248,8 +297,13 @@ def test_interaction_points(ctx_factory, ambient_dim, factor, visualize=False): 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()) @@ -300,6 +354,7 @@ def test_interaction_points(ctx_factory, ambient_dim, factor, visualize=False): vis.write_vtk_file(filename, [ ("marker", marker_dev), ]) + # }}} if __name__ == "__main__": diff --git a/test/test_matrix.py b/test/test_matrix.py index ce429642df74870fa2553fffac0e548ee5bb4061..291cefc784ae67e2474294760f17efb7a5d71f2d 100644 --- a/test/test_matrix.py +++ b/test/test_matrix.py @@ -31,171 +31,40 @@ 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 ( # noqa - ellipse, NArmedStarfish, make_curve_mesh, generate_torus) +from meshmode.mesh.generation import ellipse, NArmedStarfish import pytest from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) -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 +import extra_matrix_data as extra +import logging +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO) @pytest.mark.skipif(USE_SYMENGINE, reason="https://gitlab.tiker.net/inducer/sumpy/issues/25") @pytest.mark.parametrize("k", [0, 42]) -@pytest.mark.parametrize("curve_f", [ +@pytest.mark.parametrize("curve_fn", [ partial(ellipse, 3), NArmedStarfish(5, 0.25)]) -@pytest.mark.parametrize("lpot_id", [2, 3]) -def test_matrix_build(ctx_factory, k, curve_f, lpot_id, visualize=False): +@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. + """ + cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) @@ -204,48 +73,63 @@ def test_matrix_build(ctx_factory, k, curve_f, lpot_id, visualize=False): from sympy.core.cache import clear_cache clear_cache() - 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) + 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) from pytential.qbx.refinement import refine_geometry_collection - places = GeometryCollection(qbx) + places = GeometryCollection(qbx, auto_where=case.name) places = refine_geometry_collection(places, kernel_length_scale=(5 / k if k else None)) - source = places.auto_source.to_stage1() - density_discr = places.get_discretization(source.geometry) + 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) - op, u_sym, knl_kwargs = _build_op(lpot_id, k=k, - source=places.auto_source, - target=places.auto_target) - bound_op = bind(places, op) + # }}} + + # {{{ dense matrix from pytential.symbolic.execution import build_matrix - mat = build_matrix(actx, places, op, u_sym).get() + 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 if visualize: from sumpy.tools import build_matrix as build_matrix_via_matvec - mat2 = bound_op.scipy_op(actx, "u", dtype=mat.dtype, **knl_kwargs) + mat2 = bound_op.scipy_op(actx, "u", dtype=mat.dtype, + **case.knl_concrete_kwargs) mat2 = build_matrix_via_matvec(mat2) - print(la.norm((mat - mat2).real, "fro") / la.norm(mat2.real, "fro"), - la.norm((mat - mat2).imag, "fro") / la.norm(mat2.imag, "fro")) + + 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")) pt.subplot(121) pt.imshow(np.log10(np.abs(1.0e-20 + (mat - mat2).real))) @@ -254,6 +138,7 @@ def test_matrix_build(ctx_factory, k, curve_f, lpot_id, 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) @@ -263,36 +148,47 @@ def test_matrix_build(ctx_factory, k, curve_f, lpot_id, 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 is_obj_array(u_sym): + if isinstance(sym_u, np.ndarray): u = make_obj_array([ np.random.randn(density_discr.ndofs) - for _ in range(len(u_sym)) + for _ in range(len(sym_u)) ]) 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)) - ) + res_matvec = np.hstack(flatten_to_numpy(actx, + bound_op(actx, u=u_dev, **case.knl_concrete_kwargs) + )) 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) - print("AbsErr {:.5e} RelErr {:.5e}".format(abs_err, rel_err)) - assert rel_err < 1e-13, 'iteration: {}'.format(i) + logger.info("AbsErr {:.5e} RelErr {:.5e}".format(abs_err, rel_err)) + assert rel_err < 1.0e-13, 'iteration: {}'.format(i) + + # }}} -@pytest.mark.parametrize("factor", [1.0, 0.6]) @pytest.mark.parametrize("ambient_dim", [2, 3]) -@pytest.mark.parametrize("lpot_id", [1, 2]) -def test_p2p_block_builder(ctx_factory, factor, ambient_dim, lpot_id, - visualize=False): +@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.""" + ctx = ctx_factory() queue = cl.CommandQueue(ctx) actx = PyOpenCLArrayContext(queue) @@ -301,53 +197,85 @@ def test_p2p_block_builder(ctx_factory, factor, ambient_dim, lpot_id, from sympy.core.cache import clear_cache clear_cache() - 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 + 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}") - 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) + logger.info("\n%s", case) - dd = places.auto_source + # {{{ geometry + + dd = sym.DOFDescriptor(case.name, discr_stage=sym.QBX_SOURCE_STAGE2) + qbx = case.get_layer_potential(actx, case.resolutions[-1], case.target_order) + + places = GeometryCollection(qbx, auto_where=(dd, dd.to_stage1())) 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) + + 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) from pytential.symbolic.execution import _prepare_expr - expr = _prepare_expr(places, op) + sym_prep_op = _prepare_expr(places, sym_op) - 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={}, - exclude_self=True) - mat = mbuilder(expr) + # }}} - from pytential.symbolic.matrix import FarFieldBlockBuilder - mbuilder = FarFieldBlockBuilder(actx, - dep_expr=u_sym, + # {{{ matrix + + index_set = case.get_block_indices(actx, density_discr) + kwargs = dict( + dep_expr=sym_u, other_dep_exprs=[], dep_source=places.get_geometry(dd.geometry), - dep_discr=places.get_discretization(dd.geometry, dd.discr_stage), + dep_discr=density_discr, places=places, - index_set=index_set, - context={}, - exclude_self=True) - blk = mbuilder(expr) + 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}'") + + 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 index_set = index_set.get(actx.queue) if visualize and ambient_dim == 2: @@ -363,19 +291,27 @@ def test_p2p_block_builder(ctx_factory, factor, ambient_dim, lpot_id, _, (ax1, ax2) = pt.subplots(1, 2, figsize=(10, 8), dpi=300, constrained_layout=True) ax1.imshow(blk_full) - ax1.set_title('FarFieldBlockBuilder') + ax1.set_title(type(BlockMatrixBuilder).__name__) ax2.imshow(mat_full) - ax2.set_title('P2PMatrixBuilder') - pt.savefig("test_p2p_block_{}d_{:.1f}.png".format(ambient_dim, factor)) + ax2.set_title(type(MatrixBuilder).__name__) - assert _max_block_error(mat, blk, index_set) < 1.0e-14 + filename = f"matrix_block_{block_builder_type}_{ambient_dim}d" + pt.savefig(filename) + assert extra.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) @@ -384,170 +320,91 @@ def test_qbx_block_builder(ctx_factory, factor, ambient_dim, lpot_id, from sympy.core.cache import clear_cache clear_cache() - 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") + case = extra.CurveTestCase( + name="starfish", + curve_fn=NArmedStarfish(5, 0.25), - from pytential.symbolic.execution import _prepare_expr - expr = _prepare_expr(places, op) + target_order=4, + resolutions=[32], - 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) + index_sparsity_factor=0.6, + op_type="scalar", + tree_kind=None, + ) - 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) + logger.info("\n%s", case) - 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) + # {{{ geometry - index_set = index_set.get(queue) - if visualize: - blk_full = np.zeros_like(mat) - mat_full = np.zeros_like(mat) + dd = sym.DOFDescriptor(case.name) + qbx = case.get_layer_potential(actx, case.resolutions[-1], case.target_order) - for i in range(index_set.nblocks): - itgt, isrc = index_set.block_indices(i) + places = GeometryCollection({case.name: qbx}, + auto_where=( + dd.copy(discr_stage=source_discr_stage), + dd.copy(discr_stage=target_discr_stage))) - blk_full[np.ix_(itgt, isrc)] = index_set.block_take(blk, i) - mat_full[np.ix_(itgt, isrc)] = index_set.take(mat, i) + dd = places.auto_source + density_discr = places.get_discretization(dd.geometry, dd.discr_stage) - _, (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) + # }}} - assert _max_block_error(mat, blk, index_set) < 1.0e-14 + # {{{ symbolic + if source_discr_stage is target_discr_stage: + qbx_forced_limit = -1 + else: + qbx_forced_limit = None -@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) + sym_u, sym_op = case.get_operator(places.ambient_dim, qbx_forced_limit) - # prevent cache explosion - from sympy.core.cache import clear_cache - clear_cache() + from pytential.symbolic.execution import _prepare_expr + sym_prep_op = _prepare_expr(places, sym_op) - qbx_forced_limit = -1 - place_ids = ( - sym.DOFDescriptor( - geometry=sym.DEFAULT_SOURCE, - discr_stage=source_discr_stage), - sym.DOFDescriptor(geometry=sym.DEFAULT_TARGET) - ) + # }}} - # 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) + # {{{ check - index_set = _build_block_index(actx, source_discr, factor=0.6) - index_set = MatrixBlockIndexRanges(ctx, index_set, index_set) + source_discr = places.get_discretization(case.name, source_discr_stage) + target_discr = places.get_discretization(case.name, target_discr_stage) - from pytential.symbolic.execution import _prepare_expr - op = _prepare_expr(places, op) + logger.info("nelements: %d", density_discr.mesh.nelements) + logger.info("ndofs: %d", source_discr.ndofs) + logger.info("ndofs: %d", target_discr.ndofs) - # 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) + 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) - # build full p2p matrix - from pytential.symbolic.matrix import P2PMatrixBuilder - mbuilder = P2PMatrixBuilder(actx, - dep_expr=u_sym, + kwargs = dict( + dep_expr=sym_u, other_dep_exprs=[], - dep_source=places.get_geometry(dd.geometry), - dep_discr=places.get_discretization(dd.geometry, dd.discr_stage), + dep_source=places.get_geometry(case.name), + dep_discr=density_discr, places=places, - context={}) - p2p_mat = mbuilder(op) + context=case.knl_concrete_kwargs, + ) - assert p2p_mat.shape == (target_discr.ndofs, source_discr.ndofs) + # 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 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 + assert mat.shape == (target_discr.ndofs, source_discr.ndofs) + assert extra.max_block_error(mat, blk, index_set.get(queue)) < 1.0e-14 + + # 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) + + assert mat.shape == (target_discr.ndofs, source_discr.ndofs) + assert extra.max_block_error(mat, blk, index_set.get(queue)) < 1.0e-14 + + # }}} if __name__ == "__main__":