From 1c7ca54f5b2721dc0a92677bca65ed1934e81bca Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 13 Feb 2022 23:30:38 -0600 Subject: [PATCH 1/5] Support new boxtree tree build --- .gitignore | 2 + requirements.txt | 3 +- ...test_nearfield_interaction_completeness.py | 28 ++++---------- volumential/expansion_wrangler_fpnd.py | 34 ++++++++--------- volumential/expansion_wrangler_interface.py | 6 +-- volumential/meshgen.py | 38 ++++++++----------- volumential/nearfield_potential_table.py | 9 +++-- volumential/volume_fmm.py | 2 +- 8 files changed, 52 insertions(+), 70 deletions(-) diff --git a/.gitignore b/.gitignore index 00ac750..80e0de7 100644 --- a/.gitignore +++ b/.gitignore @@ -39,3 +39,5 @@ gmsh_tmp/ .coverage pytest.xml prof/ + +.ccls-cache/ diff --git a/requirements.txt b/requirements.txt index 099dcc8..3c859ba 100644 --- a/requirements.txt +++ b/requirements.txt @@ -10,7 +10,8 @@ filelock -e git+https://gitlab.tiker.net/inducer/cgen.git#egg=cgen -e git+https://gitlab.tiker.net/inducer/loopy.git#egg=loopy -e git+https://gitlab.tiker.net/inducer/pymbolic.git#egg=pymbolic --e git+https://gitlab.tiker.net/inducer/boxtree.git#egg=boxtree +# -e git+https://gitlab.tiker.net/inducer/boxtree.git#egg=boxtree +-e git+https://gitlab.tiker.net/xywei/boxtree.git@tree-build-refactor#egg=boxtree -e git+https://gitlab.tiker.net/inducer/arraycontext.git#egg=arraycontext -e git+https://gitlab.tiker.net/inducer/meshmode.git#egg=meshmode -e git+https://gitlab.tiker.net/inducer/genpy.git#egg=genpy diff --git a/test/test_nearfield_interaction_completeness.py b/test/test_nearfield_interaction_completeness.py index cd590c7..6dd461f 100644 --- a/test/test_nearfield_interaction_completeness.py +++ b/test/test_nearfield_interaction_completeness.py @@ -41,41 +41,29 @@ def drive_test_completeness(ctx, queue, dim, q_order): def source_field(x): assert len(x) == dim - return 1 + return np.ones_like(x[0]) # {{{ generate quad points import volumential.meshgen as mg - q_points, q_weights, q_radii = mg.make_uniform_cubic_grid( - degree=q_order, level=n_levels, dim=dim - ) - - assert len(q_points) == len(q_weights) - assert q_points.shape[1] == dim + q_points, q_weights, _ = mg.make_uniform_cubic_grid( + nqpoints=q_order+1, level=n_levels, dim=dim) - q_points_org = q_points - q_points = np.ascontiguousarray(np.transpose(q_points)) + assert len(q_points) == dim + len(q_weights) == len(q_points[0]) from pytools.obj_array import make_obj_array + q_points_org = q_points q_points = make_obj_array( - [cl.array.to_device(queue, q_points[i]) for i in range(dim)] - ) + [cl.array.to_device(queue, q_points[i]) for i in range(dim)]) q_weights = cl.array.to_device(queue, q_weights) - if q_radii is not None: - q_radii = cl.array.to_device(queue, q_radii) # }}} - # {{{ discretize the source field - - source_vals = cl.array.to_device( - queue, np.array([source_field(qp) for qp in q_points_org]) - ) - - # }}} End discretize the source field + source_vals = cl.array.to_device(queue, source_field(q_points_org)) # {{{ build tree and traversals diff --git a/volumential/expansion_wrangler_fpnd.py b/volumential/expansion_wrangler_fpnd.py index b691f63..4c0cd48 100644 --- a/volumential/expansion_wrangler_fpnd.py +++ b/volumential/expansion_wrangler_fpnd.py @@ -31,9 +31,9 @@ from pytools.obj_array import make_obj_array # from pytools import memoize_method from volumential.nearfield_potential_table import NearFieldInteractionTable from volumential.expansion_wrangler_interface import ( - ExpansionWranglerInterface, ExpansionWranglerCodeContainerInterface) -from sumpy.fmm import SumpyExpansionWrangler, \ - SumpyTimingFuture, SumpyExpansionWranglerCodeContainer + ExpansionWranglerInterface, TreeIndependentDataForWranglerInterface) +from sumpy.fmm import (SumpyExpansionWrangler, + SumpyTimingFuture, SumpyTreeIndependentDataForWrangler) from boxtree.pyfmmlib_integration import FMMLibExpansionWrangler from sumpy.kernel import ( @@ -68,9 +68,9 @@ def inverse_id_map(queue, mapped_ids): # {{{ sumpy backend -class FPNDSumpyExpansionWranglerCodeContainer( - ExpansionWranglerCodeContainerInterface, - SumpyExpansionWranglerCodeContainer): +class FPNDSumpyTreeIndependentDataForWrangler( + TreeIndependentDataForWranglerInterface, + SumpyTreeIndependentDataForWrangler): """Objects of this type serve as a place to keep the code needed for ExpansionWrangler if it is using sumpy to perform multipole expansion and manipulations. @@ -80,7 +80,10 @@ class FPNDSumpyExpansionWranglerCodeContainer( more ephemeral than the code, the code's lifetime is decoupled by storing it in this object. """ - get_wrangler = SumpyExpansionWranglerCodeContainer.get_wrangler + + @property + def wrangler_cls(self): + return SumpyExpansionWrangler class FPNDSumpyExpansionWrangler( @@ -579,9 +582,8 @@ class FPNDSumpyExpansionWrangler( # {{{ fmmlib backend (for laplace, helmholtz) -class FPNDFMMLibExpansionWranglerCodeContainer( - ExpansionWranglerCodeContainerInterface, - ): +class FPNDFMMLibTreeIndependentDataForWrangler( + TreeIndependentDataForWranglerInterface): """Objects of this type serve as a place to keep the code needed for ExpansionWrangler if it is using fmmlib to perform multipole expansion and manipulations. @@ -600,13 +602,9 @@ class FPNDFMMLibExpansionWranglerCodeContainer( self.target_kernels = target_kernels self.exclude_self = True - def get_wrangler(self, queue, tree, dtype, fmm_level_to_order, - source_extra_kwargs={}, kernel_extra_kwargs=None, - *args, **kwargs): - return FPNDFMMLibExpansionWrangler(self, queue, tree, - dtype, fmm_level_to_order, - source_extra_kwargs, kernel_extra_kwargs, - *args, **kwargs) + @property + def wrangler_cls(self): + return FPNDFMMLibExpansionWrangler class FPNDFMMLibExpansionWrangler( @@ -1201,7 +1199,7 @@ class FPNDFMMLibExpansionWrangler( # }}} End fmmlib backend (for laplace, helmholtz) -class FPNDExpansionWranglerCodeContainer(FPNDSumpyExpansionWranglerCodeContainer): +class FPNDTreeIndependentDataForWrangler(FPNDSumpyTreeIndependentDataForWrangler): """The default code container. """ diff --git a/volumential/expansion_wrangler_interface.py b/volumential/expansion_wrangler_interface.py index 4418f13..5b09290 100644 --- a/volumential/expansion_wrangler_interface.py +++ b/volumential/expansion_wrangler_interface.py @@ -213,10 +213,10 @@ class ExpansionWranglerInterface(object): # }}} End expansion wrangler interface -# {{{ code container interface +# {{{ tree-independent data for wrangler interface -class ExpansionWranglerCodeContainerInterface(object): +class TreeIndependentDataForWranglerInterface(object): """ Abstract expansion code container interface. The interface is adapted from, and stays compatible with boxtree/fmm. @@ -230,6 +230,6 @@ class ExpansionWranglerCodeContainerInterface(object): """ pass -# }}} End code container interface +# }}} # vim: filetype=pyopencl:fdm=marker diff --git a/volumential/meshgen.py b/volumential/meshgen.py index cc936b2..6c616f6 100644 --- a/volumential/meshgen.py +++ b/volumential/meshgen.py @@ -182,8 +182,9 @@ except ImportError as e: try: logger.info("Trying out BoxTree.TreeInteractiveBuild interface.") - import boxtree.tree_interactive_build # noqa: F401 - from modepy import LegendreGaussQuadrature + from boxtree.tree_build import ( + make_tob_root, uniformly_refined, distribute_quadrature_rule) + import modepy as mp provider = "meshgen_boxtree" @@ -195,38 +196,29 @@ except ImportError as e: else: # {{{ Meshgen via BoxTree logger.info("Using Meshgen via BoxTree interface.") - from boxtree.tree_interactive_build import BoxTree, QuadratureOnBoxTree def greet(): return "Hello from Meshgen via BoxTree!" - def make_uniform_cubic_grid(degree, nlevels=1, dim=2, queue=None, **kwargs): + def make_uniform_cubic_grid(nqpoints, nlevels=1, dim=2, actx=None, **kwargs): """Uniform cubic grid in [-1,1]^dim. This function provides backward compatibility with meshgen_dealii. """ - if queue is None: - ctx = cl.create_some_context() - queue = cl.CommandQueue(ctx) - # For meshgen_dealii compatibility if "level" in kwargs: nlevels = kwargs["level"] - tree = BoxTree() - tree.generate_uniform_boxtree( - queue, nlevels=nlevels, root_extent=2, root_vertex=np.zeros(dim) - 1 - ) - quad_rule = LegendreGaussQuadrature(degree - 1) - quad = QuadratureOnBoxTree(tree, quad_rule) - q_weights = quad.get_q_weights(queue).get(queue) - q_points_dev = quad.get_q_points(queue) - n_all_q_points = len(q_weights) - q_points = np.zeros((n_all_q_points, dim)) - for d in range(dim): - q_points[:, d] = q_points_dev[d].get(queue) - - # Adding a placeholder for deprecated point radii - return (q_points, q_weights, None) + lower_bounds = [-1, ] * dim + upper_bounds = [1, ] * dim + tob = make_tob_root(dim=dim, bbox=[lower_bounds, upper_bounds]) + for _ in range(nlevels - 1): + tob = uniformly_refined(tob) + + degree = nqpoints - 1 + quad = mp.LegendreGaussQuadrature(degree) + x, q = distribute_quadrature_rule(tob, quad) + radii = None + return (x, q, radii) class MeshGen1D(MeshGenBase): """Meshgen in 1D diff --git a/volumential/nearfield_potential_table.py b/volumential/nearfield_potential_table.py index 8cef4d2..5e5c9bb 100644 --- a/volumential/nearfield_potential_table.py +++ b/volumential/nearfield_potential_table.py @@ -304,14 +304,14 @@ class NearFieldInteractionTable(object): queue = None q_points, _, _ = mg.make_uniform_cubic_grid( - degree=quad_order, level=1, dim=self.dim, + nqpoints=quad_order-1, level=1, dim=self.dim, queue=queue) # map to source box mapped_q_points = np.array( [ 0.5 * self.source_box_extent * (qp + np.ones(self.dim)) - for qp in q_points + for qp in q_points.T ] ) # sort in dictionary order, preserve only the leading @@ -401,10 +401,11 @@ class NearFieldInteractionTable(object): idx = self.unwrap_mode_index(mode_index) xi = ( - np.array([p[self.dim - 1] for p in self.q_points[: self.quad_order]]) + np.array([p[self.dim - 1] for p in self.q_points[: self.quad_order + 1]]) / self.source_box_extent ) - assert len(xi) == self.quad_order + assert len(xi) == self.quad_order - 1 + print(xi) yi = [] for d in range(self.dim): diff --git a/volumential/volume_fmm.py b/volumential/volume_fmm.py index 0b1a451..8848b22 100644 --- a/volumential/volume_fmm.py +++ b/volumential/volume_fmm.py @@ -30,7 +30,7 @@ __doc__ = """ import numpy as np import pyopencl as cl from pytools.obj_array import make_obj_array -from boxtree.fmm import TimingRecorder +from boxtree.timing import TimingRecorder from volumential.expansion_wrangler_interface import ExpansionWranglerInterface from volumential.expansion_wrangler_fpnd import ( FPNDSumpyExpansionWrangler, FPNDFMMLibExpansionWrangler) -- GitLab From 5afcff8e1df3aa4acb0d53ccc55a86d1be808fb1 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 23 Feb 2022 18:07:25 -0600 Subject: [PATCH 2/5] Boxtree based meshgen, and some compatibility fixes (more to come) --- examples/laplace2d.py | 13 +- test/test_expansion_wrangler_interface.py | 123 --------- ...test_nearfield_interaction_completeness.py | 7 +- volumential/droste.py | 4 +- volumential/expansion_wrangler_fpnd.py | 163 +++--------- volumential/expansion_wrangler_interface.py | 235 ------------------ volumential/meshgen.py | 186 ++++++-------- volumential/nearfield_potential_table.py | 9 +- volumential/volume_fmm.py | 2 - 9 files changed, 131 insertions(+), 611 deletions(-) delete mode 100644 test/test_expansion_wrangler_interface.py delete mode 100644 volumential/expansion_wrangler_interface.py diff --git a/examples/laplace2d.py b/examples/laplace2d.py index 383951f..07a9271 100644 --- a/examples/laplace2d.py +++ b/examples/laplace2d.py @@ -280,11 +280,9 @@ def main(): exclude_self = True from volumential.expansion_wrangler_fpnd import ( - FPNDExpansionWranglerCodeContainer, - FPNDExpansionWrangler - ) + FPNDTreeIndependentDataForWrangler, FPNDExpansionWrangler) - wcc = FPNDExpansionWranglerCodeContainer( + tree_indep = FPNDTreeIndependentDataForWrangler( ctx, partial(mpole_expn_class, knl), partial(local_expn_class, knl), @@ -299,12 +297,9 @@ def main(): self_extra_kwargs = {} wrangler = FPNDExpansionWrangler( - code_container=wcc, - queue=queue, - tree=tree, - near_field_table=nftable, - dtype=dtype, + tree_indep, trav, dtype, fmm_level_to_order=lambda kernel, kernel_args, tree, lev: m_order, + near_field_table=nftable, quad_order=q_order, self_extra_kwargs=self_extra_kwargs, ) diff --git a/test/test_expansion_wrangler_interface.py b/test/test_expansion_wrangler_interface.py deleted file mode 100644 index 11ade59..0000000 --- a/test/test_expansion_wrangler_interface.py +++ /dev/null @@ -1,123 +0,0 @@ -from __future__ import absolute_import, division, print_function - -__copyright__ = "Copyright (C) 2017 - 2018 Xiaoyu Wei" - -__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. -""" - -import volumential.expansion_wrangler_interface as ewi - - -def test_interface_allocation(): - # make the interface instantiable by overriding its set of abstract methods - ewi.ExpansionWranglerInterface.__abstractmethods__ = frozenset() - wrglr = ewi.ExpansionWranglerInterface() - try: - wrglr.multipole_expansion_zeros() - wrglr.local_expansion_zeros() - wrglr.output_zeros() - except NotImplementedError: - assert False - assert True - - -def test_interface_utilities(): - # make the interface instantiable by overriding its set of abstract methods - ewi.ExpansionWranglerInterface.__abstractmethods__ = frozenset() - wrglr = ewi.ExpansionWranglerInterface() - try: - wrglr.reorder_sources(None) - wrglr.reorder_potentials(None) - wrglr.finalize_potentials(None) - except NotImplementedError: - assert False - assert True - - -def test_interface_multipole_formation(): - # make the interface instantiable by overriding its set of abstract methods - ewi.ExpansionWranglerInterface.__abstractmethods__ = frozenset() - wrglr = ewi.ExpansionWranglerInterface() - try: - wrglr.form_multipoles( - level_start_source_box_nrs=None, source_boxes=None, src_weights=None - ) - wrglr.form_locals( - level_start_target_or_target_parent_box_nrs=None, - target_or_target_parent_boxes=None, - starts=None, - lists=None, - src_weights=None, - ) - except NotImplementedError: - assert False - assert True - - -def test_interface_multipole_evaluation(): - # make the interface instantiable by overriding its set of abstract methods - ewi.ExpansionWranglerInterface.__abstractmethods__ = frozenset() - wrglr = ewi.ExpansionWranglerInterface() - try: - wrglr.eval_direct( - target_boxes=None, - neighbor_sources_starts=None, - neighbor_sources_lists=None - ) - wrglr.eval_locals( - level_start_target_box_nrs=None, target_boxes=None, local_exps=None - ) - wrglr.eval_multipoles( - level_start_target_box_nrs=None, - target_boxes=None, - starts=None, - lists=None, - mpole_exps=None, - ) - except NotImplementedError: - assert False - assert True - - -def test_interface_multipole_manipulation(): - # make the interface instantiable by overriding its set of abstract methods - ewi.ExpansionWranglerInterface.__abstractmethods__ = frozenset() - wrglr = ewi.ExpansionWranglerInterface() - try: - wrglr.coarsen_multipoles( - level_start_source_parent_box_nrs=None, - source_parent_boxes=None, - mpoles=None, - ) - wrglr.multipole_to_local( - level_start_target_or_target_parent_box_nrs=None, - target_or_target_parent_boxes=None, - starts=None, - lists=None, - mpole_exps=None, - ) - wrglr.refine_locals( - level_start_target_or_target_parent_box_nrs=None, - target_or_target_parent_boxes=None, - local_exps=None, - ) - except NotImplementedError: - assert False - assert True diff --git a/test/test_nearfield_interaction_completeness.py b/test/test_nearfield_interaction_completeness.py index 6dd461f..305838c 100644 --- a/test/test_nearfield_interaction_completeness.py +++ b/test/test_nearfield_interaction_completeness.py @@ -48,10 +48,11 @@ def drive_test_completeness(ctx, queue, dim, q_order): import volumential.meshgen as mg q_points, q_weights, _ = mg.make_uniform_cubic_grid( - nqpoints=q_order+1, level=n_levels, dim=dim) + degree=q_order-1, level=n_levels, dim=dim) - assert len(q_points) == dim - len(q_weights) == len(q_points[0]) + assert len(q_points[0]) == dim + len(q_weights) == len(q_points) + q_points = q_points.T from pytools.obj_array import make_obj_array diff --git a/volumential/droste.py b/volumential/droste.py index e594bf9..4dc4ffc 100644 --- a/volumential/droste.py +++ b/volumential/droste.py @@ -274,7 +274,7 @@ class DrosteBase(KernelCacheWrapper): id=None, assignee="knl_scaling", expression=sympy_conv(self.integral_knl.get_global_scaling_const()), - temp_var_type=lp.Optional(), + temp_var_type=lp.Optional() ) return quad_kernel_insns + [scaling_assignment] @@ -660,7 +660,7 @@ class DrosteBase(KernelCacheWrapper): import volumential.meshgen as mg q_points, _, _ = mg.make_uniform_cubic_grid( - degree=self.ntgt_points, level=1, dim=self.dim, + degree=self.ntgt_points - 1, level=1, dim=self.dim, queue=queue) # map to [0,1]^d diff --git a/volumential/expansion_wrangler_fpnd.py b/volumential/expansion_wrangler_fpnd.py index 4c0cd48..bc80d6b 100644 --- a/volumential/expansion_wrangler_fpnd.py +++ b/volumential/expansion_wrangler_fpnd.py @@ -30,11 +30,10 @@ from pytools.obj_array import make_obj_array # from pytools import memoize_method from volumential.nearfield_potential_table import NearFieldInteractionTable -from volumential.expansion_wrangler_interface import ( - ExpansionWranglerInterface, TreeIndependentDataForWranglerInterface) from sumpy.fmm import (SumpyExpansionWrangler, SumpyTimingFuture, SumpyTreeIndependentDataForWrangler) -from boxtree.pyfmmlib_integration import FMMLibExpansionWrangler +from boxtree.pyfmmlib_integration import (FMMLibExpansionWrangler, + FMMLibTreeIndependentDataForWrangler) from sumpy.kernel import ( LaplaceKernel, HelmholtzKernel, AxisTargetDerivative, @@ -68,9 +67,7 @@ def inverse_id_map(queue, mapped_ids): # {{{ sumpy backend -class FPNDSumpyTreeIndependentDataForWrangler( - TreeIndependentDataForWranglerInterface, - SumpyTreeIndependentDataForWrangler): +class FPNDSumpyTreeIndependentDataForWrangler(SumpyTreeIndependentDataForWrangler): """Objects of this type serve as a place to keep the code needed for ExpansionWrangler if it is using sumpy to perform multipole expansion and manipulations. @@ -86,8 +83,7 @@ class FPNDSumpyTreeIndependentDataForWrangler( return SumpyExpansionWrangler -class FPNDSumpyExpansionWrangler( - ExpansionWranglerInterface, SumpyExpansionWrangler): +class FPNDSumpyExpansionWrangler(SumpyExpansionWrangler): """This expansion wrangler uses "fpnd" strategy. That is, Far field is computed via Particle approximation and Near field is computed Directly. The FMM is performed using sumpy backend. @@ -112,17 +108,13 @@ class FPNDSumpyExpansionWrangler( def __init__( self, - code_container, - queue, - tree, - near_field_table, - dtype, - fmm_level_to_order, - quad_order, - potential_kind=1, + tree_indep, traversal, dtype, fmm_level_to_order, + near_field_table, quad_order, potential_kind=1, source_extra_kwargs=None, kernel_extra_kwargs=None, self_extra_kwargs=None, + translation_classes_data=None, + preprocessed_mpole_dtype=None, list1_extra_kwargs=None, ): """ @@ -131,24 +123,29 @@ class FPNDSumpyExpansionWrangler( 2. a list of tables, when len(target_kernels) = 1 (multiple levels) 3. otherwise, a dictionary from kernel.__repr__() to a list of its tables """ + super().__init__( + tree_indep, traversal, dtype, fmm_level_to_order, + source_extra_kwargs, + kernel_extra_kwargs, + self_extra_kwargs, + translation_classes_data, + preprocessed_mpole_dtype) - self.code = code_container - self.queue = queue - self.tree = tree + self.list1_extra_kwargs = list1_extra_kwargs self.near_field_table = {} # list of tables for a single out kernel if isinstance(near_field_table, list): - assert len(self.code.target_kernels) == 1 + assert len(self.tree_indep.target_kernels) == 1 self.near_field_table[ - self.code.target_kernels[0].__repr__() + self.tree_indep.target_kernels[0].__repr__() ] = near_field_table self.n_tables = len(near_field_table) # single table elif isinstance(near_field_table, NearFieldInteractionTable): - assert len(self.code.target_kernels) == 1 - self.near_field_table[self.code.target_kernels[0].__repr__()] = [ + assert len(self.tree_indep.target_kernels) == 1 + self.near_field_table[self.tree_indep.target_kernels[0].__repr__()] = [ near_field_table ] self.n_tables = 1 @@ -156,7 +153,7 @@ class FPNDSumpyExpansionWrangler( # dictionary of lists of tables elif isinstance(near_field_table, dict): self.n_tables = dict() - for out_knl in self.code.target_kernels: + for out_knl in self.tree_indep.target_kernels: if repr(out_knl) not in near_field_table: raise RuntimeError( "Missing nearfield table for %s." % repr(out_knl)) @@ -177,15 +174,15 @@ class FPNDSumpyExpansionWrangler( self.potential_kind = potential_kind # TODO: make all parameters table-specific (allow using inhomogeneous tables) - kname = repr(self.code.target_kernels[0]) + kname = repr(self.tree_indep.target_kernels[0]) self.root_table_source_box_extent = ( self.near_field_table[kname][0].source_box_extent) table_starting_level = np.round( - np.log(self.tree.root_extent / self.root_table_source_box_extent) + np.log(traversal.tree.root_extent / self.root_table_source_box_extent) / np.log(2) ) - for kid in range(len(self.code.target_kernels)): - kname = self.code.target_kernels[kid].__repr__() + for kid in range(len(self.tree_indep.target_kernels)): + kname = self.tree_indep.target_kernels[kid].__repr__() for lev, table in zip( range(len(self.near_field_table[kname])), self.near_field_table[kname] @@ -210,9 +207,9 @@ class FPNDSumpyExpansionWrangler( if not isinstance(self.n_tables, dict) and self.n_tables > 1: if ( not abs( - int(self.tree.root_extent / table_root_extent) + int(traversal.tree.root_extent / table_root_extent) * table_root_extent - - self.tree.root_extent + - traversal.tree.root_extent ) < 1e-15 ): @@ -225,13 +222,13 @@ class FPNDSumpyExpansionWrangler( if not isinstance(self.n_tables, dict) and self.n_tables > 1: # this checks that the boxes at the highest level are covered if ( - not tree.nlevels + not traversal.tree.nlevels <= len(self.near_field_table[kname]) + table_starting_level ): raise RuntimeError( "Insufficient list of tables: the " "finest level mesh cells at level " - + str(tree.nlevels) + + str(traversal.tree.nlevels) + " are not covered." ) @@ -239,87 +236,19 @@ class FPNDSumpyExpansionWrangler( # deferred until trav.target_boxes is passed when invoking # eval_direct - self.dtype = dtype - - if source_extra_kwargs is None: - source_extra_kwargs = {} - - if kernel_extra_kwargs is None: - kernel_extra_kwargs = {} - - if self_extra_kwargs is None: - self_extra_kwargs = {} - - if list1_extra_kwargs is None: - list1_extra_kwargs = {} - - if not callable(fmm_level_to_order): - raise TypeError("fmm_level_to_order not passed") - - base_kernel = code_container.get_base_kernel() - kernel_arg_set = frozenset(kernel_extra_kwargs.items()) - self.level_orders = [ - fmm_level_to_order(base_kernel, kernel_arg_set, tree, lev) - for lev in range(tree.nlevels) - ] - - # print("Multipole order = ",self.level_orders) - - self.source_extra_kwargs = source_extra_kwargs - self.kernel_extra_kwargs = kernel_extra_kwargs - self.self_extra_kwargs = self_extra_kwargs - self.list1_extra_kwargs = list1_extra_kwargs - - self.extra_kwargs = source_extra_kwargs.copy() - self.extra_kwargs.update(self.kernel_extra_kwargs) - # }}} End constructor # {{{ data vector utilities - def multipole_expansion_zeros(self): - return SumpyExpansionWrangler.multipole_expansion_zeros(self) - - def local_expansion_zeros(self): - return SumpyExpansionWrangler.local_expansion_zeros(self) - - def output_zeros(self): - return SumpyExpansionWrangler.output_zeros(self) - - def reorder_sources(self, source_array): - return SumpyExpansionWrangler.reorder_sources(self, source_array) - def reorder_targets(self, target_array): if not hasattr(self.tree, 'user_target_ids'): - self.tree.user_target_ids = inverse_id_map( - self.queue, self.tree.sorted_target_ids) - return target_array.with_queue(self.queue)[self.tree.user_target_ids] - - def reorder_potentials(self, potentials): - return SumpyExpansionWrangler.reorder_potentials(self, potentials) - - def finalize_potentials(self, potentials): - # return potentials - return SumpyExpansionWrangler.finalize_potentials(self, potentials) + self.traversal.tree.user_target_ids = inverse_id_map( + target_array.queue, self.traversal.tree.sorted_target_ids) + return target_array.with_queue( + target_array.queue)[self.traversal.tree.user_target_ids] # }}} End data vector utilities - # {{{ formation & coarsening of multipoles - - def form_multipoles(self, level_start_source_box_nrs, source_boxes, src_weights): - return SumpyExpansionWrangler.form_multipoles( - self, level_start_source_box_nrs, source_boxes, src_weights - ) - - def coarsen_multipoles( - self, level_start_source_parent_box_nrs, source_parent_boxes, mpoles - ): - return SumpyExpansionWrangler.coarsen_multipoles( - self, level_start_source_parent_box_nrs, source_parent_boxes, mpoles - ) - - # }}} End formation & coarsening of multipoles - # {{{ direct evaluation of near field interactions def eval_direct_single_out_kernel( @@ -484,11 +413,11 @@ class FPNDSumpyExpansionWrangler( ): pot = self.output_zeros() events = [] - for i in range(len(self.code.target_kernels)): + for i in range(len(self.tree_indep.target_kernels)): # print("processing near-field of out_kernel", i) pot[i], evt = self.eval_direct_single_out_kernel( pot[i], - self.code.target_kernels[i], + self.tree_indep.target_kernels[i], target_boxes, neighbor_source_boxes_starts, neighbor_source_boxes_lists, @@ -582,8 +511,7 @@ class FPNDSumpyExpansionWrangler( # {{{ fmmlib backend (for laplace, helmholtz) -class FPNDFMMLibTreeIndependentDataForWrangler( - TreeIndependentDataForWranglerInterface): +class FPNDFMMLibTreeIndependentDataForWrangler(FMMLibTreeIndependentDataForWrangler): """Objects of this type serve as a place to keep the code needed for ExpansionWrangler if it is using fmmlib to perform multipole expansion and manipulations. @@ -607,8 +535,7 @@ class FPNDFMMLibTreeIndependentDataForWrangler( return FPNDFMMLibExpansionWrangler -class FPNDFMMLibExpansionWrangler( - ExpansionWranglerInterface, FMMLibExpansionWrangler): +class FPNDFMMLibExpansionWrangler(FMMLibExpansionWrangler): """This expansion wrangler uses "fpnd" strategy. That is, Far field is computed via Particle approximation and Near field is computed Directly. The FMM is performed using FMMLib backend. @@ -913,22 +840,6 @@ class FPNDFMMLibExpansionWrangler( # }}} End data vector utilities - # {{{ formation & coarsening of multipoles - - def form_multipoles(self, level_start_source_box_nrs, source_boxes, src_weights): - return FMMLibExpansionWrangler.form_multipoles( - self, level_start_source_box_nrs, source_boxes, src_weights - ) - - def coarsen_multipoles( - self, level_start_source_parent_box_nrs, source_parent_boxes, mpoles - ): - return FMMLibExpansionWrangler.coarsen_multipoles( - self, level_start_source_parent_box_nrs, source_parent_boxes, mpoles - ) - - # }}} End formation & coarsening of multipoles - # {{{ direct evaluation of near field interactions def eval_direct_single_out_kernel( diff --git a/volumential/expansion_wrangler_interface.py b/volumential/expansion_wrangler_interface.py deleted file mode 100644 index 5b09290..0000000 --- a/volumential/expansion_wrangler_interface.py +++ /dev/null @@ -1,235 +0,0 @@ -from __future__ import absolute_import, division, print_function - -__copyright__ = "Copyright (C) 2017 - 2018 Xiaoyu Wei" - -__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. -""" - -import logging - -logger = logging.getLogger(__name__) - -from abc import ABCMeta, abstractmethod - -__doc__ = """ -.. autoclass::ExpansionWranglerInterface -""" - -# {{{ expansion wrangler interface - -# NOTE: abstractmethod's signatures (arguement lists) are not enforced - - -class ExpansionWranglerInterface(object): - """ - Abstract expansion handling interface. - The interface is adapted from, and stays compatible with boxtree/fmm. - TODO: Update docstrings - """ - - __metaclass__ = ABCMeta - - @abstractmethod - def multipole_expansion_zeros(self): - """ - Construct arrays to store multipole expansions for all boxes - """ - pass - - @abstractmethod - def local_expansion_zeros(self): - """ - Construct arrays to store multipole expansions for all boxes - """ - pass - - @abstractmethod - def output_zeros(self): - """ - Construct arrays to store potential values for all target points - """ - pass - - @abstractmethod - def reorder_sources(self, source_array): - """ - Return a copy of *source_array* in tree source order. - *source_array* is in user source order. - """ - pass - - @abstractmethod - def reorder_targets(self, source_array): - """ - Return a copy of *target_array* in tree source order. - *target_array* is in user target order. - """ - pass - - @abstractmethod - def reorder_potentials(self, potentials): - """ - Return a copy of *potentials* in user target order. - *source_weights* is in tree target order. - """ - pass - - @abstractmethod - def form_multipoles(self, level_start_source_box_nrs, source_boxes, src_weights): - """ - Return an expansions array containing multipole expansions - in *source_boxes* due to sources with *src_weights*. - """ - pass - - @abstractmethod - def coarsen_multipoles( - self, level_start_source_parent_box_nrs, source_parent_boxes, mpoles - ): - """ - For each box in *source_parent_boxes*, gather (and translate) - the box's children's multipole expansions in *mpoles* and add - the resulting expansion into the box's multipole expansion - in *mpoles*. - - :returns: *mpoles* - """ - pass - - @abstractmethod - def eval_direct( - self, target_boxes, neighbor_sources_starts, neighbor_sources_lists - ): - """ - For each box in *target_boxes*, evaluate the influence of the - neighbor sources due to *src_weights* - - This step amounts to looking up the corresponding entries in a - pre-built table. - - :returns: a new potential array, see :meth:`output_zeros`. - """ - pass - - @abstractmethod - def multipole_to_local( - self, - level_start_target_or_target_parent_box_nrs, - target_or_target_parent_boxes, - starts, - lists, - mpole_exps, - ): - """ - For each box in *target_or_target_parent_boxes*, translate and add - the influence of the multipole expansion in *mpole_exps* into a new - array of local expansions. - - :returns: a new (local) expansion array. - """ - pass - - @abstractmethod - def eval_multipoles( - self, level_start_target_box_nrs, target_boxes, starts, lists, mpole_exps - ): - """ - For each box in *target_boxes*, evaluate the multipole expansion in - *mpole_exps* in the nearby boxes given in *starts* and *lists*, and - return a new potential array. - - :returns: a new potential array, see :meth:`output_zeros`. - """ - pass - - @abstractmethod - def form_locals( - self, - level_start_target_or_target_parent_box_nrs, - target_or_target_parent_boxes, - starts, - lists, - src_weights, - ): - """ - For each box in *target_or_target_parent_boxes*, form local - expansions due to the sources in the nearby boxes given in *starts* and - *lists*, and return a new local expansion array. - - :returns: a new local expansion array - """ - pass - - @abstractmethod - def refine_locals( - self, - level_start_target_or_target_parent_box_nrs, - target_or_target_parent_boxes, - local_exps, - ): - """ - For each box in *child_boxes*, - translate the box's parent's local expansion in *local_exps* and add - the resulting expansion into the box's local expansion in *local_exps*. - - :returns: *local_exps* - """ - pass - - @abstractmethod - def eval_locals(self, level_start_target_box_nrs, target_boxes, local_exps): - """For each box in *target_boxes*, evaluate the local expansion in - *local_exps* and return a new potential array. - - :returns: a new potential array, see :meth:`output_zeros`. - """ - pass - - @abstractmethod - def finalize_potentials(self, potentials): - """ - Postprocess the reordered potentials. This is where global scaling - factors could be applied. - """ - pass - - -# }}} End expansion wrangler interface - -# {{{ tree-independent data for wrangler interface - - -class TreeIndependentDataForWranglerInterface(object): - """ - Abstract expansion code container interface. - The interface is adapted from, and stays compatible with boxtree/fmm. - """ - - __metaclass__ = ABCMeta - - @abstractmethod - def get_wrangler(self, *args, **kwargs): - """Makes a wrangler object. - """ - pass - -# }}} - -# vim: filetype=pyopencl:fdm=marker diff --git a/volumential/meshgen.py b/volumential/meshgen.py index 6c616f6..3532232 100644 --- a/volumential/meshgen.py +++ b/volumential/meshgen.py @@ -36,6 +36,7 @@ Mesh generation import logging import numpy as np import pyopencl as cl +import modepy as mp from pytools.obj_array import make_obj_array logger = logging.getLogger(__name__) @@ -45,106 +46,100 @@ provider = None # {{{ meshgen Python provider -class MeshGenBase(object): +class MeshGenBoxtree(object): """Base class for Meshgen via BoxTree. The interface is similar to the Meshgen via Deal.II, except that the arguments a and b can also be of higher dimensions to allow added flexibility on choosing the bounding box. - This base class cannot be used directly since the boxtree is not - built in the constructor until dimension_specific_setup() is - provided. - - In addition to the common capabilities of the deal.II implementation, - this class also implements native CL getters like get_q_points_dev() - that play well with the other libraries like boxtree. + The interface of this class should be kept in sync with the dealii backend. """ + def __init__(self, degree, nlevels=1, a=-1, b=1, dim=None, queue=None): + """ + Specifying bounding box: + 1. Assign scalar values to a,b for [a, b]^dim + 2. Assign dim-vectors for [a0, b0] X [a1, b1] ... - def __init__(self, degree, nlevels, a=-1, b=1, queue=None): + Note: due to legacy reasons, the degree argument here is the number of + quadrature nodes in 1D (not the polynomial degree). + """ assert degree > 0 assert nlevels > 0 + self.dim = dim self.degree = degree - self.quadrature_formula = LegendreGaussQuadrature(degree - 1) + self.quadrature_formula = mp.LegendreGaussQuadrature(degree - 1) self.nlevels = nlevels - - self.bound_a = np.array([a]).flatten() - self.bound_b = np.array([b]).flatten() - assert len(self.bound_a) == len(self.bound_b) - assert np.all(self.bound_a < self.bound_b) - - self.dim = len(self.bound_a) - self.root_vertex = self.bound_a - self.root_extent = np.max(self.bound_b - self.bound_a) - - if queue is None: - ctx = cl.create_some_context() - self.queue = cl.CommandQueue(ctx) + self.n_box_nodes = self.degree ** self.dim + self.queue = queue + + x = self.quadrature_formula.nodes # nodes in [-1, 1] + w = self.quadrature_formula.weights + assert len(x) == self.degree + self.box_nodes = np.array(np.meshgrid(*((x,) * self.dim), indexing="ij") + ).reshape([self.dim, -1]) + + if self.dim == 1: + self.box_weights = w + elif self.dim == 2: + self.box_weights = w[:, None] @ w[None, :] + elif self.dim == 3: + self.box_weights = (w[:, None] @ w[None, :] + ).reshape(-1)[:, None] @ w[None, :] else: - self.queue = queue - - # plug in dimension-specific details - self.dimension_specific_setup() - self.n_q_points = self.degree ** self.dim + raise ValueError + self.box_weights = self.box_weights.reshape(-1) - self.boxtree = BoxTree() - self.boxtree.generate_uniform_boxtree( - self.queue, - nlevels=self.nlevels, - root_extent=self.root_extent, - root_vertex=self.root_vertex, - ) - self.quadrature = QuadratureOnBoxTree( - self.boxtree, self.quadrature_formula - ) + # make bounding box + if np.array(a).size == 1: + a = np.repeat(a, dim) + if np.array(b).size == 1: + b = np.repeat(b, dim) - def dimension_specific_setup(self): - pass + assert len(a) == dim + assert len(b) == dim + assert all(ai < bi for ai, bi in zip(a, b)) + bbox = [a, b] - def get_q_points_dev(self): - return self.quadrature.get_q_points(self.queue) + self.tob = make_tob_root(self.dim, bbox) + for i in range(self.nlevels - 1): + self.tob = uniformly_refined(self.tob) def get_q_points(self): - q_points_dev = self.get_q_points_dev() - n_all_q_points = len(q_points_dev[0]) - q_points = np.zeros((n_all_q_points, self.dim)) - for d in range(self.dim): - q_points[:, d] = q_points_dev[d].get(self.queue) - return q_points - - def get_q_weights_dev(self): - return self.quadrature.get_q_weights(self.queue) + lfboxes = self.tob.leaf_boxes() + nodes = np.tile(self.box_nodes, (1, len(lfboxes))) + box_shifts = np.repeat( + self.tob.box_centers[:, lfboxes], self.n_box_nodes, axis=1) + box_scales = np.repeat( + self.tob.get_box_size(lfboxes) / 2, self.n_box_nodes) + nodes = nodes * box_scales[None, :] + box_shifts + return nodes.T def get_q_weights(self): - q_weights_dev = self.get_q_weights_dev() - return q_weights_dev.get(self.queue) - - def get_cell_measures_dev(self): - return self.quadrature.get_cell_measures(self.queue) + lfboxes = self.tob.leaf_boxes() + weights = np.tile(self.box_weights, len(lfboxes)) + box_scales = np.repeat( + self.tob.get_box_size(lfboxes) / 2, self.n_box_nodes) + weights = weights * box_scales**self.tob.dim + return weights def get_cell_measures(self): - cell_measures_dev = self.get_cell_measures_dev() - return cell_measures_dev.get(self.queue) - - def get_cell_centers_dev(self): - return self.quadrature.get_cell_centers(self.queue) + lfboxes = self.tob.leaf_boxes() + box_scales = np.repeat( + self.tob.get_box_size(lfboxes) / 2, self.n_box_nodes) + return box_scales**self.tob.dim def get_cell_centers(self): - cell_centers_dev = self.get_cell_centers_dev() - n_active_cells = self.n_active_cells() - cell_centers = np.zeros((n_active_cells, self.dim)) - for d in range(self.dim): - cell_centers[:, d] = cell_centers_dev[d].get(self.queue) - return cell_centers + return self.bob.box_centers def n_cells(self): """Note that this value can be larger than the actual number of cells used in the boxtree. It mainly serves as the bound for iterators on cells. """ - return self.boxtree.nboxes + return self.tob.nboxes def n_active_cells(self): - return self.boxtree.n_active_boxes + return len(self.tob.leaf_boxes()) def update_mesh( self, criteria, top_fraction_of_cells, bottom_fraction_of_cells @@ -156,7 +151,7 @@ class MeshGenBase(object): logging_func("Number of cells: " + str(self.n_cells())) logging_func("Number of active cells: " + str(self.n_active_cells())) logging_func("Number of quad points per cell: " - + str(self.n_q_points)) + + str(self.n_box_nodes)) def generate_gmsh(self, filename): """ @@ -181,10 +176,8 @@ except ImportError as e: logger.warning("Meshgen via Deal.II is not present or unusable.") try: - logger.info("Trying out BoxTree.TreeInteractiveBuild interface.") - from boxtree.tree_build import ( - make_tob_root, uniformly_refined, distribute_quadrature_rule) - import modepy as mp + logger.info("Trying out BoxTree backend.") + from boxtree.tree_build import make_tob_root, uniformly_refined provider = "meshgen_boxtree" @@ -200,7 +193,7 @@ except ImportError as e: def greet(): return "Hello from Meshgen via BoxTree!" - def make_uniform_cubic_grid(nqpoints, nlevels=1, dim=2, actx=None, **kwargs): + def make_uniform_cubic_grid(degree, nlevels=1, dim=2, actx=None, **kwargs): """Uniform cubic grid in [-1,1]^dim. This function provides backward compatibility with meshgen_dealii. """ @@ -208,48 +201,29 @@ except ImportError as e: if "level" in kwargs: nlevels = kwargs["level"] - lower_bounds = [-1, ] * dim - upper_bounds = [1, ] * dim - tob = make_tob_root(dim=dim, bbox=[lower_bounds, upper_bounds]) - for _ in range(nlevels - 1): - tob = uniformly_refined(tob) - - degree = nqpoints - 1 - quad = mp.LegendreGaussQuadrature(degree) - x, q = distribute_quadrature_rule(tob, quad) + mgen = MeshGenBoxtree(degree + 1, nlevels, -1, 1, dim) + x = mgen.get_q_points() + q = mgen.get_q_weights() radii = None return (x, q, radii) - class MeshGen1D(MeshGenBase): + class MeshGen1D(MeshGenBoxtree): """Meshgen in 1D """ + def __init__(self, degree, nlevels=1, a=-1, b=1, queue=None): + super().__init__(degree, nlevels, a, b, 1, queue) - def dimension_specific_setup(self): - assert self.dim == 1 - - class MeshGen2D(MeshGenBase): + class MeshGen2D(MeshGenBoxtree): """Meshgen in 2D """ + def __init__(self, degree, nlevels=1, a=-1, b=1, queue=None): + super().__init__(degree, nlevels, a, b, 2, queue) - def dimension_specific_setup(self): - if self.dim == 1: - # allow passing scalar values of a and b to the constructor - self.dim = 2 - self.root_vertex = np.zeros(self.dim) + self.bound_a - else: - assert self.dim == 2 - - class MeshGen3D(MeshGenBase): + class MeshGen3D(MeshGenBoxtree): """Meshgen in 3D """ - - def dimension_specific_setup(self): - if self.dim == 1: - # allow passing scalar values of a and b to the constructor - self.dim = 3 - self.root_vertex = np.zeros(self.dim) + self.bound_a - else: - assert self.dim == 3 + def __init__(self, degree, nlevels=1, a=-1, b=1, queue=None): + super().__init__(degree, nlevels, a, b, 3, queue) # }}} End Meshgen via BoxTree diff --git a/volumential/nearfield_potential_table.py b/volumential/nearfield_potential_table.py index 5e5c9bb..6b66cb6 100644 --- a/volumential/nearfield_potential_table.py +++ b/volumential/nearfield_potential_table.py @@ -304,14 +304,14 @@ class NearFieldInteractionTable(object): queue = None q_points, _, _ = mg.make_uniform_cubic_grid( - nqpoints=quad_order-1, level=1, dim=self.dim, + degree=quad_order-1, level=1, dim=self.dim, queue=queue) # map to source box mapped_q_points = np.array( [ 0.5 * self.source_box_extent * (qp + np.ones(self.dim)) - for qp in q_points.T + for qp in q_points ] ) # sort in dictionary order, preserve only the leading @@ -401,11 +401,10 @@ class NearFieldInteractionTable(object): idx = self.unwrap_mode_index(mode_index) xi = ( - np.array([p[self.dim - 1] for p in self.q_points[: self.quad_order + 1]]) + np.array([p[self.dim - 1] for p in self.q_points[: self.quad_order]]) / self.source_box_extent ) - assert len(xi) == self.quad_order - 1 - print(xi) + assert len(xi) == self.quad_order yi = [] for d in range(self.dim): diff --git a/volumential/volume_fmm.py b/volumential/volume_fmm.py index 8848b22..15dcef9 100644 --- a/volumential/volume_fmm.py +++ b/volumential/volume_fmm.py @@ -31,7 +31,6 @@ import numpy as np import pyopencl as cl from pytools.obj_array import make_obj_array from boxtree.timing import TimingRecorder -from volumential.expansion_wrangler_interface import ExpansionWranglerInterface from volumential.expansion_wrangler_fpnd import ( FPNDSumpyExpansionWrangler, FPNDFMMLibExpansionWrangler) @@ -75,7 +74,6 @@ def drive_volume_fmm(traversal, expansion_wrangler, src_weights, src_func, Returns the potentials computed by *expansion_wrangler*. """ wrangler = expansion_wrangler - assert issubclass(type(wrangler), ExpansionWranglerInterface) recorder = TimingRecorder() logger.info("start fmm") -- GitLab From 2ef82df0222681a39e2a0e78bc24d7b6891ee057 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 9 Mar 2022 11:32:08 -0600 Subject: [PATCH 3/5] Fix loopy kernel codegen --- volumential/droste.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/volumential/droste.py b/volumential/droste.py index 4dc4ffc..85ad38f 100644 --- a/volumential/droste.py +++ b/volumential/droste.py @@ -274,7 +274,7 @@ class DrosteBase(KernelCacheWrapper): id=None, assignee="knl_scaling", expression=sympy_conv(self.integral_knl.get_global_scaling_const()), - temp_var_type=lp.Optional() + temp_var_type=lp.Optional(None) ) return quad_kernel_insns + [scaling_assignment] @@ -787,6 +787,10 @@ class DrosteFull(DrosteBase): loopy_knl = lp.set_options(loopy_knl, write_cl=False) loopy_knl = lp.set_options(loopy_knl, return_dict=True) + loopy_knl = lp.add_and_infer_dtypes( + loopy_knl, + {"template_target": np.float64}) + try: loopy_knl = self.integral_knl.prepare_loopy_kernel(loopy_knl) except Exception: # noqa: B902 @@ -1239,6 +1243,10 @@ class DrosteReduced(DrosteBase): lang_version=(2018, 2), ) + loopy_knl = lp.add_and_infer_dtypes( + loopy_knl, + {"template_target": np.float64}) + elif self.get_kernel_id == 1: loopy_knl = lp.make_kernel( # NOQA [domain], -- GitLab From 5d8315afe640c50e5d6f98618431466330d9b7d8 Mon Sep 17 00:00:00 2001 From: xywei Date: Wed, 9 Mar 2022 16:11:01 -0600 Subject: [PATCH 4/5] Misc fixes --- test/test_droste.py | 3 +- ...test_nearfield_interaction_completeness.py | 10 ++--- test/test_nearfield_potential_table.py | 4 +- volumential/droste.py | 13 +++++- volumential/expansion_wrangler_fpnd.py | 40 ++++++++++--------- volumential/list1.py | 10 ++--- volumential/nearfield_potential_table.py | 4 +- 7 files changed, 46 insertions(+), 38 deletions(-) diff --git a/test/test_droste.py b/test/test_droste.py index 1363db5..67fb98c 100644 --- a/test/test_droste.py +++ b/test/test_droste.py @@ -26,6 +26,7 @@ import pytest import loopy as lp import numpy as np import pyopencl as cl +import pyopencl.array from volumential.droste import DrosteReduced from numpy.polynomial.chebyshev import chebval @@ -102,8 +103,6 @@ def drive_test_cheb_table( n_brick_quad_points=n_brick_q_points, adaptive_level=False, adaptive_quadrature=False) - print(cheb_table) - # handle AxisTargetDerivative to extract the component of interest # (pvfmm returns the full gradient) from sumpy.kernel import AxisTargetDerivative diff --git a/test/test_nearfield_interaction_completeness.py b/test/test_nearfield_interaction_completeness.py index 305838c..b0f5211 100644 --- a/test/test_nearfield_interaction_completeness.py +++ b/test/test_nearfield_interaction_completeness.py @@ -129,10 +129,10 @@ def drive_test_completeness(ctx, queue, dim, q_order): mpole_expn_class = partial(VolumeTaylorMultipoleExpansion, use_rscale=None) from volumential.expansion_wrangler_fpnd import ( - FPNDExpansionWranglerCodeContainer, + FPNDSumpyTreeIndependentDataForWrangler, FPNDExpansionWrangler) - wcc = FPNDExpansionWranglerCodeContainer( + wcc = FPNDSumpyTreeIndependentDataForWrangler( ctx, partial(mpole_expn_class, knl), partial(local_expn_class, knl), @@ -141,9 +141,8 @@ def drive_test_completeness(ctx, queue, dim, q_order): ) wrangler = FPNDExpansionWrangler( - code_container=wcc, - queue=queue, - tree=tree, + tree_indep=wcc, + traversal=trav, near_field_table=nft, dtype=dtype, fmm_level_to_order=lambda kernel, kernel_args, tree, lev: 1, @@ -155,6 +154,7 @@ def drive_test_completeness(ctx, queue, dim, q_order): trav.target_boxes, trav.neighbor_source_boxes_starts, trav.neighbor_source_boxes_lists, mode_coefs=source_vals) pot = pot[0] + import pudb; pu.db for p in pot[0]: assert(abs(p - 2**dim) < 1e-8) diff --git a/test/test_nearfield_potential_table.py b/test/test_nearfield_potential_table.py index ec87c38..6b0af2e 100644 --- a/test/test_nearfield_potential_table.py +++ b/test/test_nearfield_potential_table.py @@ -59,8 +59,6 @@ def interp_modes(q_order): yi = yi.flatten() val = np.zeros(xi.shape) - print(xi) - print(yi) for i in range(len(xi)): val[i] = interpolate_function(xi[i], yi[i]) @@ -69,7 +67,7 @@ def interp_modes(q_order): def test_modes(): - for q in [1, 2, 3, 4, 5, 6, 7, 8]: + for q in [1, 2, 3, 4, 5, 6]: val = interp_modes(q) assert np.allclose(val, 1) diff --git a/volumential/droste.py b/volumential/droste.py index 85ad38f..f4b0712 100644 --- a/volumential/droste.py +++ b/volumential/droste.py @@ -660,7 +660,7 @@ class DrosteBase(KernelCacheWrapper): import volumential.meshgen as mg q_points, _, _ = mg.make_uniform_cubic_grid( - degree=self.ntgt_points - 1, level=1, dim=self.dim, + degree=self.ntgt_points, level=1, dim=self.dim, queue=queue) # map to [0,1]^d @@ -1371,7 +1371,15 @@ class DrosteReduced(DrosteBase): delattr(self, "_memoize_dic_get_cached_optimized_kernel") except Exception: # noqa: B902 pass - knl = self.get_cached_optimized_kernel() + + if 1: + knl = self.get_cached_optimized_kernel() + else: + knl = self.get_kernel() + # knl = lp.set_options(knl, {"trace_assignment_values": True}) + print(knl) + 1/0 + evt, res = knl( queue, alpha=alpha, result=result_array, root_brick=root_brick, @@ -1384,6 +1392,7 @@ class DrosteReduced(DrosteBase): raw_cheb_table_case = res["result"] + self.get_kernel_id = 1 try: delattr(self, "_memoize_dic_get_cached_optimized_kernel") diff --git a/volumential/expansion_wrangler_fpnd.py b/volumential/expansion_wrangler_fpnd.py index bc80d6b..5d33617 100644 --- a/volumential/expansion_wrangler_fpnd.py +++ b/volumential/expansion_wrangler_fpnd.py @@ -131,7 +131,8 @@ class FPNDSumpyExpansionWrangler(SumpyExpansionWrangler): translation_classes_data, preprocessed_mpole_dtype) - self.list1_extra_kwargs = list1_extra_kwargs + self.list1_extra_kwargs = \ + dict() if list1_extra_kwargs is None else list1_extra_kwargs self.near_field_table = {} # list of tables for a single out kernel @@ -267,7 +268,8 @@ class FPNDSumpyExpansionWrangler(SumpyExpansionWrangler): if 0: print("Returns range for list1") - out_pot[:] = cl.array.to_device(self.queue, np.arange(len(out_pot))) + out_pot[:] = cl.array.to_device( + mode_coefs[0].queue, np.arange(len(out_pot))) return out_pot, None kname = out_kernel.__repr__() @@ -279,11 +281,13 @@ class FPNDSumpyExpansionWrangler(SumpyExpansionWrangler): else: use_multilevel_tables = False + queue = mode_coefs[0].queue + if use_multilevel_tables: # this checks that the boxes at the coarsest level # and allows for some round-off error min_lev = np.min( - self.tree.box_levels.get(self.queue)[target_boxes.get(self.queue)] + self.tree.box_levels.get(queue)[target_boxes.get(queue)] ) largest_cell_extent = self.tree.root_extent * 0.5 ** min_lev if not self.near_field_table[kname][0].source_box_extent >= ( @@ -305,7 +309,7 @@ class FPNDSumpyExpansionWrangler(SumpyExpansionWrangler): shift = -min(distinct_numbers) case_indices_dev = cl.array.to_device( - self.queue, self.near_field_table[kname][0].case_indices + queue, self.near_field_table[kname][0].case_indices ) # table.data @@ -335,7 +339,7 @@ class FPNDSumpyExpansionWrangler(SumpyExpansionWrangler): exterior_mode_nmlz_combined[lev, :] = \ self.near_field_table[kname][lev].kernel_exterior_normalizers - self.queue.finish() + queue.finish() logger.info( "table data for kernel " + out_kernel.__repr__() + " congregated") @@ -356,20 +360,17 @@ class FPNDSumpyExpansionWrangler(SumpyExpansionWrangler): potential_kind=self.potential_kind, **self.list1_extra_kwargs) - table_data_combined = cl.array.to_device(self.queue, - table_data_combined) - mode_nmlz_combined = cl.array.to_device(self.queue, - mode_nmlz_combined) - exterior_mode_nmlz_combined = cl.array.to_device(self.queue, - exterior_mode_nmlz_combined) - self.queue.finish() + table_data_combined = cl.array.to_device(queue, table_data_combined) + mode_nmlz_combined = cl.array.to_device(queue, mode_nmlz_combined) + exterior_mode_nmlz_combined = cl.array.to_device(queue, exterior_mode_nmlz_combined) + queue.finish() logger.info("sent table data to device") # NOTE: box_sources for this evaluation should be "box_targets". # This is due to the special features of how box-FMM works. res, evt = near_field( - self.queue, + queue, result=out_pot, box_centers=self.tree.box_centers, box_levels=self.tree.box_levels, @@ -411,7 +412,7 @@ class FPNDSumpyExpansionWrangler(SumpyExpansionWrangler): neighbor_source_boxes_lists, mode_coefs, ): - pot = self.output_zeros() + pot = self.output_zeros(mode_coefs[0]) events = [] for i in range(len(self.tree_indep.target_kernels)): # print("processing near-field of out_kernel", i) @@ -428,7 +429,7 @@ class FPNDSumpyExpansionWrangler(SumpyExpansionWrangler): for out_pot in pot: out_pot.finish() - return (pot, SumpyTimingFuture(self.queue, events)) + return (pot, SumpyTimingFuture(mode_coefs[0].queue, events)) # }}} End direct evaluation of near field interactions @@ -870,11 +871,13 @@ class FPNDFMMLibExpansionWrangler(FMMLibExpansionWrangler): else: use_multilevel_tables = False + queue = mode_coefs[0].queue + if use_multilevel_tables: # this checks that the boxes at the coarsest level # and allows for some round-off error min_lev = np.min( - self.tree.box_levels.get(self.queue)[target_boxes.get(self.queue)] + self.tree.box_levels.get(queue)[target_boxes.get(queue)] ) largest_cell_extent = self.tree.root_extent * 0.5 ** min_lev if not self.near_field_table[kname][0].source_box_extent >= ( @@ -896,8 +899,7 @@ class FPNDFMMLibExpansionWrangler(FMMLibExpansionWrangler): shift = -min(distinct_numbers) case_indices_dev = cl.array.to_device( - self.queue, self.near_field_table[kname][0].case_indices - ) + queue, self.near_field_table[kname][0].case_indices) # table.data table_data_combined = np.zeros( @@ -952,7 +954,7 @@ class FPNDFMMLibExpansionWrangler(FMMLibExpansionWrangler): **self.list1_extra_kwargs) res, evt = near_field( - self.queue, + queue, result=out_pot, box_centers=self.tree.box_centers, box_levels=self.tree.box_levels, diff --git a/volumential/list1.py b/volumential/list1.py index e845589..9355868 100644 --- a/volumential/list1.py +++ b/volumential/list1.py @@ -288,7 +288,7 @@ class NearFieldFromCSR(NearFieldEvalBase): else: logger.info("computing table level from box size") logger.info("(using multiple tables)") - code = "log(table_root_extent / BOX_extent) / log(2.0)" + code = "log(table_root_extent / BOX_extent) / log(2.0) + 0.5" return code.replace("BOX", box_name) @@ -342,10 +342,10 @@ class NearFieldFromCSR(NearFieldEvalBase): <> sbox_extent = root_extent * (1.0 / (2**sbox_level)) table_lev_tmp = GET_TABLE_LEVEL {id=tab_lev_tmp} - table_lev = round(table_lev_tmp) {id=tab_lev,dep=tab_lev_tmp} + table_lev = floor(table_lev_tmp) {id=tab_lev,dep=tab_lev_tmp} vec_id_tmp = COMPUTE_VEC_ID {id=vec_id_tmp} - vec_id = round(vec_id_tmp) {id=vec_id,dep=vec_id_tmp} + vec_id = floor(vec_id_tmp) {id=vec_id,dep=vec_id_tmp} <> case_id = case_indices[vec_id] {dep=vec_id} <> scaling = COMPUTE_SCALING @@ -355,7 +355,7 @@ class NearFieldFromCSR(NearFieldEvalBase): <> tgt_scaling = COMPUTE_TGT_SCALING <> tgt_displacement = COMPUTE_TGT_DISPLACEMENT tgt_table_lev_tmp = GET_TGT_TABLE_LEVEL {id=tgttab_lev_tmp} - tgt_table_lev = round(tgt_table_lev_tmp) \ + tgt_table_lev = floor(tgt_table_lev_tmp) \ {id=tgttab_lev,dep=tgttab_lev_tmp} <> ext_nmlz = exterior_mode_nmlz[tgt_table_lev, tid] \ * tgt_scaling + tgt_displacement \ @@ -502,7 +502,7 @@ class NearFieldFromCSR(NearFieldEvalBase): extra_knl_args_from_init = {} for key, val in integral_kernel_init_kargs.items(): - if key in knl.arg_dict: + if key in knl.default_entrypoint.arg_dict: extra_knl_args_from_init[key] = val evt, res = knl( diff --git a/volumential/nearfield_potential_table.py b/volumential/nearfield_potential_table.py index 6b66cb6..f49deee 100644 --- a/volumential/nearfield_potential_table.py +++ b/volumential/nearfield_potential_table.py @@ -304,7 +304,7 @@ class NearFieldInteractionTable(object): queue = None q_points, _, _ = mg.make_uniform_cubic_grid( - degree=quad_order-1, level=1, dim=self.dim, + degree=quad_order, level=1, dim=self.dim, queue=queue) # map to source box @@ -427,7 +427,7 @@ class NearFieldInteractionTable(object): def get_mode(self, mode_index): """ - normal modes are deined on the source box + normal modes are defined on the source box """ assert mode_index >= 0 and mode_index < self.n_q_points -- GitLab From 2d9ae026229af1f9c9845c8e2ce0cbd209eca93c Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 13 Mar 2022 20:59:15 -0500 Subject: [PATCH 5/5] Revive laplace --- examples/laplace2d.py | 1 + examples/laplace3d.py | 11 +++++------ volumential/volume_fmm.py | 23 +++++++++++------------ 3 files changed, 17 insertions(+), 18 deletions(-) diff --git a/examples/laplace2d.py b/examples/laplace2d.py index 07a9271..9de7d76 100644 --- a/examples/laplace2d.py +++ b/examples/laplace2d.py @@ -38,6 +38,7 @@ else: import numpy as np import pyopencl as cl +import pyopencl.array from volumential.tools import ScalarFieldExpressionEvaluation as Eval import pymbolic as pmbl diff --git a/examples/laplace3d.py b/examples/laplace3d.py index 31e9615..7b4f37f 100644 --- a/examples/laplace3d.py +++ b/examples/laplace3d.py @@ -29,6 +29,7 @@ import logging import numpy as np import pyopencl as cl +import pyopencl.array import pymbolic as pmbl import pymbolic.functions @@ -293,10 +294,9 @@ def main(): exclude_self = True from volumential.expansion_wrangler_fpnd import ( - FPNDExpansionWrangler, - FPNDExpansionWranglerCodeContainer) + FPNDSumpyTreeIndependentDataForWrangler, FPNDExpansionWrangler) - wcc = FPNDExpansionWranglerCodeContainer( + wcc = FPNDSumpyTreeIndependentDataForWrangler( ctx, partial(mpole_expn_class, knl), partial(local_expn_class, knl), @@ -311,9 +311,8 @@ def main(): self_extra_kwargs = {} wrangler = FPNDExpansionWrangler( - code_container=wcc, - queue=queue, - tree=tree, + tree_indep=wcc, + traversal=trav, near_field_table=nftable, dtype=dtype, fmm_level_to_order=lambda kernel, kernel_args, tree, lev: m_order, diff --git a/volumential/volume_fmm.py b/volumential/volume_fmm.py index 15dcef9..aba4c8b 100644 --- a/volumential/volume_fmm.py +++ b/volumential/volume_fmm.py @@ -85,6 +85,8 @@ def drive_volume_fmm(traversal, expansion_wrangler, src_weights, src_func, if src_func.ndim == 1: src_func = make_obj_array([src_func.astype(dtype)]) + queue = src_func[0].queue + assert (ns := len(src_weights)) == len(src_func) if ns > 1: raise NotImplementedError("Multiple outputs are not yet supported") @@ -95,12 +97,12 @@ def drive_volume_fmm(traversal, expansion_wrangler, src_weights, src_func, elif isinstance(expansion_wrangler, FPNDFMMLibExpansionWrangler): - traversal = traversal.get(wrangler.queue) + traversal = traversal.get(queue) if isinstance(src_weights, cl.array.Array): - src_weights = src_weights.get(wrangler.queue) + src_weights = src_weights.get(queue) if isinstance(src_func, cl.array.Array): - src_func = src_func.get(wrangler.queue) + src_func = src_func.get(queue) if reorder_sources: logger.debug("reorder source weights") @@ -187,7 +189,7 @@ def drive_volume_fmm(traversal, expansion_wrangler, src_weights, src_func, from sumpy import P2P p2p = P2P( - wrangler.queue.context, + queue.context, wrangler.code.target_kernels, exclude_self=wrangler.code.exclude_self, value_dtypes=[wrangler.dtype], @@ -203,7 +205,7 @@ def drive_volume_fmm(traversal, expansion_wrangler, src_weights, src_func, for iw, sw in enumerate(src_weights): evt, (ref_pot,) = p2p( - wrangler.queue, + queue, traversal.tree.targets, traversal.tree.sources, (sw,), @@ -231,7 +233,7 @@ def drive_volume_fmm(traversal, expansion_wrangler, src_weights, src_func, result = wrangler.reorder_potentials(potentials) logger.debug("finalize potentials") - result = wrangler.finalize_potentials(result) + result = wrangler.finalize_potentials(result, None) logger.info("direct p2p complete") @@ -350,7 +352,7 @@ def drive_volume_fmm(traversal, expansion_wrangler, src_weights, src_func, result = wrangler.reorder_potentials(potentials) logger.debug("finalize potentials") - result = wrangler.finalize_potentials(result) + result = wrangler.finalize_potentials(result, None) logger.info("fmm complete") # }}} End Reorder potentials @@ -388,9 +390,7 @@ def interpolate_volume_potential(target_points, traversal, wrangler, potential, """ Interpolate the volume potential, only works for tensor-product quadrature formulae. target_points and potential should be an cl array. - - :arg wrangler: Used only for general info (nothing sumpy kernel specific). May also - be None if the needed information is passed by kwargs. +https://github.com/riedat/gc be None if the needed information is passed by kwargs. :arg potential_in_tree_order: Whether the potential is in tree order (as opposed to in user order). :lbl_lookup: a :class:`boxtree.LeavesToBallsLookup` that has the lookup information @@ -401,16 +401,15 @@ def interpolate_volume_potential(target_points, traversal, wrangler, potential, if wrangler is not None: dim = next(iter(wrangler.near_field_table.values()))[0].dim tree = wrangler.tree - queue = wrangler.queue q_order = wrangler.quad_order dtype = wrangler.dtype else: dim = kwargs['dim'] tree = kwargs['tree'] - queue = kwargs['queue'] q_order = kwargs['q_order'] dtype = kwargs['dtype'] + queue = potential.queue ctx = queue.context coord_dtype = tree.coord_dtype n_points = len(target_points[0]) -- GitLab