From 0b9c3ec3eb9722f955f688d92a2bfb0c478cedbb Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 15:34:55 -0500 Subject: [PATCH 01/11] 3D point-in-simplex tests --- volumential/interpolation.py | 144 ++++++++++++++++++++++++++++++----- 1 file changed, 123 insertions(+), 21 deletions(-) diff --git a/volumential/interpolation.py b/volumential/interpolation.py index 4a737ce..9508d55 100644 --- a/volumential/interpolation.py +++ b/volumential/interpolation.py @@ -172,6 +172,114 @@ 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 = \ + f"s0 * s1 >= {tol} and s1 * s2 >= {tol}" + + # }}} 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]] + <> Ay = 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]] + <> By = 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]] + <> Cy = 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 = \ + f"s0 * s1 >= {tol} and s1 * s2 >= {tol} and s2 * s3 >= {tol}" + + # }}} 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 +293,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 +304,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,21 +313,18 @@ 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]) 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"), @@ -288,20 +386,24 @@ 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), 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] -- GitLab From aa21965fb8b369c1f1a3ae7a5e31f64d33f55fc2 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 15:53:15 -0500 Subject: [PATCH 02/11] Refactor interpolation tests --- test/test_interpolation.py | 127 ++++++++++--------------------------- 1 file changed, 32 insertions(+), 95 deletions(-) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index d38cd85..2d8c25b 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -70,51 +70,10 @@ def eval_func_on_discr_nodes(queue, discr, func): # }}} End test data -def drive_test_from_meshmode_exact_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) - - # algebraically exact interpolation - func = random_polynomial_func(dim, degree, seed) - - 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)) - - 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): + a=-0.5, b=0.5, seed=0, test_case='exact'): """ meshmode mesh control: nel_1d, degree volumential mesh control: n_levels, q_order @@ -140,9 +99,15 @@ def drive_test_from_meshmode_interpolation_2d( cl_ctx, tree=boxgeo.tree, discr=discr) lookup, evt = lookup_fac(queue) - def func(pts): - x, y = pts - return np.sin(x + y) * x + if test_case == 'exact': + # algebraically exact interpolation + func = random_polynomial_func(dim, degree, seed) + elif test_case == 'non-exact': + def func(pts): + x, y = pts + return np.sin(x + y) * x + else: + raise ValueError() dof_vec = eval_func_on_discr_nodes(queue, discr, func) res = interpolate_from_meshmode(queue, dof_vec, lookup).get(queue) @@ -150,14 +115,17 @@ def drive_test_from_meshmode_interpolation_2d( 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( +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): + a=-0.5, b=0.5, seed=0, test_case='exact'): """ meshmode mesh control: nel_1d, degree volumential mesh control: n_levels, q_order @@ -183,55 +151,24 @@ 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': + def func(pts): + x, y = pts + return np.sin(x + y) * x + 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) + if test_case == 'exact': + return np.allclose(ref, res) - 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 - - 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 @@ -243,8 +180,8 @@ 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_2d( + cl_ctx, queue, *params, test_case='exact') @pytest.mark.parametrize("params", [ @@ -254,7 +191,7 @@ 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 + cl_ctx, queue, *params, test_case='non-exact') < 1e-3 @pytest.mark.parametrize("params", [ @@ -264,8 +201,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_2d( + cl_ctx, queue, *params, test_case='exact') @pytest.mark.parametrize("params", [ @@ -275,7 +212,7 @@ 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 + cl_ctx, queue, *params, test_case='non-exact') < 1e-3 if __name__ == '__main__': -- GitLab From 8d5d66af3f83867be998cffeeb147826e390ff6d Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 16:02:54 -0500 Subject: [PATCH 03/11] More refactoring --- test/test_interpolation.py | 70 ++++++++++++++++++++++---------------- 1 file changed, 40 insertions(+), 30 deletions(-) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index 2d8c25b..76784d8 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_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, 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, @@ -103,9 +101,16 @@ def drive_test_from_meshmode_interpolation_2d( # algebraically exact interpolation func = random_polynomial_func(dim, degree, seed) elif test_case == 'non-exact': - def func(pts): - x, y = pts - return np.sin(x + y) * x + 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() @@ -122,26 +127,24 @@ def drive_test_from_meshmode_interpolation_2d( return resid -def drive_test_to_meshmode_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, 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, @@ -155,9 +158,16 @@ def drive_test_to_meshmode_interpolation_2d( # algebraically exact interpolation func = random_polynomial_func(dim, degree, seed) elif test_case == 'non-exact': - def func(pts): - x, y = pts - return np.sin(x + y) * x + 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() @@ -180,8 +190,8 @@ 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_interpolation_2d( - cl_ctx, queue, *params, test_case='exact') + assert drive_test_from_meshmode_interpolation( + cl_ctx, queue, 2, *params, test_case='exact') @pytest.mark.parametrize("params", [ @@ -190,8 +200,8 @@ def test_from_meshmode_interpolation_2d_exact(ctx_factory, params): 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, test_case='non-exact') < 1e-3 + assert drive_test_from_meshmode_interpolation( + cl_ctx, queue, 2, *params, test_case='non-exact') < 1e-3 @pytest.mark.parametrize("params", [ @@ -201,8 +211,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_interpolation_2d( - cl_ctx, queue, *params, test_case='exact') + assert drive_test_to_meshmode_interpolation( + cl_ctx, queue, 2, *params, test_case='exact') @pytest.mark.parametrize("params", [ @@ -211,8 +221,8 @@ 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, test_case='non-exact') < 1e-3 + assert drive_test_to_meshmode_interpolation( + cl_ctx, queue, 2, *params, test_case='non-exact') < 1e-3 if __name__ == '__main__': -- GitLab From 6b7344cbde03c16765fda384828ad106b2beadfb Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 17:02:34 -0500 Subject: [PATCH 04/11] Various fixes --- contrib/meshgen11_dealii/CMakeLists.txt | 5 ++-- test/test_interpolation.py | 37 ++++++++++++++++++++++++- volumential/interpolation.py | 7 ++--- volumential/meshgen.py | 4 +-- volumential/volume_fmm.py | 4 ++- 5 files changed, 46 insertions(+), 11 deletions(-) diff --git a/contrib/meshgen11_dealii/CMakeLists.txt b/contrib/meshgen11_dealii/CMakeLists.txt index b7c9ea9..3073e50 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 76784d8..6601a02 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -183,6 +183,8 @@ def drive_test_to_meshmode_interpolation( 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], @@ -224,8 +226,41 @@ def test_to_meshmode_interpolation_2d_nonexact(ctx_factory, params): 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, 5, 3, 2], [2, 7, 3, 3], [3, 7, 3, 4], + ]) +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) + + # FIXME: the following set of parameters produce random failures + for iter in range(10): + resid = drive_test_to_meshmode_interpolation( + cl_ctx, queue, + dim=3, degree=1, nel_1d=5, n_levels=4, q_order=2, + test_case="exact") + print(resid) diff --git a/volumential/interpolation.py b/volumential/interpolation.py index 9508d55..e3571ba 100644 --- a/volumential/interpolation.py +++ b/volumential/interpolation.py @@ -237,13 +237,13 @@ class ElementsToSourcesLookupBuilder: """ <> Ax = mesh_vertices_0[mesh_vertex_indices[iel, 0]] <> Ay = mesh_vertices_1[mesh_vertex_indices[iel, 0]] - <> Ay = mesh_vertices_2[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]] - <> By = mesh_vertices_2[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]] - <> Cy = mesh_vertices_2[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]] @@ -704,7 +704,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 0b835af..cc936b2 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 409f904..40e66f4 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 @@ -592,11 +592,13 @@ 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.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 -- GitLab From 01de2fdec17242c76b169efb8cd9ab3fb0378111 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 17:52:02 -0500 Subject: [PATCH 05/11] Use atomic operations for interpolation --- test/test_interpolation.py | 7 ++++++- volumential/volume_fmm.py | 13 ++++++++++++- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index 6601a02..0fcd266 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -177,6 +177,10 @@ def drive_test_to_meshmode_interpolation( ref = eval_func_on_discr_nodes(queue, discr, func).get(queue) if test_case == 'exact': + if 1: + print(np.linalg.norm(ref - res, ord=np.inf), '&') + print(np.where(np.abs(ref - res) > 1e-10), '&') + return np.allclose(ref, res) resid = np.linalg.norm(ref - res, ord=np.inf) @@ -259,8 +263,9 @@ if __name__ == '__main__': # FIXME: the following set of parameters produce random failures for iter in range(10): + print("-----") resid = drive_test_to_meshmode_interpolation( cl_ctx, queue, dim=3, degree=1, nel_1d=5, n_levels=4, q_order=2, test_case="exact") - print(resid) + print(resid, '&') diff --git a/volumential/volume_fmm.py b/volumential/volume_fmm.py index 40e66f4..6de305b 100644 --- a/volumential/volume_fmm.py +++ b/volumential/volume_fmm.py @@ -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 @@ -593,6 +593,7 @@ def interpolate_volume_potential(target_points, traversal, wrangler, potential, 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), @@ -626,6 +627,12 @@ def interpolate_volume_potential(target_points, traversal, wrangler, potential, # fetching from user_source_ids converts potential to tree order user_mode_ids = tree.user_source_ids + if 0: + import sys + np.set_printoptions(threshold=sys.maxsize) + print(balls_near_box_lists.with_queue(queue).get(), '&') + print(balls_near_box_starts.with_queue(queue).get(), '&') + lpknl = loopy.set_options(lpknl, return_dict=True) lpknl = loopy.fix_parameters(lpknl, dim=int(dim), q_order=int(q_order)) lpknl = loopy.split_iname(lpknl, "tbox", 128, outer_tag="g.0", inner_tag="l.0") @@ -653,6 +660,10 @@ def interpolate_volume_potential(target_points, traversal, wrangler, potential, pout.add_event(evt) multiplicity.add_event(evt) + if 0: + print(pout.get(), '&') + print(multiplicity.get(), '&') + return pout / multiplicity # }}} End free form interpolation of potentials -- GitLab From aead16ee4138857425758f6b83b13470cabb9c19 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 17:53:01 -0500 Subject: [PATCH 06/11] Clean up --- test/test_interpolation.py | 18 +++++------------- volumential/volume_fmm.py | 10 ---------- 2 files changed, 5 insertions(+), 23 deletions(-) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index 0fcd266..8903b91 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -177,10 +177,6 @@ def drive_test_to_meshmode_interpolation( ref = eval_func_on_discr_nodes(queue, discr, func).get(queue) if test_case == 'exact': - if 1: - print(np.linalg.norm(ref - res, ord=np.inf), '&') - print(np.where(np.abs(ref - res) > 1e-10), '&') - return np.allclose(ref, res) resid = np.linalg.norm(ref - res, ord=np.inf) @@ -260,12 +256,8 @@ def test_to_meshmode_interpolation_3d_nonexact(ctx_factory, params): if __name__ == '__main__': cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - - # FIXME: the following set of parameters produce random failures - for iter in range(10): - print("-----") - resid = drive_test_to_meshmode_interpolation( - cl_ctx, queue, - dim=3, degree=1, nel_1d=5, n_levels=4, q_order=2, - test_case="exact") - print(resid, '&') + resid = drive_test_to_meshmode_interpolation( + cl_ctx, queue, + dim=3, degree=1, nel_1d=5, n_levels=4, q_order=2, + test_case="exact") + print(resid) diff --git a/volumential/volume_fmm.py b/volumential/volume_fmm.py index 6de305b..cb28481 100644 --- a/volumential/volume_fmm.py +++ b/volumential/volume_fmm.py @@ -627,12 +627,6 @@ def interpolate_volume_potential(target_points, traversal, wrangler, potential, # fetching from user_source_ids converts potential to tree order user_mode_ids = tree.user_source_ids - if 0: - import sys - np.set_printoptions(threshold=sys.maxsize) - print(balls_near_box_lists.with_queue(queue).get(), '&') - print(balls_near_box_starts.with_queue(queue).get(), '&') - lpknl = loopy.set_options(lpknl, return_dict=True) lpknl = loopy.fix_parameters(lpknl, dim=int(dim), q_order=int(q_order)) lpknl = loopy.split_iname(lpknl, "tbox", 128, outer_tag="g.0", inner_tag="l.0") @@ -660,10 +654,6 @@ def interpolate_volume_potential(target_points, traversal, wrangler, potential, pout.add_event(evt) multiplicity.add_event(evt) - if 0: - print(pout.get(), '&') - print(multiplicity.get(), '&') - return pout / multiplicity # }}} End free form interpolation of potentials -- GitLab From 425299ae060c6f4cc4af43baccc87f485056c22f Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 19:47:19 -0500 Subject: [PATCH 07/11] Fix point-in-simplex test when the point lies exactly on a hyperplane --- test/test_interpolation.py | 7 ++-- volumential/interpolation.py | 74 ++++++++++++++++++++++++++++++++---- 2 files changed, 69 insertions(+), 12 deletions(-) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index 8903b91..67c00d5 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -197,7 +197,7 @@ def test_from_meshmode_interpolation_2d_exact(ctx_factory, params): @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() @@ -256,8 +256,7 @@ def test_to_meshmode_interpolation_3d_nonexact(ctx_factory, params): if __name__ == '__main__': cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - resid = drive_test_to_meshmode_interpolation( + resid = drive_test_from_meshmode_interpolation( cl_ctx, queue, - dim=3, degree=1, nel_1d=5, n_levels=4, q_order=2, + dim=3, degree=1, nel_1d=3, n_levels=4, q_order=2, test_case="exact") - print(resid) diff --git a/volumential/interpolation.py b/volumential/interpolation.py index e3571ba..b56f37e 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 @@ -224,8 +225,9 @@ class ElementsToSourcesLookupBuilder: <> s1 = {code_s1} <> s2 = {code_s2} """ - code_measures_have_common_sign = \ - f"s0 * s1 >= {tol} and s1 * s2 >= {tol}" + code_measures_have_common_sign = " and ".join([ + f"s{c1} * s{c2} >= {tol}" + for c1, c2 in itertools.combinations(["0", "1"], 2)]) # }}} End 2d @@ -266,8 +268,9 @@ class ElementsToSourcesLookupBuilder: <> s2 = {code_s2} <> s3 = {code_s3} """ - code_measures_have_common_sign = \ - f"s0 * s1 >= {tol} and s1 * s2 >= {tol} and s2 * s3 >= {tol}" + 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 @@ -320,7 +323,7 @@ class ElementsToSourcesLookupBuilder: result[source_id] = if( {code_measures_have_common_sign}, iel, - result[source_id]) + result[source_id]) {{atomic}} end end end @@ -330,7 +333,7 @@ class ElementsToSourcesLookupBuilder: 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), @@ -349,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 @@ -397,7 +400,7 @@ class ElementsToSourcesLookupBuilder: 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, box_source_starts=self.tree.box_source_starts, box_source_counts_cumul=self.tree.box_source_counts_cumul, @@ -406,8 +409,63 @@ class ElementsToSourcesLookupBuilder: wait_for=wait_for, **mesh_vertices_kwargs, **source_points_kwargs) source_to_element_lookup, = res + + if 0: + from mpl_toolkits.mplot3d import Axes3D + import matplotlib.pyplot as plt + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + if 0: + ax.scatter( + self.tree.sources[0].get(queue), + self.tree.sources[1].get(queue), + self.tree.sources[2].get(queue), + c=source_to_element_lookup.get(queue), + cmap='tab20c', alpha=0.5) + + highlight = [362, 365, 874, 877, 1386, 1389, 1898, 1901] + for h in highlight[:1]: + print("=====") + src0 = self.tree.sources[0].get(queue)[h] + src1 = self.tree.sources[1].get(queue)[h] + src2 = self.tree.sources[2].get(queue)[h] + print(f"A bad source: ({src0}, {src1}, {src2})") + + e = source_to_element_lookup.get(queue)[h] + tet = self.discr.mesh.groups[0].nodes[:, e, :] + tetvs = [0, 1, 2, 0, 3, 1, 3, 2] + print(f"It is said to be in element {e}: {tet}") + + ax.plot(tet[0, tetvs], tet[1, tetvs], tet[2, tetvs]) + ax.scatter( + self.tree.sources[0].get(queue)[h], + self.tree.sources[1].get(queue)[h], + self.tree.sources[2].get(queue)[h], + # c=source_to_element_lookup.get(queue)[highlight], + cmap='tab20c', alpha=1, s=50) + + + # plt.show() + Ax, Ay, Az = tet[:, 0] + Bx, By, Bz = tet[:, 1] + Cx, Cy, Cz = tet[:, 2] + Dx, Dy, Dz = tet[:, 3] + Px, Py, Pz = (src0, src1, src2) + + s3 = Ax * By * Cz + -1.0 * Ax * By * Pz + -1.0 * Ax * Bz * Cy + Ax * Bz * Py + Ax * Cy * Pz + -1.0 * Ax * Cz * Py + -1.0 * Ay * Bx * Cz + Ay * Bx * Pz + Ay * Bz * Cx + -1.0 * Ay * Bz * Px + -1.0 * Ay * Cx * Pz + Ay * Cz * Px + Az * Bx * Cy + -1.0 * Az * Bx * Py + -1.0 * Az * By * Cx + Az * By * Px + Az * Cx * Py + -1.0 * Az * Cy * Px + -1.0 * Bx * Cy * Pz + Bx * Cz * Py + By * Cx * Pz + -1.0 * By * Cz * Px + -1.0 * Bz * Cx * Py + Bz * Cy * Px + s2 = -1.0 * Ax * By * Dz + Ax * By * Pz + Ax * Bz * Dy + -1.0 * Ax * Bz * Py + -1.0 * Ax * Dy * Pz + Ax * Dz * Py + Ay * Bx * Dz + -1.0 * Ay * Bx * Pz + -1.0 * Ay * Bz * Dx + Ay * Bz * Px + Ay * Dx * Pz + -1.0 * Ay * Dz * Px + -1.0 * Az * Bx * Dy + Az * Bx * Py + Az * By * Dx + -1.0 * Az * By * Px + -1.0 * Az * Dx * Py + Az * Dy * Px + Bx * Dy * Pz + -1.0 * Bx * Dz * Py + -1.0 * By * Dx * Pz + By * Dz * Px + Bz * Dx * Py + -1.0 * Bz * Dy * Px + s1 = Ax * Cy * Dz + -1.0 * Ax * Cy * Pz + -1.0 * Ax * Cz * Dy + Ax * Cz * Py + Ax * Dy * Pz + -1.0 * Ax * Dz * Py + -1.0 * Ay * Cx * Dz + Ay * Cx * Pz + Ay * Cz * Dx + -1.0 * Ay * Cz * Px + -1.0 * Ay * Dx * Pz + Ay * Dz * Px + Az * Cx * Dy + -1.0 * Az * Cx * Py + -1.0 * Az * Cy * Dx + Az * Cy * Px + Az * Dx * Py + -1.0 * Az * Dy * Px + -1.0 * Cx * Dy * Pz + Cx * Dz * Py + Cy * Dx * Pz + -1.0 * Cy * Dz * Px + -1.0 * Cz * Dx * Py + Cz * Dy * Px + s0 = -1.0 * Bx * Cy * Dz + Bx * Cy * Pz + Bx * Cz * Dy + -1.0 * Bx * Cz * Py + -1.0 * Bx * Dy * Pz + Bx * Dz * Py + By * Cx * Dz + -1.0 * By * Cx * Pz + -1.0 * By * Cz * Dx + By * Cz * Px + By * Dx * Pz + -1.0 * By * Dz * Px + -1.0 * Bz * Cx * Dy + Bz * Cx * Py + Bz * Cy * Dx + -1.0 * Bz * Cy * Px + -1.0 * Bz * Dx * Py + Bz * Dy * Px + Cx * Dy * Pz + -1.0 * Cx * Dz * Py + -1.0 * Cy * Dx * Pz + Cy * Dz * Px + Cz * Dx * Py + -1.0 * Cz * Dy * Px + + print(s0, s1, s2, s3) + print(s0 * s1) + + 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 -- GitLab From ea615234b3eb59f859c08c212fb07eaac9f725ca Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 19:47:51 -0500 Subject: [PATCH 08/11] Clean up --- volumential/interpolation.py | 51 ------------------------------------ 1 file changed, 51 deletions(-) diff --git a/volumential/interpolation.py b/volumential/interpolation.py index b56f37e..85191f9 100644 --- a/volumential/interpolation.py +++ b/volumential/interpolation.py @@ -410,57 +410,6 @@ class ElementsToSourcesLookupBuilder: source_to_element_lookup, = res - if 0: - from mpl_toolkits.mplot3d import Axes3D - import matplotlib.pyplot as plt - fig = plt.figure() - ax = fig.add_subplot(111, projection='3d') - if 0: - ax.scatter( - self.tree.sources[0].get(queue), - self.tree.sources[1].get(queue), - self.tree.sources[2].get(queue), - c=source_to_element_lookup.get(queue), - cmap='tab20c', alpha=0.5) - - highlight = [362, 365, 874, 877, 1386, 1389, 1898, 1901] - for h in highlight[:1]: - print("=====") - src0 = self.tree.sources[0].get(queue)[h] - src1 = self.tree.sources[1].get(queue)[h] - src2 = self.tree.sources[2].get(queue)[h] - print(f"A bad source: ({src0}, {src1}, {src2})") - - e = source_to_element_lookup.get(queue)[h] - tet = self.discr.mesh.groups[0].nodes[:, e, :] - tetvs = [0, 1, 2, 0, 3, 1, 3, 2] - print(f"It is said to be in element {e}: {tet}") - - ax.plot(tet[0, tetvs], tet[1, tetvs], tet[2, tetvs]) - ax.scatter( - self.tree.sources[0].get(queue)[h], - self.tree.sources[1].get(queue)[h], - self.tree.sources[2].get(queue)[h], - # c=source_to_element_lookup.get(queue)[highlight], - cmap='tab20c', alpha=1, s=50) - - - # plt.show() - Ax, Ay, Az = tet[:, 0] - Bx, By, Bz = tet[:, 1] - Cx, Cy, Cz = tet[:, 2] - Dx, Dy, Dz = tet[:, 3] - Px, Py, Pz = (src0, src1, src2) - - s3 = Ax * By * Cz + -1.0 * Ax * By * Pz + -1.0 * Ax * Bz * Cy + Ax * Bz * Py + Ax * Cy * Pz + -1.0 * Ax * Cz * Py + -1.0 * Ay * Bx * Cz + Ay * Bx * Pz + Ay * Bz * Cx + -1.0 * Ay * Bz * Px + -1.0 * Ay * Cx * Pz + Ay * Cz * Px + Az * Bx * Cy + -1.0 * Az * Bx * Py + -1.0 * Az * By * Cx + Az * By * Px + Az * Cx * Py + -1.0 * Az * Cy * Px + -1.0 * Bx * Cy * Pz + Bx * Cz * Py + By * Cx * Pz + -1.0 * By * Cz * Px + -1.0 * Bz * Cx * Py + Bz * Cy * Px - s2 = -1.0 * Ax * By * Dz + Ax * By * Pz + Ax * Bz * Dy + -1.0 * Ax * Bz * Py + -1.0 * Ax * Dy * Pz + Ax * Dz * Py + Ay * Bx * Dz + -1.0 * Ay * Bx * Pz + -1.0 * Ay * Bz * Dx + Ay * Bz * Px + Ay * Dx * Pz + -1.0 * Ay * Dz * Px + -1.0 * Az * Bx * Dy + Az * Bx * Py + Az * By * Dx + -1.0 * Az * By * Px + -1.0 * Az * Dx * Py + Az * Dy * Px + Bx * Dy * Pz + -1.0 * Bx * Dz * Py + -1.0 * By * Dx * Pz + By * Dz * Px + Bz * Dx * Py + -1.0 * Bz * Dy * Px - s1 = Ax * Cy * Dz + -1.0 * Ax * Cy * Pz + -1.0 * Ax * Cz * Dy + Ax * Cz * Py + Ax * Dy * Pz + -1.0 * Ax * Dz * Py + -1.0 * Ay * Cx * Dz + Ay * Cx * Pz + Ay * Cz * Dx + -1.0 * Ay * Cz * Px + -1.0 * Ay * Dx * Pz + Ay * Dz * Px + Az * Cx * Dy + -1.0 * Az * Cx * Py + -1.0 * Az * Cy * Dx + Az * Cy * Px + Az * Dx * Py + -1.0 * Az * Dy * Px + -1.0 * Cx * Dy * Pz + Cx * Dz * Py + Cy * Dx * Pz + -1.0 * Cy * Dz * Px + -1.0 * Cz * Dx * Py + Cz * Dy * Px - s0 = -1.0 * Bx * Cy * Dz + Bx * Cy * Pz + Bx * Cz * Dy + -1.0 * Bx * Cz * Py + -1.0 * Bx * Dy * Pz + Bx * Dz * Py + By * Cx * Dz + -1.0 * By * Cx * Pz + -1.0 * By * Cz * Dx + By * Cz * Px + By * Dx * Pz + -1.0 * By * Dz * Px + -1.0 * Bz * Cx * Dy + Bz * Cx * Py + Bz * Cy * Dx + -1.0 * Bz * Cy * Px + -1.0 * Bz * Dx * Py + Bz * Dy * Px + Cx * Dy * Pz + -1.0 * Cx * Dz * Py + -1.0 * Cy * Dx * Pz + Cy * Dz * Px + Cz * Dx * Py + -1.0 * Cz * Dy * Px - - print(s0, s1, s2, s3) - print(s0 * s1) - - wait_for = [evt] # elements = source_to_element_lookup.get() -- GitLab From 2fbfa784843a6af5a406ce3cc76f75b177b153ba Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 19:54:17 -0500 Subject: [PATCH 09/11] Add 3d tests --- test/test_interpolation.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index 67c00d5..de558f7 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -231,8 +231,30 @@ def test_to_meshmode_interpolation_2d_nonexact(ctx_factory, params): # {{{ 3d tests +@pytest.mark.parametrize("params", [ + [1, 8, 3, 1], [2, 7, 4, 2], [3, 3, 4, 3], + [4, 16, 3, 1], [5, 32, 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", [ + [2, 64, 5, 1], [4, 32, 4, 3], [8, 16, 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, 32, 3, 5], [8, 32, 3, 9], [9, 7, 2, 10], ]) def test_to_meshmode_interpolation_3d_exact(ctx_factory, params): cl_ctx = ctx_factory() -- GitLab From 0d5653d4139240bdc5502e752fac96b64b61d899 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 21:09:17 -0500 Subject: [PATCH 10/11] Make test sizes smaller --- test/test_interpolation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index de558f7..fd66e12 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -243,7 +243,7 @@ def test_from_meshmode_interpolation_3d_exact(ctx_factory, params): @pytest.mark.parametrize("params", [ - [2, 64, 5, 1], [4, 32, 4, 3], [8, 16, 3, 4] + [2, 12, 5, 1], [4, 15, 4, 3], [8, 16, 3, 4] ]) def test_from_meshmode_interpolation_3d_nonexact(ctx_factory, params): cl_ctx = ctx_factory() @@ -254,7 +254,7 @@ def test_from_meshmode_interpolation_3d_nonexact(ctx_factory, params): @pytest.mark.parametrize("params", [ [1, 5, 3, 2], [2, 7, 3, 3], [3, 7, 3, 4], - [4, 32, 3, 5], [8, 32, 3, 9], [9, 7, 2, 10], + [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() -- GitLab From 096cee02621e1b8c03921a350aca232a131422b7 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 30 Jul 2020 21:37:59 -0500 Subject: [PATCH 11/11] Make from-meshmode tests smaller but higher order --- test/test_interpolation.py | 10 +++++----- volumential/geometry.py | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/test/test_interpolation.py b/test/test_interpolation.py index fd66e12..9135c21 100644 --- a/test/test_interpolation.py +++ b/test/test_interpolation.py @@ -232,8 +232,8 @@ def test_to_meshmode_interpolation_2d_nonexact(ctx_factory, params): # {{{ 3d tests @pytest.mark.parametrize("params", [ - [1, 8, 3, 1], [2, 7, 4, 2], [3, 3, 4, 3], - [4, 16, 3, 1], [5, 32, 2, 2], [8, 4, 3, 3], + [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() @@ -243,7 +243,7 @@ def test_from_meshmode_interpolation_3d_exact(ctx_factory, params): @pytest.mark.parametrize("params", [ - [2, 12, 5, 1], [4, 15, 4, 3], [8, 16, 3, 4] + [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() @@ -278,7 +278,7 @@ def test_to_meshmode_interpolation_3d_nonexact(ctx_factory, params): if __name__ == '__main__': cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) - resid = drive_test_from_meshmode_interpolation( + resid = drive_test_to_meshmode_interpolation( cl_ctx, queue, - dim=3, degree=1, nel_1d=3, n_levels=4, q_order=2, + 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 e454ef8..0250c91 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") -- GitLab