diff --git a/.gitignore b/.gitignore index 00ac7501afc884f227b151e60977a1f9e352365c..80e0de7a0c5959778f5ea374fde26993bc637c58 100644 --- a/.gitignore +++ b/.gitignore @@ -39,3 +39,5 @@ gmsh_tmp/ .coverage pytest.xml prof/ + +.ccls-cache/ diff --git a/examples/laplace2d.py b/examples/laplace2d.py index 383951f7341681a4a88de9332602bb7c54d435cd..9de7d76275202f926e10aa348fc1d819463accd3 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 @@ -280,11 +281,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 +298,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/examples/laplace3d.py b/examples/laplace3d.py index 31e96153e26258578217512b0c79d21d6a304678..7b4f37fd9bc25ad8f814f38f61f4f5e25cb1fcfa 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/requirements.txt b/requirements.txt index 099dcc81e227b79587a25dc924f8d9f73d9f9d47..3c859ba2c92106db388c6c746566ddcee4c01fce 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_droste.py b/test/test_droste.py index 1363db5938804a9b300f0a3f3cfb5c596d2f5c5c..67fb98c1031ac8625a64113a359dc995759e8a11 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_expansion_wrangler_interface.py b/test/test_expansion_wrangler_interface.py deleted file mode 100644 index 11ade594bf059556fcfc4c163cca5a5453b35c45..0000000000000000000000000000000000000000 --- 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 cd590c7c803710b5481656072a1c102949fa6698..b0f52119ac501381aec23b141a71d8d44a28cfd4 100644 --- a/test/test_nearfield_interaction_completeness.py +++ b/test/test_nearfield_interaction_completeness.py @@ -41,41 +41,30 @@ 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( + degree=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[0]) == dim + len(q_weights) == len(q_points) + q_points = q_points.T 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 @@ -140,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), @@ -152,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, @@ -166,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 ec87c3840a907a9a25711ee56313ed14ba70c0e3..6b0af2ea0ee02903c2e0c4df190bfa0d2ce29885 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 e594bf9fe5ecdda4a412a4e6decb7e03f820b5ad..f4b0712f3ac8fc485795cbf8e5c755187ed19e33 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], @@ -1363,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, @@ -1376,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 b691f63f4d7edd57cb36c2dd68f229d47f2d88a6..5d3361748ac4536f1e0596f30824a912d4704829 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, ExpansionWranglerCodeContainerInterface) -from sumpy.fmm import SumpyExpansionWrangler, \ - SumpyTimingFuture, SumpyExpansionWranglerCodeContainer -from boxtree.pyfmmlib_integration import FMMLibExpansionWrangler +from sumpy.fmm import (SumpyExpansionWrangler, + SumpyTimingFuture, SumpyTreeIndependentDataForWrangler) +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 FPNDSumpyExpansionWranglerCodeContainer( - ExpansionWranglerCodeContainerInterface, - SumpyExpansionWranglerCodeContainer): +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. @@ -80,11 +77,13 @@ 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( - 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. @@ -109,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, ): """ @@ -128,24 +123,30 @@ 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 = \ + dict() if list1_extra_kwargs is None else 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 @@ -153,7 +154,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)) @@ -174,15 +175,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] @@ -207,9 +208,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 ): @@ -222,13 +223,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." ) @@ -236,87 +237,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( @@ -335,7 +268,8 @@ class FPNDSumpyExpansionWrangler( 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__() @@ -347,11 +281,13 @@ class FPNDSumpyExpansionWrangler( 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 >= ( @@ -373,7 +309,7 @@ class FPNDSumpyExpansionWrangler( 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 @@ -403,7 +339,7 @@ class FPNDSumpyExpansionWrangler( 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") @@ -424,20 +360,17 @@ class FPNDSumpyExpansionWrangler( 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, @@ -479,13 +412,13 @@ class FPNDSumpyExpansionWrangler( neighbor_source_boxes_lists, mode_coefs, ): - pot = self.output_zeros() + pot = self.output_zeros(mode_coefs[0]) 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, @@ -496,7 +429,7 @@ class FPNDSumpyExpansionWrangler( 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 @@ -579,9 +512,7 @@ class FPNDSumpyExpansionWrangler( # {{{ fmmlib backend (for laplace, helmholtz) -class FPNDFMMLibExpansionWranglerCodeContainer( - ExpansionWranglerCodeContainerInterface, - ): +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. @@ -600,17 +531,12 @@ 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( - 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. @@ -915,22 +841,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( @@ -961,11 +871,13 @@ class FPNDFMMLibExpansionWrangler( 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 >= ( @@ -987,8 +899,7 @@ class FPNDFMMLibExpansionWrangler( 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( @@ -1043,7 +954,7 @@ class FPNDFMMLibExpansionWrangler( **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, @@ -1201,7 +1112,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 deleted file mode 100644 index 4418f13f8cb1e35344e146669a2e6ed8a40511e9..0000000000000000000000000000000000000000 --- 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 - -# {{{ code container interface - - -class ExpansionWranglerCodeContainerInterface(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 - -# }}} End code container interface - -# vim: filetype=pyopencl:fdm=marker diff --git a/volumential/list1.py b/volumential/list1.py index e845589c9bc7973dfb82dae2555e884de6070ca4..93558682161b8183a6e6c5d508ef6d9b880e1445 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/meshgen.py b/volumential/meshgen.py index cc936b2367b10446d2125dcce43e74c954252c61..3532232122de8d074ccb5092c6e288d12242d8d6 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 + raise ValueError + self.box_weights = self.box_weights.reshape(-1) - # plug in dimension-specific details - self.dimension_specific_setup() - self.n_q_points = self.degree ** self.dim + # 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) - 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 - ) + assert len(a) == dim + assert len(b) == dim + assert all(ai < bi for ai, bi in zip(a, b)) + bbox = [a, b] - def dimension_specific_setup(self): - pass - - 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,9 +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.") - import boxtree.tree_interactive_build # noqa: F401 - from modepy import LegendreGaussQuadrature + logger.info("Trying out BoxTree backend.") + from boxtree.tree_build import make_tob_root, uniformly_refined provider = "meshgen_boxtree" @@ -195,69 +189,41 @@ 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(degree, 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) - - class MeshGen1D(MeshGenBase): + 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(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 8cef4d2cf84de2e55e9d91c22124dbdfaaae9a27..f49deee716545327ce58ed8188e6d005d3c9d321 100644 --- a/volumential/nearfield_potential_table.py +++ b/volumential/nearfield_potential_table.py @@ -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 diff --git a/volumential/volume_fmm.py b/volumential/volume_fmm.py index 0b1a451ee4b973938625199acd6a30f9d2c59cb7..aba4c8b8bdf648c16af9d34dd54f0266cdffc698 100644 --- a/volumential/volume_fmm.py +++ b/volumential/volume_fmm.py @@ -30,8 +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 volumential.expansion_wrangler_interface import ExpansionWranglerInterface +from boxtree.timing import TimingRecorder 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") @@ -87,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") @@ -97,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") @@ -189,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], @@ -205,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,), @@ -233,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") @@ -352,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 @@ -390,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 @@ -403,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])