diff --git a/.test-conda-env-py3-requirements.txt b/.test-conda-env-py3-requirements.txt index 7cefac6831754070d40eb8427b09f70a91b498cd..e3fac95dc0e8084cd5f7d50ed120d53a5655797a 100644 --- a/.test-conda-env-py3-requirements.txt +++ b/.test-conda-env-py3-requirements.txt @@ -1,4 +1,4 @@ -git+https://github.com/inducer/boxtree +git+https://gitlab.tiker.net/inducer/boxtree git+https://github.com/inducer/pymbolic git+https://github.com/inducer/loopy git+https://github.com/inducer/sumpy diff --git a/examples/fmm-manyarm-starfish-nan.py b/examples/fmm-manyarm-starfish-nan.py new file mode 100644 index 0000000000000000000000000000000000000000..7a8b1d3f58148c03f71782419bb792a22c4e9148 --- /dev/null +++ b/examples/fmm-manyarm-starfish-nan.py @@ -0,0 +1,109 @@ +import numpy as np +import numpy.linalg as la # noqa +import pyopencl as cl +import pyopencl.clmath # noqa + +from meshmode.mesh.generation import ( # noqa + ellipse, cloverleaf, NArmedStarfish, drop, n_gon, qbx_peanut, + make_curve_mesh, starfish) +from pytential import bind, sym, norm + +import logging +logger = logging.getLogger(__name__) + + +def get_starfish_mesh(nelements, target_order): + return make_curve_mesh( + NArmedStarfish(20, 0.8), + np.linspace(0, 1, nelements+1), + target_order) + + +WITH_EXTENTS = True +EXPANSION_STICKOUT_FACTOR = 0.5 + + +def get_green_error( + queue, mesh_getter, nelements, fmm_order, qbx_order, k=0): + + target_order = 12 + + mesh = mesh_getter(nelements, target_order) + + d = mesh.ambient_dim + + # u_sym = sym.var("u") + dn_u_sym = sym.var("dn_u") + + from sumpy.kernel import LaplaceKernel + k_sym = LaplaceKernel(d) + zero_op = ( + sym.S(k_sym, dn_u_sym, qbx_forced_limit=-1) + # - sym.D(k_sym, u_sym, qbx_forced_limit="avg") + # - 0.5*u_sym + ) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory) + pre_density_discr = Discretization( + queue.context, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + from pytential.qbx import QBXLayerPotentialSource + lpot_source = QBXLayerPotentialSource( + pre_density_discr, 4*target_order, + qbx_order, fmm_order=fmm_order, + _expansions_in_tree_have_extent=WITH_EXTENTS, + _expansion_stick_out_factor=EXPANSION_STICKOUT_FACTOR, + ) + + lpot_source, _ = lpot_source.with_refinement() + + density_discr = lpot_source.density_discr + + # {{{ compute values of a solution to the PDE + + nodes_host = density_discr.nodes().get(queue) + normal = bind(density_discr, sym.normal(d))(queue).as_vector(np.object) + normal_host = [normal[j].get() for j in range(d)] + + center = np.array([3, 1, 2])[:d] + diff = nodes_host - center[:, np.newaxis] + dist_squared = np.sum(diff**2, axis=0) + dist = np.sqrt(dist_squared) + if d == 2: + u = np.log(dist) + grad_u = diff/dist_squared + elif d == 3: + u = 1/dist + grad_u = -diff/dist**3 + else: + assert False + + dn_u = 0 + for i in range(d): + dn_u = dn_u + normal_host[i]*grad_u[i] + + # }}} + + u_dev = cl.array.to_device(queue, u) + dn_u_dev = cl.array.to_device(queue, dn_u) + + bound_op = bind(lpot_source, zero_op) + error = bound_op( + queue, u=u_dev, dn_u=dn_u_dev, k=k) + + print(len(error), np.where(np.isnan(error.get()))) + return norm(density_discr, queue, error) + + +def main(): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + get_green_error(queue, get_starfish_mesh, 500, 20, 2, k=0) + + +if __name__ == "__main__": + main() diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 2aae56d9e16029f85e1c60b172ff188ca7193b20..ee716f53773e5372372931891268da839241fb20 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -75,8 +75,10 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): # FIXME default debug=False once everything has matured debug=True, _refined_for_global_qbx=False, - _expansions_in_tree_have_extent=False, - _expansion_stick_out_factor=0, + _expansions_in_tree_have_extent=True, + _expansion_stick_out_factor=0.5, + _well_sep_is_n_away=2, + _max_leaf_refine_weight=32, geometry_data_inspector=None, fmm_backend="sumpy", target_stick_out_factor=_not_provided): @@ -158,8 +160,15 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): self._expansions_in_tree_have_extent = \ _expansions_in_tree_have_extent self._expansion_stick_out_factor = _expansion_stick_out_factor + self._well_sep_is_n_away = _well_sep_is_n_away + self._max_leaf_refine_weight = _max_leaf_refine_weight self.geometry_data_inspector = geometry_data_inspector + # /!\ *All* parameters set here must also be set by copy() below, + # otherwise they will be reset to their default values behind your + # back if the layer potential source is ever copied. (such as + # during refinement) + def copy( self, density_discr=None, @@ -231,9 +240,12 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): _expansion_stick_out_factor if _expansion_stick_out_factor is not _not_provided else self._expansion_stick_out_factor), + _well_sep_is_n_away=self._well_sep_is_n_away, + _max_leaf_refine_weight=self._max_leaf_refine_weight, geometry_data_inspector=( geometry_data_inspector or self.geometry_data_inspector), - fmm_backend=self.fmm_backend) + fmm_backend=self.fmm_backend, + ) # }}} @@ -468,7 +480,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): def qbx_fmm_code_getter(self): from pytential.qbx.geometry import QBXFMMGeometryCodeGetter return QBXFMMGeometryCodeGetter(self.cl_context, self.ambient_dim, - debug=self.debug) + debug=self.debug, _well_sep_is_n_away=self._well_sep_is_n_away) # {{{ fmm-based execution diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index dc88d64b277e700ee82df42c41e7673b866eb40c..d03affe3504b0c7be81f9f739c0ce8c393c49f53 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -237,7 +237,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), wait_for = multipole_exps.events - for isrc_level, ssn in enumerate(traversal.sep_smaller_by_level): + for isrc_level, ssn in enumerate(traversal.from_sep_smaller_by_level): m2qbxl = self.code.m2qbxl( self.level_orders[isrc_level], self.qbx_order) @@ -411,8 +411,8 @@ def drive_fmm(expansion_wrangler, src_weights): local_exps = wrangler.multipole_to_local( traversal.level_start_target_or_target_parent_box_nrs, traversal.target_or_target_parent_boxes, - traversal.sep_siblings_starts, - traversal.sep_siblings_lists, + traversal.from_sep_siblings_starts, + traversal.from_sep_siblings_lists, mpole_exps) # }}} @@ -427,11 +427,11 @@ def drive_fmm(expansion_wrangler, src_weights): non_qbx_potentials = non_qbx_potentials + wrangler.eval_multipoles( traversal.level_start_target_box_nrs, traversal.target_boxes, - traversal.sep_smaller_by_level, + traversal.from_sep_smaller_by_level, mpole_exps) # assert that list 3 close has been merged into list 1 - assert traversal.sep_close_smaller_starts is None + assert traversal.from_sep_close_smaller_starts is None # }}} @@ -442,12 +442,12 @@ def drive_fmm(expansion_wrangler, src_weights): local_exps = local_exps + wrangler.form_locals( traversal.level_start_target_or_target_parent_box_nrs, traversal.target_or_target_parent_boxes, - traversal.sep_bigger_starts, - traversal.sep_bigger_lists, + traversal.from_sep_bigger_starts, + traversal.from_sep_bigger_lists, src_weights) # assert that list 4 close has been merged into list 1 - assert traversal.sep_close_bigger_starts is None + assert traversal.from_sep_close_bigger_starts is None # }}} @@ -624,7 +624,7 @@ def assemble_performance_data(geo_data, uses_pde_expansions, def process_list2(): nm2l = 0 for itgt_box, tgt_ibox in enumerate(traversal.target_or_target_parent_boxes): - start, end = traversal.sep_siblings_starts[itgt_box:itgt_box+2] + start, end = traversal.from_sep_siblings_starts[itgt_box:itgt_box+2] nm2l += (end-start) @@ -641,7 +641,7 @@ def assemble_performance_data(geo_data, uses_pde_expansions, for itgt_box, tgt_ibox in enumerate(traversal.target_boxes): ntargets = box_target_counts_nonchild[tgt_ibox] - for sep_smaller_list in traversal.sep_smaller_by_level: + for sep_smaller_list in traversal.from_sep_smaller_by_level: start, end = sep_smaller_list.starts[itgt_box:itgt_box+2] nmp_eval += ntargets * (end-start) @@ -658,9 +658,9 @@ def assemble_performance_data(geo_data, uses_pde_expansions, nform_local = 0 for itgt_box, tgt_ibox in enumerate(traversal.target_or_target_parent_boxes): - start, end = traversal.sep_bigger_starts[itgt_box:itgt_box+2] + start, end = traversal.from_sep_bigger_starts[itgt_box:itgt_box+2] - for src_ibox in traversal.sep_bigger_lists[start:end]: + for src_ibox in traversal.from_sep_bigger_lists[start:end]: nsources = tree.box_source_counts_nonchild[src_ibox] nform_local += nsources @@ -724,7 +724,7 @@ def assemble_performance_data(geo_data, uses_pde_expansions, def process_m2qbxl(): nqbx_m2l = 0 - for isrc_level, ssn in enumerate(traversal.sep_smaller_by_level): + for isrc_level, ssn in enumerate(traversal.from_sep_smaller_by_level): for itgt_center, tgt_icenter in enumerate(global_qbx_centers): icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter] diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index c47bca38817369d6ed1e1e7f6de94d5b257387b9..25ecbed94c54ea342a9bbf9e287a3ec4a82961a2 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -431,7 +431,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): mploc = self.get_translation_routine("%ddmploc", vec_suffix="_imany") for isrc_level, ssn in enumerate( - geo_data.traversal().sep_smaller_by_level): + geo_data.traversal().from_sep_smaller_by_level): source_level_start_ibox, source_mpoles_view = \ self.multipole_expansions_view(multipole_exps, isrc_level) diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 9a2914ae28d5db325232166ffe3a87f7ea8d61aa..4138fbb6902afbfb52ac6b74281861ebb81e1424 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -104,10 +104,11 @@ class target_state(Enum): # noqa class QBXFMMGeometryCodeGetter(object): - def __init__(self, cl_context, ambient_dim, debug): + def __init__(self, cl_context, ambient_dim, debug, _well_sep_is_n_away): self.cl_context = cl_context self.ambient_dim = ambient_dim self.debug = debug + self._well_sep_is_n_away = _well_sep_is_n_away @memoize_method def copy_targets_kernel(self): @@ -138,7 +139,9 @@ class QBXFMMGeometryCodeGetter(object): @memoize_method def build_traversal(self): from boxtree.traversal import FMMTraversalBuilder - return FMMTraversalBuilder(self.cl_context) + return FMMTraversalBuilder( + self.cl_context, + well_sep_is_n_away=self._well_sep_is_n_away) @memoize_method def qbx_center_to_target_box_lookup(self, particle_id_dtype, box_id_dtype): @@ -493,22 +496,11 @@ class QBXFMMGeometryData(object): refine_weights[:nsources] = 1 refine_weights.finish() - # NOTE: max_leaf_refine_weight has an impact on accuracy. - # For instance, if a panel contains 64*4 = 256 nodes, then - # a box with 128 sources will contain at most half a - # panel, meaning that its width will be on the order h/2, - # which means many QBX disks (diameter h) will be forced - # to cross boxes. - # So we set max_leaf_refine weight comfortably large - # to avoid having too many disks overlap more than one box. - # - # FIXME: Should investigate this further. - tree, _ = code_getter.build_tree(queue, particles=lpot_src.quad_stage2_density_discr.nodes(), targets=target_info.targets, target_radii=target_radii, - max_leaf_refine_weight=32, + max_leaf_refine_weight=lpot_src._max_leaf_refine_weight, refine_weights=refine_weights, debug=self.debug, stick_out_factor=lpot_src._expansion_stick_out_factor, diff --git a/requirements.txt b/requirements.txt index cdd2c65adaa5212ee49c73a6caa8c5884bed2125..7360a7c78e2a50cc32a45024544cdfb468da0934 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,10 +2,10 @@ numpy git+git://github.com/inducer/pymbolic sympy==1.0 git+https://github.com/inducer/modepy -git+https://github.com/pyopencl/pyopencl +git+https://github.com/inducer/pyopencl git+https://github.com/inducer/islpy git+https://github.com/inducer/loopy -git+https://github.com/inducer/boxtree +git+https://gitlab.tiker.net/inducer/boxtree git+https://github.com/inducer/meshmode git+https://github.com/inducer/sumpy git+https://github.com/inducer/pyfmmlib diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index 60bb570bdd59317bd261d0b545265f76871b0446..1e1c92135e903ff36db851e2ace0a1862116be98 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -34,6 +34,7 @@ from pyopencl.tools import ( # noqa from functools import partial from meshmode.mesh.generation import ( # noqa ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle, + NArmedStarfish, make_curve_mesh) # from sumpy.visualization import FieldPlotter from pytential import bind, sym, norm @@ -54,13 +55,6 @@ d1 = sym.Derivative() d2 = sym.Derivative() -def get_wobbly_circle_mesh(refinement_increment, target_order): - nelements = [3000, 5000, 7000][refinement_increment] - return make_curve_mesh(WobblyCircle.random(30, seed=30), - np.linspace(0, 1, nelements+1), - target_order) - - def get_sphere_mesh(refinement_increment, target_order): from meshmode.mesh.generation import generate_icosphere mesh = generate_icosphere(1, target_order) @@ -76,15 +70,38 @@ def get_sphere_mesh(refinement_increment, target_order): class StarfishGeometry(object): - mesh_name = "starfish" + def __init__(self, n_arms=5, amplitude=0.25): + self.n_arms = n_arms + self.amplitude = amplitude + + @property + def mesh_name(self): + return "%d-starfish-%s" % ( + self.n_arms, + self.amplitude) + dim = 2 resolutions = [30, 50, 70] def get_mesh(self, nelements, target_order): - return make_curve_mesh(starfish, - np.linspace(0, 1, nelements+1), - target_order) + return make_curve_mesh( + NArmedStarfish(self.n_arms, self.amplitude), + np.linspace(0, 1, nelements+1), + target_order) + + +class WobblyCircleGeometry(object): + dim = 2 + mesh_name = "wobbly-circle" + + resolutions = [2000, 3000, 4000] + + def get_mesh(self, resolution, target_order): + return make_curve_mesh( + WobblyCircle.random(30, seed=30), + np.linspace(0, 1, resolution+1), + target_order) class SphereGeometry(object): @@ -155,6 +172,32 @@ class StaticTestCase(object): pass +class StarfishGreenTest(StaticTestCase): + expr = GreenExpr() + geometry = StarfishGeometry() + k = 0 + qbx_order = 3 + fmm_order = 6 + + resolutions = [30, 50, 70] + + _expansion_stick_out_factor = 0.5 + + fmm_backend = "sumpy" + + +class WobblyCircleGreenTest(StaticTestCase): + expr = GreenExpr() + geometry = WobblyCircleGeometry() + k = 0 + qbx_order = 3 + fmm_order = 10 + + _expansion_stick_out_factor = 0.5 + + fmm_backend = "sumpy" + + class SphereGreenTest(StaticTestCase): expr = GreenExpr() geometry = SphereGeometry() @@ -164,7 +207,7 @@ class SphereGreenTest(StaticTestCase): resolutions = [0, 1] - _expansion_stick_out_factor = 0.75 + _expansion_stick_out_factor = 0.5 fmm_backend = "fmmlib" @@ -334,7 +377,7 @@ def test_identity_convergence(ctx_getter, case, visualize=False): from meshmode.discretization.visualization import make_visualizer bdry_vis = make_visualizer(queue, density_discr, target_order) - bdry_normals = bind(density_discr, sym.normal(3))(queue)\ + bdry_normals = bind(density_discr, sym.normal(mesh.ambient_dim))(queue)\ .as_vector(dtype=object) bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [