diff --git a/contrib/meshgen11_dealii/CMakeLists.txt b/contrib/meshgen11_dealii/CMakeLists.txt index b7c9ea92ec7f201ec81e7b0e7109cbb65043027a..3073e506eb612b19796d54aec217727dde3f209d 100644 --- a/contrib/meshgen11_dealii/CMakeLists.txt +++ b/contrib/meshgen11_dealii/CMakeLists.txt @@ -41,9 +41,8 @@ FIND_PACKAGE(deal.II REQUIRED DEAL_II_INITIALIZE_CACHED_VARIABLES() -# Avoid using 3.7 -find_package( PythonInterp 3.6 EXACT ) -find_package( PythonLibs 3.6 EXACT ) +find_package( PythonInterp 3.8 EXACT ) +find_package( PythonLibs 3.8 EXACT ) add_subdirectory(pybind11) pybind11_add_module(meshgen11 meshgen.cpp) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index d38cd85d601f90d12ef967e4a574299336cb752b..9135c2109d0c48c392f2fb7a503912094bf2cae1 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -70,26 +70,24 @@ def eval_func_on_discr_nodes(queue, discr, func): # }}} End test data -def drive_test_from_meshmode_exact_interpolation_2d( - cl_ctx, queue, +def drive_test_from_meshmode_interpolation( + cl_ctx, queue, dim, degree, nel_1d, n_levels, q_order, - a=-0.5, b=0.5, seed=0): + a=-0.5, b=0.5, seed=0, test_case='exact'): """ meshmode mesh control: nel_1d, degree volumential mesh control: n_levels, q_order """ - dim = 2 - mesh = generate_regular_rect_mesh( - a=(a, a), - b=(b, b), - n=(nel_1d, nel_1d)) + a=(a, ) * dim, + b=(b, ) * dim, + n=(nel_1d, ) * dim) arr_ctx = PyOpenCLArrayContext(queue) group_factory = PolynomialWarpAndBlendGroupFactory(order=degree) discr = Discretization(arr_ctx, mesh, group_factory) - bbox_fac = BoundingBoxFactory(dim=2) + bbox_fac = BoundingBoxFactory(dim=dim) boxfmm_fac = BoxFMMGeometryFactory( cl_ctx, dim=dim, order=q_order, nlevels=n_levels, bbox_getter=bbox_fac, @@ -99,8 +97,22 @@ def drive_test_from_meshmode_exact_interpolation_2d( cl_ctx, tree=boxgeo.tree, discr=discr) lookup, evt = lookup_fac(queue) - # algebraically exact interpolation - func = random_polynomial_func(dim, degree, seed) + if test_case == 'exact': + # algebraically exact interpolation + func = random_polynomial_func(dim, degree, seed) + elif test_case == 'non-exact': + if dim == 2: + def func(pts): + x, y = pts + return np.sin(x + y) * x + elif dim == 3: + def func(pts): + x, y, z = pts + return np.sin(x + y * z) * x + else: + raise ValueError() + else: + raise ValueError() dof_vec = eval_func_on_discr_nodes(queue, discr, func) res = interpolate_from_meshmode(queue, dof_vec, lookup).get(queue) @@ -108,72 +120,31 @@ def drive_test_from_meshmode_exact_interpolation_2d( tree = boxgeo.tree.get(queue) ref = func(np.vstack(tree.sources)) - return np.allclose(ref, res) - - -def drive_test_from_meshmode_interpolation_2d( - cl_ctx, queue, - degree, nel_1d, n_levels, q_order, - a=-0.5, b=0.5, seed=0): - """ - meshmode mesh control: nel_1d, degree - volumential mesh control: n_levels, q_order - """ - dim = 2 - - mesh = generate_regular_rect_mesh( - a=(a, a), - b=(b, b), - n=(nel_1d, nel_1d)) - - arr_ctx = PyOpenCLArrayContext(queue) - group_factory = PolynomialWarpAndBlendGroupFactory(order=degree) - discr = Discretization(arr_ctx, mesh, group_factory) - - bbox_fac = BoundingBoxFactory(dim=2) - boxfmm_fac = BoxFMMGeometryFactory( - cl_ctx, dim=dim, order=q_order, - nlevels=n_levels, bbox_getter=bbox_fac, - expand_to_hold_mesh=mesh, mesh_padding_factor=0.) - boxgeo = boxfmm_fac(queue) - lookup_fac = ElementsToSourcesLookupBuilder( - cl_ctx, tree=boxgeo.tree, discr=discr) - lookup, evt = lookup_fac(queue) - - def func(pts): - x, y = pts - return np.sin(x + y) * x - - dof_vec = eval_func_on_discr_nodes(queue, discr, func) - res = interpolate_from_meshmode(queue, dof_vec, lookup).get(queue) - - tree = boxgeo.tree.get(queue) - ref = func(np.vstack(tree.sources)) + if test_case == 'exact': + return np.allclose(ref, res) resid = np.linalg.norm(ref - res, ord=np.inf) return resid -def drive_test_to_meshmode_exact_interpolation_2d( - cl_ctx, queue, +def drive_test_to_meshmode_interpolation( + cl_ctx, queue, dim, degree, nel_1d, n_levels, q_order, - a=-0.5, b=0.5, seed=0): + a=-0.5, b=0.5, seed=0, test_case='exact'): """ meshmode mesh control: nel_1d, degree volumential mesh control: n_levels, q_order """ - dim = 2 - mesh = generate_regular_rect_mesh( - a=(a, a), - b=(b, b), - n=(nel_1d, nel_1d)) + a=(a, ) * dim, + b=(b, ) * dim, + n=(nel_1d, ) * dim) arr_ctx = PyOpenCLArrayContext(queue) group_factory = PolynomialWarpAndBlendGroupFactory(order=degree) discr = Discretization(arr_ctx, mesh, group_factory) - bbox_fac = BoundingBoxFactory(dim=2) + bbox_fac = BoundingBoxFactory(dim=dim) boxfmm_fac = BoxFMMGeometryFactory( cl_ctx, dim=dim, order=q_order, nlevels=n_levels, bbox_getter=bbox_fac, @@ -183,59 +154,37 @@ def drive_test_to_meshmode_exact_interpolation_2d( cl_ctx, trav=boxgeo.trav, discr=discr) lookup, evt = lookup_fac(queue) - # algebraically exact interpolation - func = random_polynomial_func(dim, degree, seed) + if test_case == 'exact': + # algebraically exact interpolation + func = random_polynomial_func(dim, degree, seed) + elif test_case == 'non-exact': + if dim == 2: + def func(pts): + x, y = pts + return np.sin(x + y) * x + elif dim == 3: + def func(pts): + x, y, z = pts + return np.sin(x + y * z) * x + else: + raise ValueError() + else: + raise ValueError() tree = boxgeo.tree.get(queue) potential = func(np.vstack(tree.sources)) res = interpolate_to_meshmode(queue, potential, lookup).get(queue) - ref = eval_func_on_discr_nodes(queue, discr, func).get(queue) - return np.allclose(ref, res) - - -def drive_test_to_meshmode_interpolation_2d( - cl_ctx, queue, - degree, nel_1d, n_levels, q_order, - a=-0.5, b=0.5, seed=0): - """ - meshmode mesh control: nel_1d, degree - volumential mesh control: n_levels, q_order - """ - dim = 2 - - mesh = generate_regular_rect_mesh( - a=(a, a), - b=(b, b), - n=(nel_1d, nel_1d)) - - arr_ctx = PyOpenCLArrayContext(queue) - group_factory = PolynomialWarpAndBlendGroupFactory(order=degree) - discr = Discretization(arr_ctx, mesh, group_factory) - - bbox_fac = BoundingBoxFactory(dim=2) - boxfmm_fac = BoxFMMGeometryFactory( - cl_ctx, dim=dim, order=q_order, - nlevels=n_levels, bbox_getter=bbox_fac, - expand_to_hold_mesh=mesh, mesh_padding_factor=0.) - boxgeo = boxfmm_fac(queue) - lookup_fac = LeavesToNodesLookupBuilder( - cl_ctx, trav=boxgeo.trav, discr=discr) - lookup, evt = lookup_fac(queue) - def func(pts): - x, y = pts - return np.sin(x + y) * x + if test_case == 'exact': + return np.allclose(ref, res) - tree = boxgeo.tree.get(queue) - potential = func(np.vstack(tree.sources)) - res = interpolate_to_meshmode(queue, potential, lookup).get(queue) - - ref = eval_func_on_discr_nodes(queue, discr, func).get(queue) resid = np.linalg.norm(ref - res, ord=np.inf) return resid +# {{{ 2d tests + @pytest.mark.parametrize("params", [ [1, 8, 3, 1], [1, 8, 5, 2], [1, 8, 7, 3], [2, 16, 3, 1], [2, 32, 5, 2], [2, 64, 7, 3], @@ -243,18 +192,18 @@ def drive_test_to_meshmode_interpolation_2d( def test_from_meshmode_interpolation_2d_exact(ctx_factory, params): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) - assert drive_test_from_meshmode_exact_interpolation_2d( - cl_ctx, queue, *params) + assert drive_test_from_meshmode_interpolation( + cl_ctx, queue, 2, *params, test_case='exact') @pytest.mark.parametrize("params", [ - [1, 32, 5, 1], [2, 64, 5, 3], [3, 64, 5, 4] + [1, 128, 5, 1], [2, 64, 5, 3], [3, 64, 5, 4] ]) def test_from_meshmode_interpolation_2d_nonexact(ctx_factory, params): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) - assert drive_test_from_meshmode_interpolation_2d( - cl_ctx, queue, *params) < 1e-3 + assert drive_test_from_meshmode_interpolation( + cl_ctx, queue, 2, *params, test_case='non-exact') < 1e-3 @pytest.mark.parametrize("params", [ @@ -264,8 +213,8 @@ def test_from_meshmode_interpolation_2d_nonexact(ctx_factory, params): def test_to_meshmode_interpolation_2d_exact(ctx_factory, params): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) - assert drive_test_to_meshmode_exact_interpolation_2d( - cl_ctx, queue, *params) + assert drive_test_to_meshmode_interpolation( + cl_ctx, queue, 2, *params, test_case='exact') @pytest.mark.parametrize("params", [ @@ -274,11 +223,62 @@ def test_to_meshmode_interpolation_2d_exact(ctx_factory, params): def test_to_meshmode_interpolation_2d_nonexact(ctx_factory, params): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) - assert drive_test_to_meshmode_interpolation_2d( - cl_ctx, queue, *params) < 1e-3 + assert drive_test_to_meshmode_interpolation( + cl_ctx, queue, 2, *params, test_case='non-exact') < 1e-3 + +# }}} End 2d tests + + +# {{{ 3d tests + +@pytest.mark.parametrize("params", [ + [1, 3, 3, 1], [2, 2, 4, 2], [3, 3, 4, 3], + [4, 4, 3, 1], [5, 2, 2, 2], [8, 4, 3, 3], + ]) +def test_from_meshmode_interpolation_3d_exact(ctx_factory, params): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + assert drive_test_from_meshmode_interpolation( + cl_ctx, queue, 3, *params, test_case='exact') + + +@pytest.mark.parametrize("params", [ + [6, 4, 5, 1], [7, 3, 4, 3], [8, 2, 3, 4] + ]) +def test_from_meshmode_interpolation_3d_nonexact(ctx_factory, params): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + assert drive_test_from_meshmode_interpolation( + cl_ctx, queue, 3, *params, test_case='non-exact') < 1e-3 + + +@pytest.mark.parametrize("params", [ + [1, 5, 3, 2], [2, 7, 3, 3], [3, 7, 3, 4], + [4, 8, 3, 5], [8, 10, 3, 9], [9, 7, 2, 10], + ]) +def test_to_meshmode_interpolation_3d_exact(ctx_factory, params): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + assert drive_test_to_meshmode_interpolation( + cl_ctx, queue, 3, *params, test_case='exact') + + +@pytest.mark.parametrize("params", [ + [2, 5, 4, 4], [3, 7, 5, 3], [4, 7, 3, 5], + ]) +def test_to_meshmode_interpolation_3d_nonexact(ctx_factory, params): + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + assert drive_test_to_meshmode_interpolation( + cl_ctx, queue, 3, *params, test_case='non-exact') < 1e-3 + +# }}} End 3d tests if __name__ == '__main__': cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - drive_test_from_meshmode_interpolation_2d(cl_ctx, queue, 1, 8, 5, 3) + resid = drive_test_to_meshmode_interpolation( + cl_ctx, queue, + dim=3, degree=9, nel_1d=7, n_levels=2, q_order=10, + test_case="exact") diff --git a/volumential/geometry.py b/volumential/geometry.py index e454ef8596d13b293daffc175dfec3d6b94fbb96..0250c9195c4725fb7bb99d79bbb6336d653873bc 100644 --- a/volumential/geometry.py +++ b/volumential/geometry.py @@ -275,7 +275,7 @@ class BoxFMMGeometryFactory(): tree, _ = tb(queue, particles=q_points, targets=q_points, bbox=self._bbox, max_particles_in_box=( - self.n_q_points_per_cell * (2**self.dim) * (2**self.dim) + (self.n_q_points_per_cell**self.dim) * (2**self.dim) - 1), kind="adaptive-level-restricted") diff --git a/volumential/interpolation.py b/volumential/interpolation.py index 4a737cecc5c5933d1af0c8056c0cb566e7bcc952..85191f93c1dd0fe7c0e7ed7988b28dd790eae994 100644 --- a/volumential/interpolation.py +++ b/volumential/interpolation.py @@ -22,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +import itertools import numpy as np import pyopencl as cl import loopy as lp @@ -172,6 +173,116 @@ class ElementsToSourcesLookupBuilder: # {{{ kernel generation + @memoize_method + def codegen_get_dimension_specific_snippets(self): + """Dimension-dependent code loopy instructions. + """ + import sympy as sp + axis_names = ["x", "y", "z"] + axis_names = axis_names[:self.dim] + + # tolerance + tol = -1e-12 + + def make_sympy_vec(comp_names): + comps = [] + for cn in comp_names: + comps.append(sp.var(cn)) + return sp.Matrix(comps) + + def get_simplex_measure(vtx_names): + mat0 = sp.ones(self.dim + 1) + for iv, v in enumerate(vtx_names): + vtx = sp.Matrix([sp.var(f"{v}{comp}") for comp in axis_names]) + mat0[iv, :-1] = vtx.T + return str(mat0.det()) + + if self.dim == 2: + + # {{{ 2d + + code_get_simplex = \ + """ + <> Ax = mesh_vertices_0[mesh_vertex_indices[iel, 0]] + <> Ay = mesh_vertices_1[mesh_vertex_indices[iel, 0]] + <> Bx = mesh_vertices_0[mesh_vertex_indices[iel, 1]] + <> By = mesh_vertices_1[mesh_vertex_indices[iel, 1]] + <> Cx = mesh_vertices_0[mesh_vertex_indices[iel, 2]] + <> Cy = mesh_vertices_1[mesh_vertex_indices[iel, 2]] + """ + code_get_point = \ + """ + <> Px = source_points_0[source_id] + <> Py = source_points_1[source_id] + """ + # simplex measures + code_s0 = get_simplex_measure(["P", "B", "C"]) + code_s1 = get_simplex_measure(["A", "P", "C"]) + code_s2 = get_simplex_measure(["A", "B", "P"]) + code_compute_simplex_measures = \ + f""" + <> s0 = {code_s0} + <> s1 = {code_s1} + <> s2 = {code_s2} + """ + code_measures_have_common_sign = " and ".join([ + f"s{c1} * s{c2} >= {tol}" + for c1, c2 in itertools.combinations(["0", "1"], 2)]) + + # }}} End 2d + + elif self.dim == 3: + + # {{{ 3d + + code_get_simplex = \ + """ + <> Ax = mesh_vertices_0[mesh_vertex_indices[iel, 0]] + <> Ay = mesh_vertices_1[mesh_vertex_indices[iel, 0]] + <> Az = mesh_vertices_2[mesh_vertex_indices[iel, 0]] + <> Bx = mesh_vertices_0[mesh_vertex_indices[iel, 1]] + <> By = mesh_vertices_1[mesh_vertex_indices[iel, 1]] + <> Bz = mesh_vertices_2[mesh_vertex_indices[iel, 1]] + <> Cx = mesh_vertices_0[mesh_vertex_indices[iel, 2]] + <> Cy = mesh_vertices_1[mesh_vertex_indices[iel, 2]] + <> Cz = mesh_vertices_2[mesh_vertex_indices[iel, 2]] + <> Dx = mesh_vertices_0[mesh_vertex_indices[iel, 3]] + <> Dy = mesh_vertices_1[mesh_vertex_indices[iel, 3]] + <> Dz = mesh_vertices_2[mesh_vertex_indices[iel, 3]] + """ + code_get_point = \ + """ + <> Px = source_points_0[source_id] + <> Py = source_points_1[source_id] + <> Pz = source_points_2[source_id] + """ + # simplex measures + code_s0 = get_simplex_measure(["P", "B", "C", "D"]) + code_s1 = get_simplex_measure(["A", "P", "C", "D"]) + code_s2 = get_simplex_measure(["A", "B", "P", "D"]) + code_s3 = get_simplex_measure(["A", "B", "C", "P"]) + code_compute_simplex_measures = \ + f""" + <> s0 = {code_s0} + <> s1 = {code_s1} + <> s2 = {code_s2} + <> s3 = {code_s3} + """ + code_measures_have_common_sign = " and ".join([ + f"s{c1} * s{c2} >= {tol}" + for c1, c2 in itertools.combinations(["0", "1", "2"], 2)]) + + # }}} End 3d + + else: + raise NotImplementedError() + + return {"code_get_simplex": code_get_simplex, + "code_get_point": code_get_point, + "code_compute_simplex_measures": code_compute_simplex_measures, + "code_measures_have_common_sign": code_measures_have_common_sign, + } + @memoize_method def get_simplex_lookup_kernel(self): """Returns a loopy kernel that computes a potential vector @@ -185,9 +296,7 @@ class ElementsToSourcesLookupBuilder: """ logger.debug("start building elements-to-sources lookup kernel") - if self.dim != 2: - raise NotImplementedError() - + snippets = self.codegen_get_dimension_specific_snippets() loopy_knl = lp.make_kernel( ["{ [ iel ]: 0 <= iel < nelements }", "{ [ ineighbor ]: nearby_leaves_beg <= ineighbor < nearby_leaves_end }", @@ -198,12 +307,7 @@ class ElementsToSourcesLookupBuilder: <> nearby_leaves_beg = leaves_near_ball_starts[iel] <> nearby_leaves_end = leaves_near_ball_starts[iel + 1] - <> Ax = mesh_vertices_0[mesh_vertex_indices[iel, 0]] - <> Ay = mesh_vertices_1[mesh_vertex_indices[iel, 0]] - <> Bx = mesh_vertices_0[mesh_vertex_indices[iel, 1]] - <> By = mesh_vertices_1[mesh_vertex_indices[iel, 1]] - <> Cx = mesh_vertices_0[mesh_vertex_indices[iel, 2]] - <> Cy = mesh_vertices_1[mesh_vertex_indices[iel, 2]] + {code_get_simplex} for ineighbor <> ileaf = leaves_near_ball_lists[ineighbor] @@ -212,27 +316,24 @@ class ElementsToSourcesLookupBuilder: for isrc <> source_id = box_source_beg + isrc - <> Px = source_points_0[source_id] - <> Py = source_points_1[source_id] - <> s0 = (Px - Ax) * (By - Ay) - (Py - Ay) * (Bx - Ax) - <> s1 = (Px - Bx) * (Cy - By) - (Py - By) * (Cx - Bx) - <> s2 = (Px - Cx) * (Ay - Cy) - (Py - Cy) * (Ax - Cx) + {code_get_point} + {code_compute_simplex_measures} result[source_id] = if( - s0 * s1 >= -1e-12 and s1 * s2 >= -1e-12, + {code_measures_have_common_sign}, iel, - result[source_id]) + result[source_id]) {{atomic}} end end end - """], + """.format(**snippets)], [lp.ValueArg("nelements, dim, nboxes, nsources", np.int32), lp.GlobalArg("mesh_vertex_indices", np.int32, "nelements, dim+1"), lp.GlobalArg("box_source_starts", np.int32, "nboxes"), lp.GlobalArg("box_source_counts_cumul", np.int32, "nboxes"), lp.GlobalArg("leaves_near_ball_lists", np.int32, None), - lp.GlobalArg("result", np.int32, "nsources"), + lp.GlobalArg("result", np.int32, "nsources", for_atomic=True), "..."], name="build_sources_in_simplex_lookup", lang_version=(2018, 2), @@ -251,7 +352,7 @@ class ElementsToSourcesLookupBuilder: raise NotImplementedError("Mixed elements not supported") melgrp = mesh.groups[0] ball_centers_host = (np.max(melgrp.nodes, axis=2) - + np.min(melgrp.nodes, axis=2)) / 2 + + np.min(melgrp.nodes, axis=2)) / 2 ball_radii_host = np.max( np.max(melgrp.nodes, axis=2) - np.min(melgrp.nodes, axis=2), axis=0) / 2 @@ -288,24 +389,32 @@ class ElementsToSourcesLookupBuilder: cl.array.to_device(queue, verts) for verts in self.discr.mesh.vertices]) + mesh_vertices_kwargs = { + f"mesh_vertices_{iaxis}": vertices_dev[iaxis] + for iaxis in range(self.dim)} + + source_points_kwargs = { + f"source_points_{iaxis}": self.tree.sources[iaxis] + for iaxis in range(self.dim)} + evt, res = element_lookup_kernel( queue, dim=self.dim, nboxes=self.tree.nboxes, nelements=self.discr.mesh.nelements, nsources=self.tree.nsources, - result=cl.array.zeros(queue, self.tree.nsources, dtype=np.int32), + result=cl.array.zeros(queue, self.tree.nsources, dtype=np.int32) - 1, mesh_vertex_indices=self.discr.mesh.groups[0].vertex_indices, - mesh_vertices_0=vertices_dev[0], - mesh_vertices_1=vertices_dev[1], - source_points_0=self.tree.sources[0], - source_points_1=self.tree.sources[1], box_source_starts=self.tree.box_source_starts, box_source_counts_cumul=self.tree.box_source_counts_cumul, leaves_near_ball_starts=balls_to_leaves_lookup.leaves_near_ball_starts, leaves_near_ball_lists=balls_to_leaves_lookup.leaves_near_ball_lists, - wait_for=wait_for) + wait_for=wait_for, **mesh_vertices_kwargs, **source_points_kwargs) source_to_element_lookup, = res + wait_for = [evt] + # elements = source_to_element_lookup.get() + # for idx in [362, 365, 874, 877, 1386, 1389, 1898, 1901]) + # ----------------------------------------------------------------- # Invert the source-to-element lookup by a key-value sort @@ -602,7 +711,6 @@ def interpolate_to_meshmode(queue, potential, leaves_to_nodes_lookup, # infer q_order from tree pts_per_box = tree.ntargets // traversal.ntarget_boxes - print(pts_per_box, traversal.ntarget_boxes, tree.ntargets) assert pts_per_box * traversal.ntarget_boxes == tree.ntargets # allow for +/- 0.25 floating point error diff --git a/volumential/meshgen.py b/volumential/meshgen.py index 0b835af15a75a3438b607667d059b240263d129d..cc936b2367b10446d2125dcce43e74c954252c61 100644 --- a/volumential/meshgen.py +++ b/volumential/meshgen.py @@ -178,7 +178,7 @@ try: except ImportError as e: logger.debug(repr(e)) - logger.warn("Meshgen via Deal.II is not present or unusable.") + logger.warning("Meshgen via Deal.II is not present or unusable.") try: logger.info("Trying out BoxTree.TreeInteractiveBuild interface.") @@ -189,7 +189,7 @@ except ImportError as e: except ImportError as ee: logger.debug(repr(ee)) - logger.warn("Meshgen via BoxTree is not present or unusable.") + logger.warning("Meshgen via BoxTree is not present or unusable.") raise RuntimeError("Cannot find a usable Meshgen implementation.") else: diff --git a/volumential/volume_fmm.py b/volumential/volume_fmm.py index 409f904a84b0747be234d300378fbe1b38668c63..cb2848191e0b045a12cb07ef0b34d9095350f089 100644 --- a/volumential/volume_fmm.py +++ b/volumential/volume_fmm.py @@ -518,7 +518,7 @@ def interpolate_volume_potential(target_points, traversal, wrangler, potential, end # Count how many times the potential is computed - multiplicity[target_point_id] = multiplicity[target_point_id] + 1 + multiplicity[target_point_id] = multiplicity[target_point_id] + 1 {atomic} # Map target point to template box for iaxis @@ -575,7 +575,7 @@ def interpolate_volume_potential(target_points, traversal, wrangler, potential, p_out[target_point_id] = p_out[target_point_id] + sum(mid, mode_coeff * prod_mode_val - ) {id=p_out,dep=pmod} + ) {id=p_out,dep=pmod,atomic} end @@ -592,11 +592,14 @@ def interpolate_volume_potential(target_points, traversal, wrangler, potential, # loopy.TemporaryVariable("diff", dtype, "dim, q_order"), loopy.GlobalArg("box_centers", None, "dim, aligned_nboxes"), loopy.GlobalArg("balls_near_box_lists", None, None), + loopy.GlobalArg("multiplicity", None, None, for_atomic=True), + loopy.GlobalArg("p_out", None, None, for_atomic=True), loopy.ValueArg("aligned_nboxes", np.int32), loopy.ValueArg("dim", np.int32), loopy.ValueArg("q_order", np.int32), "...", ], + lang_version=(2018, 2), ) # }}} End loopy kernel for interpolation