From 679880360add5dd0327d154ec02af68e1d16954b Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 12 May 2017 14:54:41 -0500 Subject: [PATCH 01/22] [ci skip] Start refactoring the refiner towards 3D refinement (#35). * Also start using the container and wrangler pattern (#34). --- pytential/qbx/refinement.py | 587 ++++++++++++------------------------ test/test_global_qbx.py | 22 +- 2 files changed, 197 insertions(+), 412 deletions(-) diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index c942b45b..f2984ea5 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -3,7 +3,7 @@ from __future__ import division, absolute_import, print_function __copyright__ = """ Copyright (C) 2013 Andreas Kloeckner -Copyright (C) 2016 Matt Wala +Copyright (C) 2016, 2017 Matt Wala """ __license__ = """ @@ -32,7 +32,6 @@ import numpy as np import pyopencl as cl from pytools import memoize_method -from pytential.qbx.utils import DiscrPlotterMixin from boxtree.area_query import AreaQueryElementwiseTemplate from pyopencl.elementwise import ElementwiseTemplate from boxtree.tools import InlineBinarySearch @@ -43,11 +42,32 @@ unwrap_args = AreaQueryElementwiseTemplate.unwrap_args import logging logger = logging.getLogger(__name__) + +MAX_LEVELS_INCREMENT = 10 + + __doc__ = """ Refinement ^^^^^^^^^^ -.. autoclass:: QBXLayerPotentialSourceRefiner +The refiner takes a layer potential source and ensures that it satisfies three +QBX refinement criteria: + + * *Condition 1* (Expansion disk undisturbed by sources) + A center must be closest to its own source. + + * *Condition 2* (Sufficient quadrature sampling from all source panels) + The quadrature contribution from each panel is as accurate + as from the center's own source panel. + + * *Condition 3* (Panel size bounded based on wavelength) + (Helmholtz only) The panel size is bounded with respect to the wavelength. + +.. autoclass:: RefinerCodeContainer +.. autoclass:: RefinerWrangler + +.. automethod:: make_empty_refine_flags +.. automethod:: refine_for_global_qbx """ # {{{ kernels @@ -110,10 +130,8 @@ TUNNEL_QUERY_DISTANCE_FINDER_TEMPLATE = ElementwiseTemplate( preamble=str(InlineBinarySearch("particle_id_t"))) -# Implements "Algorithm for triggering refinement based on Condition 1" -# -# Does not use r_max as we do not use Newton for checking center-panel closeness. -CENTER_IS_CLOSEST_TO_ORIG_PANEL_REFINER = AreaQueryElementwiseTemplate( +# Refinement checker for Condition 1. +CENTER_IS_CLOSEST_TO_ORIG_PANEL_CHECKER = AreaQueryElementwiseTemplate( extra_args=r""" /* input */ particle_id_t *box_to_panel_starts, @@ -182,119 +200,16 @@ CENTER_IS_CLOSEST_TO_ORIG_PANEL_REFINER = AreaQueryElementwiseTemplate( preamble=str(InlineBinarySearch("particle_id_t"))) -# Implements "Algorithm for triggering refinement based on Condition 2" -CENTER_IS_FAR_FROM_NONNEIGHBOR_PANEL_REFINER = AreaQueryElementwiseTemplate( - extra_args=r""" - /* input */ - particle_id_t *box_to_panel_starts, - particle_id_t *box_to_panel_lists, - particle_id_t *panel_to_source_starts, - particle_id_t *panel_to_center_starts, - particle_id_t source_offset, - particle_id_t center_offset, - particle_id_t panel_offset, - particle_id_t *sorted_target_ids, - coord_t *panel_sizes, - particle_id_t *panel_adjacency_starts, - particle_id_t *panel_adjacency_lists, - int npanels, - coord_t *tunnel_query_dists, - - /* output */ - int *panel_refine_flags, - int *found_panel_to_refine, - - /* input, dim-dependent length */ - %for ax in AXIS_NAMES[:dimensions]: - coord_t *particles_${ax}, - %endfor - """, - ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r""" - particle_id_t my_panel = bsearch(panel_to_center_starts, npanels + 1, i); - coord_vec_t my_center_coords; - - ${load_particle("INDEX_FOR_CENTER_PARTICLE(i)", "my_center_coords")} - ${load_particle("INDEX_FOR_PANEL_PARTICLE(my_panel)", ball_center)} - ${ball_radius} = tunnel_query_dists[my_panel]; - """, - leaf_found_op=QBX_TREE_MAKO_DEFS + r""" - for (particle_id_t panel_idx = box_to_panel_starts[${leaf_box_id}]; - panel_idx < box_to_panel_starts[${leaf_box_id} + 1]; - ++panel_idx) - { - particle_id_t panel = box_to_panel_lists[panel_idx]; - - bool is_self_or_adjacent = (my_panel == panel); - - for (particle_id_t adj_panel_idx = panel_adjacency_starts[my_panel]; - adj_panel_idx < panel_adjacency_starts[my_panel + 1]; - ++adj_panel_idx) - { - is_self_or_adjacent |= ( - panel_adjacency_lists[adj_panel_idx] == panel); - } - - // Skip self and adjacent panels. - if (is_self_or_adjacent) - { - continue; - } - - bool is_close = false; - - for (particle_id_t source = panel_to_source_starts[panel]; - source < panel_to_source_starts[panel + 1]; - ++source) - { - coord_vec_t source_coords; - - ${load_particle( - "INDEX_FOR_SOURCE_PARTICLE(source)", "source_coords")} - - is_close |= ( - distance(my_center_coords, source_coords) - <= panel_sizes[panel] / 2); - } - - if (is_close) - { - panel_refine_flags[my_panel] = 1; - *found_panel_to_refine = 1; - break; - } - } - """, - name="refine_center_far_from_nonneighbor_panels", - preamble=str(InlineBinarySearch("particle_id_t"))) - -# }}} - - -# {{{ lpot source refiner - -class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin): - """ - Driver for refining the QBX source grid. Follows [1]_. - - .. [1] Rachh, Manas, Andreas Klöckner, and Michael O'Neil. "Fast - algorithms for Quadrature by Expansion I: Globally valid expansions." - - .. automethod:: get_refine_flags - .. automethod:: __call__ - """ +# {{{ refiner code container - def __init__(self, context): - self.cl_context = context - from pytential.qbx.utils import TreeWithQBXMetadataBuilder - self.tree_builder = TreeWithQBXMetadataBuilder(self.cl_context) - from boxtree.area_query import PeerListFinder - self.peer_list_finder = PeerListFinder(self.cl_context) +class RefinerCodeContainer(object): - # {{{ kernels + def __init__(self, cl_context): + self.cl_context = cl_context @memoize_method - def get_tunnel_query_distance_finder(self, dimensions, coord_dtype, - particle_id_dtype): + def tunnel_query_distance_finder( + self, dimensions, coord_dtype, particle_id_dtype): from pyopencl.tools import dtype_to_ctype from boxtree.tools import AXIS_NAMES logger.info("refiner: building tunnel query distance finder kernel") @@ -317,65 +232,16 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin): return knl @memoize_method - def get_center_is_closest_to_orig_panel_refiner(self, dimensions, - coord_dtype, box_id_dtype, - peer_list_idx_dtype, - particle_id_dtype, - max_levels): - return CENTER_IS_CLOSEST_TO_ORIG_PANEL_REFINER.generate(self.cl_context, - dimensions, coord_dtype, box_id_dtype, peer_list_idx_dtype, - max_levels, - extra_type_aliases=(("particle_id_t", particle_id_dtype),)) - - @memoize_method - def get_center_is_far_from_nonneighbor_panel_refiner(self, dimensions, - coord_dtype, - box_id_dtype, - peer_list_idx_dtype, - particle_id_dtype, - max_levels): - return CENTER_IS_FAR_FROM_NONNEIGHBOR_PANEL_REFINER.generate(self.cl_context, + def center_is_closest_to_orig_panel_checker( + self, dimensions, coord_dtype, box_id_dtype, peer_list_idx_dtype, + particle_id_dtype, max_levels): + return CENTER_IS_CLOSEST_TO_ORIG_PANEL_CHECKER.generate(self.cl_context, dimensions, coord_dtype, box_id_dtype, peer_list_idx_dtype, max_levels, extra_type_aliases=(("particle_id_t", particle_id_dtype),)) @memoize_method - def get_2_to_1_panel_ratio_refiner(self): - knl = lp.make_kernel([ - "{[panel]: 0<=panel neighbor_start = panel_adjacency_starts[panel] - <> neighbor_stop = panel_adjacency_starts[panel + 1] - for ineighbor - <> neighbor = panel_adjacency_lists[ineighbor] - - <> oversize = (refine_flags_prev[panel] == 0 - and ( - (panel_sizes[panel] > 2 * panel_sizes[neighbor]) or - (panel_sizes[panel] > panel_sizes[neighbor] and - refine_flags_prev[neighbor] == 1))) - - if oversize - refine_flags[panel] = 1 - refine_flags_updated = 1 {id=write_refine_flags_updated} - end - end - end - """, [ - lp.GlobalArg("panel_adjacency_lists", shape=None), - "..." - ], - options="return_dict", - silenced_warnings="write_race(write_refine_flags_updated)", - name="refine_2_to_1_adj_panel_size_ratio") - knl = lp.split_iname(knl, "panel", 128, inner_tag="l.0", outer_tag="g.0") - return knl - - @memoize_method - def get_helmholtz_k_to_panel_ratio_refiner(self): + def helmholtz_k_to_panel_size_ratio_checker(self): knl = lp.make_kernel( "{[panel]: 0<=panel npanels_to_refine_prev: - logger.debug("refiner: found {} panel(s) to refine".format( - npanels_to_refine - npanels_to_refine_prev)) +class RefinerWrangler(object): - logger.info("refiner: done checking center is closest to orig panel") + def __init__(self, code_container, queue, refiner): + self.code_container = code_container + self.queue = queue + self.refiner = refiner - return found_panel_to_refine.get()[0] == 1 + def check_center_is_closest_to_orig_panel(self, + lpot_source, tree, peer_lists, refine_flags, + debug, wait_for=None): - def refinement_check_center_is_far_from_nonneighbor_panels(self, queue, - tree, lpot_source, peer_lists, tq_dists, refine_flags, debug, - wait_for=None): # Avoid generating too many kernels. from pytools import div_ceil - max_levels = 10 * div_ceil(tree.nlevels, 10) + max_levels = MAX_LEVELS_INCREMENT * div_ceil( + tree.nlevels, MAX_LEVELS_INCREMENT) - knl = self.get_center_is_far_from_nonneighbor_panel_refiner( + knl = self.code_container.center_is_closest_to_orig_panel_checker( tree.dimensions, tree.coord_dtype, tree.box_id_dtype, peer_lists.peer_list_starts.dtype, tree.particle_id_dtype, max_levels) - logger.info("refiner: checking center is far from nonneighbor panels") + logger.info("refiner: checking center is closest to orig panel") if debug: npanels_to_refine_prev = cl.array.sum(refine_flags).get() - found_panel_to_refine = cl.array.zeros(queue, 1, np.int32) + found_panel_to_refine = cl.array.zeros(self.queue, 1, np.int32) found_panel_to_refine.finish() - adjacency = self.get_adjacency_on_device(queue, lpot_source) - evt = knl( *unwrap_args( tree, peer_lists, @@ -481,18 +321,14 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin): tree.qbx_panel_to_center_starts, tree.qbx_user_source_slice.start, tree.qbx_user_center_slice.start, - tree.qbx_user_panel_slice.start, tree.sorted_target_ids, lpot_source.panel_sizes("npanels"), - adjacency.adjacency_starts, - adjacency.adjacency_lists, tree.nqbxpanels, - tq_dists, refine_flags, found_panel_to_refine, *tree.sources), range=slice(tree.nqbxcenters), - queue=queue) + queue=self.queue) cl.wait_for_events([evt]) @@ -502,20 +338,23 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin): logger.debug("refiner: found {} panel(s) to refine".format( npanels_to_refine - npanels_to_refine_prev)) - logger.info("refiner: done checking center is far from nonneighbor panels") + logger.info("refiner: done checking center is closest to orig panel") return found_panel_to_refine.get()[0] == 1 - def refinement_check_helmholtz_k_to_panel_size_ratio(self, queue, lpot_source, + def check_sufficient_quadrature_resolution(self): + return True + + def check_helmholtz_k_to_panel_size_ratio(self, lpot_source, helmholtz_k, refine_flags, debug, wait_for=None): - knl = self.get_helmholtz_k_to_panel_ratio_refiner() + knl = self.code_container.helmholtz_k_to_panel_size_ratio_checker() logger.info("refiner: checking helmholtz k to panel size ratio") if debug: npanels_to_refine_prev = cl.array.sum(refine_flags).get() - evt, out = knl(queue, + evt, out = knl(self.queue, panel_sizes=lpot_source.panel_sizes("npanels"), refine_flags=refine_flags, refine_flags_updated=np.array(0), @@ -534,53 +373,9 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin): return (out["refine_flags_updated"].get() == 1).all() - def refinement_check_2_to_1_panel_size_ratio(self, queue, lpot_source, - refine_flags, debug, wait_for=None): - knl = self.get_2_to_1_panel_ratio_refiner() - adjacency = self.get_adjacency_on_device(queue, lpot_source) - - refine_flags_updated = False - - logger.info("refiner: checking 2-to-1 panel size ratio") - - if debug: - npanels_to_refine_prev = cl.array.sum(refine_flags).get() - - # Iterative refinement until no more panels can be marked - while True: - evt, out = knl(queue, - npanels=lpot_source.density_discr.mesh.nelements, - panel_sizes=lpot_source.panel_sizes("npanels"), - refine_flags=refine_flags, - # It's safe to pass this here, as the resulting data - # race won't affect the final result of the - # computation. - refine_flags_prev=refine_flags, - refine_flags_updated=np.array(0), - panel_adjacency_starts=adjacency.adjacency_starts, - panel_adjacency_lists=adjacency.adjacency_lists, - wait_for=wait_for) - - cl.wait_for_events([evt]) - - if (out["refine_flags_updated"].get() == 1).all(): - refine_flags_updated = True - else: - break - - if debug: - npanels_to_refine = cl.array.sum(refine_flags).get() - if npanels_to_refine > npanels_to_refine_prev: - logger.debug("refiner: found {} panel(s) to refine".format( - npanels_to_refine - npanels_to_refine_prev)) - - logger.info("refiner: done checking 2-to-1 panel size ratio") - - return refine_flags_updated - # }}} - # {{{ other utilities + # {{{ def get_tunnel_query_dists(self, queue, tree, lpot_source): """ @@ -612,164 +407,166 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin): return tq_dists, evt - def get_adjacency_on_device(self, queue, lpot_source): - """ - Take adjacency information from the mesh and place it onto the device. - """ - from boxtree.tools import DeviceDataRecord - adjacency = lpot_source.density_discr.mesh.nodal_adjacency - adjacency_starts = cl.array.to_device(queue, adjacency.neighbors_starts) - adjacency_lists = cl.array.to_device(queue, adjacency.neighbors) - - return DeviceDataRecord( - adjacency_starts=adjacency_starts, - adjacency_lists=adjacency_lists) - - def get_refine_flags(self, queue, lpot_source): - """ - Return an array on the device suitable for use as element refine flags. - - :arg queue: An instance of :class:`pyopencl.CommandQueue`. - :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`. + # }}} - :returns: An instance of :class:`pyopencl.array.Array` suitable for - use as refine flags, initialized to zero. - """ - result = cl.array.zeros( - queue, lpot_source.density_discr.mesh.nelements, np.int32) - return result, result.events[0] + def build_tree(self, lpot_source): + tb = self.code_container.tree_builder() + return tb(self.queue, lpot_source) - # }}} + def find_peer_lists(self, tree): + plf = self.code_container.peer_list_finder() + peer_lists, evt = plf(self.queue, tree) + cl.wait_for_events([evt]) + return peer_lists - def refine(self, queue, lpot_source, refine_flags, refiner, factory, debug): + def refine(self, lpot_source, refine_flags, factory, debug): """ Refine the underlying mesh and discretization. """ if isinstance(refine_flags, cl.array.Array): - refine_flags = refine_flags.get(queue) + refine_flags = refine_flags.get(self.queue) refine_flags = refine_flags.astype(np.bool) logger.info("refiner: calling meshmode") - refiner.refine(refine_flags) + self.refiner.refine(refine_flags) from meshmode.discretization.connection import make_refinement_connection conn = make_refinement_connection( - refiner, lpot_source.density_discr, + self.refiner, lpot_source.density_discr, factory) logger.info("refiner: done calling meshmode") - new_density_discr = conn.to_discr - new_lpot_source = lpot_source.copy( - density_discr=new_density_discr, - refined_for_global_qbx=True) + density_discr=conn.to_discr) return new_lpot_source, conn - def __call__(self, lpot_source, discr_factory, helmholtz_k=None, - # FIXME: Set debug=False once everything works. - refine_flags=None, debug=True, maxiter=50): - """ - Entry point for calling the refiner. +# }}} + - :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`. +def make_empty_refine_flags(queue, lpot_source): + """Return an array on the device suitable for use as element refine flags. - :arg group_factory: An instance of - :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for - discretizing the refined mesh. + :arg queue: An instance of :class:`pyopencl.CommandQueue`. + :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`. - :arg helmholtz_k: The Helmholtz parameter, or `None` if not applicable. + :returns: A :class:`pyopencl.array.Array` suitable for use as refine flags, + initialized to zero. + """ + result = cl.array.zeros( + queue, lpot_source.density_discr.mesh.nelements, np.int32) + result.finish() + return result - :arg refine_flags: A :class:`pyopencl.array.Array` indicating which - panels should get refined initially, or `None` if no initial - refinement should be done. Should have size equal to the number of - panels. See also :meth:`get_refine_flags()`. - :returns: A tuple ``(lpot_source, conns)`` where ``lpot_source`` is the - refined layer potential source, and ``conns`` is a list of - :class:`meshmode.discretization.connection.DiscretizationConnection` - objects going from the original mesh to the refined mesh. - """ - from meshmode.mesh.refinement import Refiner - refiner = Refiner(lpot_source.density_discr.mesh) - connections = [] - - lpot_source = lpot_source.copy( - debug=debug, - refined_for_global_qbx=True) - - with cl.CommandQueue(self.cl_context) as queue: - if refine_flags: - lpot_source, conn = self.refine( - queue, lpot_source, refine_flags, refiner, discr_factory, - debug) - connections.append(conn) +def refine_for_global_qbx(lpot_source, code_container, + group_factory, helmholtz_k=None, + # FIXME: Set debug=False once everything works. + refine_flags=None, debug=True, maxiter=50): + """ + Entry point for calling the refiner. + + :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`. - done_refining = False - niter = 0 + :arg code_container: An instance of :class:`RefinerCodeContainer`. - while not done_refining: - niter += 1 + :arg group_factory: An instance of + :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for + discretizing the refined mesh. - if niter > maxiter: - logger.warning( - "Max iteration count reached in QBX layer potential source" - " refiner.") - break + :arg helmholtz_k: The Helmholtz parameter, or `None` if not applicable. - # Build tree and auxiliary data. - # FIXME: The tree should not have to be rebuilt at each iteration. - tree = self.tree_builder(queue, lpot_source) - wait_for = [] + :arg refine_flags: A :class:`pyopencl.array.Array` indicating which + panels should get refined initially, or `None` if no initial + refinement should be done. Should have size equal to the number of + panels. See also :func:`make_empty_refine_flags()`. - peer_lists, evt = self.peer_list_finder(queue, tree, wait_for) - wait_for = [evt] + :arg maxiter: The maximum number of refiner iterations. - refine_flags, evt = self.get_refine_flags(queue, lpot_source) - wait_for.append(evt) + :returns: A tuple ``(lpot_source, conns)`` where ``lpot_source`` is the + refined layer potential source, and ``conns`` is a list of + :class:`meshmode.discretization.connection.DiscretizationConnection` + objects going from the original mesh to the refined mesh. + """ + + # Algorithm: + # + # 1. Do initial refinement, if requested. + # + # 2. While not converged: + # Refine until each center is closest to its own + # source and panel sizes bounded by wavelength. + # + # 3. While not converged: + # Refine fine density discretization until + # sufficient quadrature resolution from all panels achieved. + # + # TODO: Stop doing redundant checks by avoiding panels which no longer need + # refinement. + + from meshmode.mesh.refinement import Refiner + refiner = Refiner(lpot_source.density_discr.mesh) + connections = [] + + with cl.CommandQueue(lpot_source.cl_context) as queue: + wrangler = code_container.get_wrangler(queue, refiner) + + # Do initial refinement. + if refine_flags is not None: + lpot_source, conn = wrangler.refine(lpot_source, refine_flags, + group_factory, debug) + connections.append(conn) + + # {{{ first stage refinement + + must_refine = True + niter = 0 + + while must_refine: + must_refine = False + niter += 1 + + if niter > maxiter: + logger.warning( + "Max iteration count reached in QBX layer potential source" + " refiner (first stage).") + break - tq_dists, evt = self.get_tunnel_query_dists(queue, tree, lpot_source) - wait_for.append(evt) + # Build tree and auxiliary data. + # FIXME: The tree should not have to be rebuilt at each iteration. + tree = wrangler.build_tree(lpot_source) + peer_lists = wrangler.find_peer_lists(tree) + refine_flags = make_empty_refine_flags(queue, lpot_source) - # Run refinement checkers. - must_refine = False + # Check condition 1. + must_refine |= wrangler.check_center_is_closest_to_orig_panel( + lpot_source, tree, peer_lists, refine_flags, debug) - must_refine |= \ - self.refinement_check_center_is_closest_to_orig_panel( - queue, tree, lpot_source, peer_lists, refine_flags, - debug, wait_for) + # Check condition 3. + if helmholtz_k is not None: + must_refine |= wrangler.check_helmholtz_k_to_panel_size_ratio( + lpot_source, helmholtz_k, refine_flags, debug) - must_refine |= \ - self.refinement_check_center_is_far_from_nonneighbor_panels( - queue, tree, lpot_source, peer_lists, tq_dists, - refine_flags, debug, wait_for) + if must_refine: + lpot_source, conn = wrangler.refine(lpot_source, refine_flags, + group_factory, debug) + connections.append(conn) - must_refine |= \ - self.refinement_check_2_to_1_panel_size_ratio( - queue, lpot_source, refine_flags, debug, wait_for) + del tree + del refine_flags + del peer_lists - if helmholtz_k: - must_refine |= \ - self.refinement_check_helmholtz_k_to_panel_size_ratio( - queue, lpot_source, helmholtz_k, refine_flags, debug, - wait_for) + # }}} - if must_refine: - lpot_source, conn = self.refine( - queue, lpot_source, refine_flags, refiner, discr_factory, - debug) - connections.append(conn) + # {{{ second stage refinement - del tree - del peer_lists - del tq_dists - del refine_flags - done_refining = not must_refine + # }}} - return lpot_source, connections + lpot_source = lpot_source.copy(debug=debug, refined_for_global_qbx=True) + + return lpot_source, connections -# }}} # vim: foldmethod=marker:filetype=pyopencl diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index d7e51451..34e878ab 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -112,13 +112,15 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): discr = Discretization(cl_ctx, mesh, factory) - from pytential.qbx.refinement import QBXLayerPotentialSourceRefiner + from pytential.qbx.refinement import ( + RefinerCodeContainer, refine_for_global_qbx) lpot_source = QBXLayerPotentialSource(discr, order) del discr - refiner = QBXLayerPotentialSourceRefiner(cl_ctx) - lpot_source, conn = refiner(lpot_source, factory, helmholtz_k) + lpot_source, conn = refine_for_global_qbx( + lpot_source, RefinerCodeContainer(cl_ctx), + factory, helmholtz_k) discr_nodes = lpot_source.density_discr.nodes().get(queue) int_centers = lpot_source.centers(-1) @@ -132,11 +134,6 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): # {{{ check if satisfying criteria def check_panel(panel): - # Check 2-to-1 panel to neighbor size ratio. - for neighbor in panel.neighbors: - assert panel_sizes[panel.element_nr] / panel_sizes[neighbor] <= 2, \ - (panel_sizes[panel.element_nr], panel_sizes[neighbor]) - # Check wavenumber to panel size ratio. assert panel_sizes[panel.element_nr] * helmholtz_k <= 5 @@ -165,15 +162,6 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): assert dist >= h_1 / 2, (dist, h_1, panel_1.element_nr, panel_2.element_nr) - # Criterion 2: - # A center cannot be closer to another panel than that panel's - # centers - unless the panels are adjacent, to allow for refinement. - - if panel_2.element_nr in panel_1.neighbors: - return - - assert dist >= h_2 / 2, (dist, h_2, panel_1.element_nr, panel_2.element_nr) - for panel_1 in iter_elements(lpot_source.density_discr): check_panel(panel_1) for panel_2 in iter_elements(lpot_source.density_discr): -- GitLab From c7772c3624ddbaef53055b98ae58f61e95567e14 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 14 May 2017 19:35:10 -0500 Subject: [PATCH 02/22] [ci skip] Finish implementing new refinement criteria. --- doc/qbx.rst | 5 + pytential/qbx/__init__.py | 88 ++++++-- pytential/qbx/refinement.py | 401 ++++++++++++++++++++-------------- pytential/qbx/target_assoc.py | 16 +- pytential/qbx/utils.py | 373 ++++++++++++++++--------------- test/test_global_qbx.py | 67 ++++-- 6 files changed, 565 insertions(+), 385 deletions(-) diff --git a/doc/qbx.rst b/doc/qbx.rst index d297db5e..46244da1 100644 --- a/doc/qbx.rst +++ b/doc/qbx.rst @@ -3,6 +3,11 @@ QBX internals ============= +Refinement +---------- + +.. automodule:: pytential.qbx.refinement + Data structures describing geometry ----------------------------------- diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index a18196f4..ea96e87c 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -29,7 +29,6 @@ import loopy as lp import numpy as np from pytools import memoize_method from meshmode.discretization import Discretization -from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory from pytential.qbx.target_assoc import QBXTargetAssociationFailedException from pytential.source import PotentialSource @@ -141,6 +140,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_order=None, fmm_level_to_order=None, target_stick_out_factor=1e-10, + fine_density_discr=None, + resampler=None, # begin undocumented arguments # FIXME default debug=False once everything works @@ -150,10 +151,21 @@ class QBXLayerPotentialSource(LayerPotentialSource): """ :arg fine_order: The total degree to which the (upsampled) underlying quadrature is exact. + :arg fine_density_discr: A discretization or *None.* + If non-*None*, should also supply *resampler.* + :arg resampler: A connection used for resampling. + If non-*None*, should also supply *fine_density_discr*. :arg fmm_order: `False` for direct calculation. ``None`` will set a reasonable(-ish?) default. """ + if fine_density_discr is None and resampler is not None: + raise ValueError( + "resampler is supplied; must also supply fine_density_discr") + if fine_density_discr is not None and resampler is None: + raise ValueError( + "fine_density_discr is supplied; must also supply resampler") + if fmm_level_to_order is None: if fmm_order is None and qbx_order is not None: fmm_order = qbx_order + 1 @@ -174,6 +186,10 @@ class QBXLayerPotentialSource(LayerPotentialSource): self.fmm_level_to_order = fmm_level_to_order self.target_stick_out_factor = target_stick_out_factor + # Default values are lazily provided if these are None + self._fine_density_discr = fine_density_discr + self._resampler = resampler + self.debug = debug self.refined_for_global_qbx = refined_for_global_qbx self.performance_data_file = performance_data_file @@ -185,6 +201,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): qbx_order=None, fmm_level_to_order=None, target_stick_out_factor=None, + fine_density_discr=None, + resampler=None, debug=None, refined_for_global_qbx=None, @@ -200,6 +218,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_level_to_order or self.fmm_level_to_order), target_stick_out_factor=( target_stick_out_factor or self.target_stick_out_factor), + fine_density_discr=fine_density_discr or self._fine_density_discr, + resampler=resampler or self._resampler, debug=( debug if debug is not None else self.debug), @@ -211,13 +231,23 @@ class QBXLayerPotentialSource(LayerPotentialSource): @property @memoize_method def fine_density_discr(self): + if self._fine_density_discr is not None: + return self._fine_density_discr + + from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory) + return Discretization( self.density_discr.cl_context, self.density_discr.mesh, - QuadratureSimplexGroupFactory(self.fine_order), self.real_dtype) + InterpolatoryQuadratureSimplexGroupFactory(self.fine_order), + self.real_dtype) @property @memoize_method def resampler(self): + if self._resampler is not None: + return self._resampler + from meshmode.discretization.connection import make_same_mesh_connection return make_same_mesh_connection( self.fine_density_discr, self.density_discr) @@ -238,23 +268,35 @@ class QBXLayerPotentialSource(LayerPotentialSource): global_array.shape[:-1] + (group.nelements,)) + @property + @memoize_method + def refiner_code_container(self): + from pytential.qbx.refinement import RefinerCodeContainer + return RefinerCodeContainer(self.cl_context) + @memoize_method - def with_refinement(self, target_order=None, maxiter=3): + def with_refinement(self, target_order=None, maxiter=10): """ :returns: a tuple ``(lpot_src, cnx)``, where ``lpot_src`` is a :class:`QBXLayerPotentialSource` and ``cnx`` is a :class:`meshmode.discretization.connection.DiscretizationConnection` from the originally given to the refined geometry. """ - from pytential.qbx.refinement import QBXLayerPotentialSourceRefiner - refiner = QBXLayerPotentialSourceRefiner(self.cl_context) + from pytential.qbx.refinement import refine_for_global_qbx + from meshmode.discretization.poly_element import ( - InterpolatoryQuadratureSimplexGroupFactory) + InterpolatoryQuadratureSimplexGroupFactory) + if target_order is None: target_order = self.density_discr.groups[0].order - lpot, connection = refiner(self, + + lpot, connection = refine_for_global_qbx( + self, + self.refiner_code_container, InterpolatoryQuadratureSimplexGroupFactory(target_order), + InterpolatoryQuadratureSimplexGroupFactory(self.fine_order), maxiter=maxiter) + return lpot, connection @property @@ -284,8 +326,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): def complex_dtype(self): return self.density_discr.complex_dtype - @memoize_method - def panel_centers_of_mass(self): + def _centers_of_mass_for_discr(self, discr): knl = lp.make_kernel( """{[dim,k,i]: 0<=dim npanels_to_refine_prev: + logger.debug("refiner: found {} panel(s) to refine".format( + npanels_to_refine - npanels_to_refine_prev)) + + logger.info("refiner: done checking adequate quadrature resolution") + + return found_panel_to_refine.get()[0] == 1 def check_helmholtz_k_to_panel_size_ratio(self, lpot_source, - helmholtz_k, refine_flags, debug, wait_for=None): + helmholtz_k, refine_flags, debug, wait_for=None): knl = self.code_container.helmholtz_k_to_panel_size_ratio_checker() logger.info("refiner: checking helmholtz k to panel size ratio") @@ -373,45 +418,11 @@ class RefinerWrangler(object): return (out["refine_flags_updated"].get() == 1).all() - # }}} - - # {{{ - - def get_tunnel_query_dists(self, queue, tree, lpot_source): - """ - Compute radii for the tubular neighborhood around each panel center of mass. - """ - nqbxpanels = lpot_source.density_discr.mesh.nelements - # atomic_max only works on float32 - tq_dists = cl.array.zeros(queue, nqbxpanels, np.float32) - tq_dists.finish() - - knl = self.get_tunnel_query_distance_finder(tree.dimensions, - tree.coord_dtype, tree.particle_id_dtype) - - evt = knl(tree.qbx_user_source_slice.start, - tree.qbx_user_panel_slice.start, - nqbxpanels, - tree.qbx_panel_to_source_starts, - tree.sorted_target_ids, - lpot_source.panel_sizes("npanels"), - tq_dists, - *tree.sources, - queue=queue, - range=slice(tree.nqbxsources)) - - cl.wait_for_events([evt]) - - if tree.coord_dtype != tq_dists.dtype: - tq_dists = tq_dists.astype(tree.coord_dtype) - - return tq_dists, evt - - # }}} - - def build_tree(self, lpot_source): + def build_tree(self, lpot_source, use_fine_discr=False): tb = self.code_container.tree_builder() - return tb(self.queue, lpot_source) + from pytential.qbx.utils import build_tree_with_qbx_metadata + return build_tree_with_qbx_metadata( + self.queue, tb, lpot_source, use_fine_discr=use_fine_discr) def find_peer_lists(self, tree): plf = self.code_container.peer_list_finder() @@ -419,7 +430,7 @@ class RefinerWrangler(object): cl.wait_for_events([evt]) return peer_lists - def refine(self, lpot_source, refine_flags, factory, debug): + def refine(self, density_discr, refiner, refine_flags, factory, debug): """ Refine the underlying mesh and discretization. """ @@ -429,24 +440,18 @@ class RefinerWrangler(object): logger.info("refiner: calling meshmode") - self.refiner.refine(refine_flags) + refiner.refine(refine_flags) from meshmode.discretization.connection import make_refinement_connection - - conn = make_refinement_connection( - self.refiner, lpot_source.density_discr, - factory) + conn = make_refinement_connection(refiner, density_discr, factory) logger.info("refiner: done calling meshmode") - new_lpot_source = lpot_source.copy( - density_discr=conn.to_discr) - - return new_lpot_source, conn + return conn # }}} -def make_empty_refine_flags(queue, lpot_source): +def make_empty_refine_flags(queue, lpot_source, use_fine_discr=False): """Return an array on the device suitable for use as element refine flags. :arg queue: An instance of :class:`pyopencl.CommandQueue`. @@ -455,14 +460,16 @@ def make_empty_refine_flags(queue, lpot_source): :returns: A :class:`pyopencl.array.Array` suitable for use as refine flags, initialized to zero. """ - result = cl.array.zeros( - queue, lpot_source.density_discr.mesh.nelements, np.int32) + discr = (lpot_source.fine_density_discr + if use_fine_discr + else lpot_source.density_discr) + result = cl.array.zeros(queue, discr.mesh.nelements, np.int32) result.finish() return result def refine_for_global_qbx(lpot_source, code_container, - group_factory, helmholtz_k=None, + group_factory, fine_group_factory, helmholtz_k=None, # FIXME: Set debug=False once everything works. refine_flags=None, debug=True, maxiter=50): """ @@ -474,7 +481,11 @@ def refine_for_global_qbx(lpot_source, code_container, :arg group_factory: An instance of :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for - discretizing the refined mesh. + discretizing the coarse refined mesh. + + :arg fine_group_factory: An instance of + :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for + discretizing the fine refined mesh. :arg helmholtz_k: The Helmholtz parameter, or `None` if not applicable. @@ -485,10 +496,10 @@ def refine_for_global_qbx(lpot_source, code_container, :arg maxiter: The maximum number of refiner iterations. - :returns: A tuple ``(lpot_source, conns)`` where ``lpot_source`` is the - refined layer potential source, and ``conns`` is a list of + :returns: A tuple ``(lpot_source, *conn*)`` where ``lpot_source`` is the + refined layer potential source, and ``conn`` is a :class:`meshmode.discretization.connection.DiscretizationConnection` - objects going from the original mesh to the refined mesh. + going from the original mesh to the refined mesh. """ # Algorithm: @@ -507,17 +518,21 @@ def refine_for_global_qbx(lpot_source, code_container, # refinement. from meshmode.mesh.refinement import Refiner - refiner = Refiner(lpot_source.density_discr.mesh) - connections = [] + from meshmode.discretization.connection import ChainedDiscretizationConnection with cl.CommandQueue(lpot_source.cl_context) as queue: - wrangler = code_container.get_wrangler(queue, refiner) + wrangler = code_container.get_wrangler(queue) + + refiner = Refiner(lpot_source.density_discr.mesh) + connections = [] # Do initial refinement. if refine_flags is not None: - lpot_source, conn = wrangler.refine(lpot_source, refine_flags, - group_factory, debug) + conn = wrangler.refine( + lpot_source.density_discr, refiner, refine_flags, group_factory, + debug) connections.append(conn) + lpot_source = lpot_source.copy(density_discr=conn.to_discr) # {{{ first stage refinement @@ -541,7 +556,7 @@ def refine_for_global_qbx(lpot_source, code_container, refine_flags = make_empty_refine_flags(queue, lpot_source) # Check condition 1. - must_refine |= wrangler.check_center_is_closest_to_orig_panel( + must_refine |= wrangler.check_expansion_disks_undisturbed_by_sources( lpot_source, tree, peer_lists, refine_flags, debug) # Check condition 3. @@ -550,9 +565,11 @@ def refine_for_global_qbx(lpot_source, code_container, lpot_source, helmholtz_k, refine_flags, debug) if must_refine: - lpot_source, conn = wrangler.refine(lpot_source, refine_flags, + conn = wrangler.refine( + lpot_source.density_discr, refiner, refine_flags, group_factory, debug) connections.append(conn) + lpot_source = lpot_source.copy(density_discr=conn.to_discr) del tree del refine_flags @@ -562,11 +579,59 @@ def refine_for_global_qbx(lpot_source, code_container, # {{{ second stage refinement + del refiner + + must_refine = True + niter = 0 + fine_refiner = Refiner(lpot_source.fine_density_discr.mesh) + fine_connections = [lpot_source.resampler] + + while must_refine: + must_refine = False + niter += 1 + + if niter > maxiter: + logger.warning( + "Max iteration count reached in QBX layer potential source" + " refiner (second stage).") + break + + # Build tree and auxiliary data. + # FIXME: The tree should not have to be rebuilt at each iteration. + tree = wrangler.build_tree(lpot_source, use_fine_discr=True) + peer_lists = wrangler.find_peer_lists(tree) + refine_flags = make_empty_refine_flags( + queue, lpot_source, use_fine_discr=True) + + must_refine |= wrangler.check_sufficient_source_quadrature_resolution( + lpot_source, tree, peer_lists, refine_flags, debug) + + if must_refine: + conn = wrangler.refine( + lpot_source.fine_density_discr, + fine_refiner, refine_flags, fine_group_factory, debug) + fine_connections.append(conn) + lpot_source = lpot_source.copy( + fine_density_discr=conn.to_discr, + resampler=ChainedDiscretizationConnection(fine_connections)) + + del tree + del refine_flags + del peer_lists + # }}} lpot_source = lpot_source.copy(debug=debug, refined_for_global_qbx=True) - return lpot_source, connections + if len(connections) == 0: + from meshmode.discretization.connection import make_same_mesh_connection + connection = make_same_mesh_connection( + lpot_source.density_discr, + lpot_source.density_discr) + else: + connection = ChainedDiscretizationConnection(connections) + + return lpot_source, connection # vim: foldmethod=marker:filetype=pyopencl diff --git a/pytential/qbx/target_assoc.py b/pytential/qbx/target_assoc.py index 8890d992..b5f6b153 100644 --- a/pytential/qbx/target_assoc.py +++ b/pytential/qbx/target_assoc.py @@ -36,7 +36,7 @@ from boxtree.tools import DeviceDataRecord from boxtree.area_query import AreaQueryElementwiseTemplate from boxtree.tools import InlineBinarySearch from pytential.qbx.utils import ( - QBX_TREE_C_PREAMBLE, QBX_TREE_MAKO_DEFS, DiscrPlotterMixin) + QBX_TREE_C_PREAMBLE, QBX_TREE_MAKO_DEFS) unwrap_args = AreaQueryElementwiseTemplate.unwrap_args @@ -318,11 +318,11 @@ class QBXTargetAssociation(DeviceDataRecord): pass -class QBXTargetAssociator(DiscrPlotterMixin): +class QBXTargetAssociator(object): def __init__(self, cl_context): - from pytential.qbx.utils import TreeWithQBXMetadataBuilder - self.tree_builder = TreeWithQBXMetadataBuilder(cl_context) + from boxtree.tree_build import TreeBuilder + self.tree_builder = TreeBuilder(cl_context) self.cl_context = cl_context from boxtree.area_query import PeerListFinder, SpaceInvaderQueryBuilder self.peer_list_finder = PeerListFinder(cl_context) @@ -654,8 +654,14 @@ class QBXTargetAssociator(DiscrPlotterMixin): """ with cl.CommandQueue(self.cl_context) as queue: - tree = self.tree_builder(queue, lpot_source, + from pytential.qbx.utils import build_tree_with_qbx_metadata + + tree = build_tree_with_qbx_metadata( + queue, + self.tree_builder, + lpot_source, [discr for discr, _ in target_discrs_and_qbx_sides]) + peer_lists, evt = self.peer_list_finder(queue, tree, wait_for) wait_for = [evt] diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index 60da4b5b..076b4098 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -118,203 +118,220 @@ def get_interleaved_centers(queue, lpot_source): # }}} -# {{{ discr plotter mixin - -class DiscrPlotterMixin(object): - - def plot_discr(self, lpot_source, outfilename="discr.pdf"): - with cl.CommandQueue(self.cl_context) as queue: - tree = self.tree_builder(queue, lpot_source).get(queue=queue) - from boxtree.visualization import TreePlotter - - import matplotlib - matplotlib.use('Agg') - import matplotlib.pyplot as plt - tp = TreePlotter(tree) - tp.draw_tree() - sources = (tree.sources[0], tree.sources[1]) - sti = tree.sorted_target_ids - plt.plot(sources[0][sti[tree.qbx_user_source_slice]], - sources[1][sti[tree.qbx_user_source_slice]], - lw=0, marker=".", markersize=1, label="sources") - plt.plot(sources[0][sti[tree.qbx_user_center_slice]], - sources[1][sti[tree.qbx_user_center_slice]], - lw=0, marker=".", markersize=1, label="centers") - plt.plot(sources[0][sti[tree.qbx_user_target_slice]], - sources[1][sti[tree.qbx_user_target_slice]], - lw=0, marker=".", markersize=1, label="targets") - plt.axis("equal") - plt.legend() - plt.savefig(outfilename) +# {{{ discr plotter + +def plot_discr(lpot_source, outfilename="discr.pdf"): + with cl.CommandQueue(lpot_source.cl_context) as queue: + from boxtree.tree_builder import TreeBuilder + tree_builder = TreeBuilder(lpot_source.cl_context) + tree = tree_builder(queue, lpot_source).get(queue=queue) + from boxtree.visualization import TreePlotter + + import matplotlib + matplotlib.use('Agg') + import matplotlib.pyplot as plt + tp = TreePlotter(tree) + tp.draw_tree() + sources = (tree.sources[0], tree.sources[1]) + sti = tree.sorted_target_ids + plt.plot(sources[0][sti[tree.qbx_user_source_slice]], + sources[1][sti[tree.qbx_user_source_slice]], + lw=0, marker=".", markersize=1, label="sources") + plt.plot(sources[0][sti[tree.qbx_user_center_slice]], + sources[1][sti[tree.qbx_user_center_slice]], + lw=0, marker=".", markersize=1, label="centers") + plt.plot(sources[0][sti[tree.qbx_user_target_slice]], + sources[1][sti[tree.qbx_user_target_slice]], + lw=0, marker=".", markersize=1, label="targets") + plt.axis("equal") + plt.legend() + plt.savefig(outfilename) # }}} # {{{ tree creation -class TreeWithQBXMetadataBuilder(object): - MAX_REFINE_WEIGHT = 64 +class TreeWithQBXMetadata(Tree): + """ + .. attribute:: nqbxpanels + .. attribuet:: nqbxsources + .. attribute:: nqbxcenters + .. attribute:: nqbxtargets - class TreeWithQBXMetadata(Tree): - """ - .. attribute:: nqbxpanels - .. attribuet:: nqbxsources - .. attribute:: nqbxcenters - .. attribute:: nqbxtargets + .. attribute:: box_to_qbx_panel_starts + .. attribute:: box_to_qbx_panel_lists - .. attribute:: box_to_qbx_panel_starts - .. attribute:: box_to_qbx_panel_lists + .. attribute:: box_to_qbx_source_starts + .. attribute:: box_to_qbx_source_lists - .. attribute:: box_to_qbx_source_starts - .. attribute:: box_to_qbx_source_lists + .. attribute:: box_to_qbx_center_starts + .. attribute:: box_to_qbx_center_lists - .. attribute:: box_to_qbx_center_starts - .. attribute:: box_to_qbx_center_lists + .. attribute:: box_to_qbx_target_starts + .. attribute:: box_to_qbx_target_lists - .. attribute:: box_to_qbx_target_starts - .. attribute:: box_to_qbx_target_lists + .. attribute:: qbx_panel_to_source_starts + .. attribute:: qbx_panel_to_center_starts - .. attribute:: qbx_panel_to_source_starts - .. attribute:: qbx_panel_to_center_starts + .. attribute:: qbx_user_source_slice + .. attribute:: qbx_user_center_slice + .. attribute:: qbx_user_panel_slice + .. attribute:: qbx_user_target_slice + """ + pass - .. attribute:: qbx_user_source_slice - .. attribute:: qbx_user_center_slice - .. attribute:: qbx_user_panel_slice - .. attribute:: qbx_user_target_slice - """ - pass - def __init__(self, context): - self.context = context - from boxtree.tree_build import TreeBuilder - self.tree_builder = TreeBuilder(self.context) +MAX_REFINE_WEIGHT = 64 - def __call__(self, queue, lpot_source, targets_list=()): - """ - Return a :class:`TreeWithQBXMetadata` built from the given layer - potential source. - :arg queue: An instance of :class:`pyopencl.CommandQueue` +def build_tree_with_qbx_metadata( + queue, tree_builder, lpot_source, targets_list=(), + use_fine_discr=False): + """ + Return a :class:`TreeWithQBXMetadata` built from the given layer + potential source. - :arg lpot_source: An instance of - :class:`pytential.qbx.NewQBXLayerPotentialSource`. + :arg queue: An instance of :class:`pyopencl.CommandQueue` - :arg targets_list: A list of :class:`pytential.target.TargetBase` - """ - # The ordering of particles is as follows: - # - sources go first - # - then centers - # - then panels (=centers of mass) - # - then targets - - logger.info("start building tree with qbx metadata") - - sources = lpot_source.density_discr.nodes() - centers = get_interleaved_centers(queue, lpot_source) - centers_of_mass = lpot_source.panel_centers_of_mass() - targets = (tgt.nodes() for tgt in targets_list) - - particles = tuple( - cl.array.concatenate(dim_coords, queue=queue) - for dim_coords in zip(sources, centers, centers_of_mass, *targets)) - - # Counts - nparticles = len(particles[0]) - npanels = len(centers_of_mass[0]) - nsources = len(sources[0]) - ncenters = len(centers[0]) - # Each source gets an interior / exterior center. - assert 2 * nsources == ncenters - ntargets = sum(tgt.nnodes for tgt in targets_list) - - # Slices - qbx_user_source_slice = slice(0, nsources) - panel_slice_start = 3 * nsources - qbx_user_center_slice = slice(nsources, panel_slice_start) - target_slice_start = panel_slice_start + npanels - qbx_user_panel_slice = slice(panel_slice_start, panel_slice_start + npanels) - qbx_user_target_slice = slice(target_slice_start, - target_slice_start + ntargets) - - # Build tree with sources, centers, and centers of mass. Split boxes - # only because of sources. - refine_weights = cl.array.zeros(queue, nparticles, np.int32) - refine_weights[:nsources].fill(1) - - refine_weights.finish() - - tree, evt = self.tree_builder(queue, particles, - max_leaf_refine_weight=self.MAX_REFINE_WEIGHT, - refine_weights=refine_weights) - - # Compute box => particle class relations - flags = refine_weights - del refine_weights - particle_classes = {} - - for class_name, particle_slice, fixup in ( - ("box_to_qbx_source", qbx_user_source_slice, 0), - ("box_to_qbx_target", qbx_user_target_slice, -target_slice_start), - ("box_to_qbx_center", qbx_user_center_slice, -nsources), - ("box_to_qbx_panel", qbx_user_panel_slice, -panel_slice_start)): - flags.fill(0) - flags[particle_slice].fill(1) - flags.finish() - - from boxtree.tree import filter_target_lists_in_user_order - box_to_class = ( - filter_target_lists_in_user_order(queue, tree, flags) - .with_queue(queue)) - - if fixup: - box_to_class.target_lists += fixup - particle_classes[class_name + "_starts"] = box_to_class.target_starts - particle_classes[class_name + "_lists"] = box_to_class.target_lists - - del flags - del box_to_class - - # Compute panel => source relation - qbx_panel_to_source_starts = cl.array.empty( - queue, npanels + 1, dtype=tree.particle_id_dtype) - el_offset = 0 - for group in lpot_source.density_discr.groups: - qbx_panel_to_source_starts[el_offset:el_offset + group.nelements] = \ - cl.array.arange(queue, group.node_nr_base, - group.node_nr_base + group.nnodes, - group.nunit_nodes, - dtype=tree.particle_id_dtype) - el_offset += group.nelements - qbx_panel_to_source_starts[-1] = nsources - - # Compute panel => center relation - qbx_panel_to_center_starts = 2 * qbx_panel_to_source_starts - - # Transfer all tree attributes. - tree_attrs = {} - for attr_name in tree.__class__.fields: - try: - tree_attrs[attr_name] = getattr(tree, attr_name) - except AttributeError: - pass - - tree_attrs.update(particle_classes) - - logger.info("done building tree with qbx metadata") - - return self.TreeWithQBXMetadata( - qbx_panel_to_source_starts=qbx_panel_to_source_starts, - qbx_panel_to_center_starts=qbx_panel_to_center_starts, - qbx_user_source_slice=qbx_user_source_slice, - qbx_user_panel_slice=qbx_user_panel_slice, - qbx_user_center_slice=qbx_user_center_slice, - qbx_user_target_slice=qbx_user_target_slice, - nqbxpanels=npanels, - nqbxsources=nsources, - nqbxcenters=ncenters, - nqbxtargets=ntargets, - **tree_attrs).with_queue(None) + :arg lpot_source: An instance of + :class:`pytential.qbx.NewQBXLayerPotentialSource`. + + :arg targets_list: A list of :class:`pytential.target.TargetBase` + + :arg use_fine_discr: If *True*, + """ + # The ordering of particles is as follows: + # - sources go first + # - then centers + # - then panels (=centers of mass) + # - then targets + + logger.info("start building tree with qbx metadata") + + sources = ( + lpot_source.density_discr.nodes() + if not use_fine_discr + else lpot_source.fine_density_discr.nodes()) + + centers = get_interleaved_centers(queue, lpot_source) + + centers_of_mass = ( + lpot_source.panel_centers_of_mass() + if not use_fine_discr + else lpot_source.fine_panel_centers_of_mass()) + + targets = (tgt.nodes() for tgt in targets_list) + + particles = tuple( + cl.array.concatenate(dim_coords, queue=queue) + for dim_coords in zip(sources, centers, centers_of_mass, *targets)) + + # Counts + nparticles = len(particles[0]) + npanels = len(centers_of_mass[0]) + nsources = len(sources[0]) + ncenters = len(centers[0]) + # Each source gets an interior / exterior center. + assert 2 * nsources == ncenters or use_fine_discr + ntargets = sum(tgt.nnodes for tgt in targets_list) + + # Slices + qbx_user_source_slice = slice(0, nsources) + panel_slice_start = nsources + ncenters + qbx_user_center_slice = slice(nsources, panel_slice_start) + target_slice_start = panel_slice_start + npanels + qbx_user_panel_slice = slice(panel_slice_start, panel_slice_start + npanels) + qbx_user_target_slice = slice(target_slice_start, + target_slice_start + ntargets) + + # Build tree with sources, centers, and centers of mass. Split boxes + # only because of sources. + refine_weights = cl.array.zeros(queue, nparticles, np.int32) + refine_weights[:nsources].fill(1) + + refine_weights.finish() + + tree, evt = tree_builder(queue, particles, + max_leaf_refine_weight=MAX_REFINE_WEIGHT, + refine_weights=refine_weights) + + # Compute box => particle class relations + flags = refine_weights + del refine_weights + particle_classes = {} + + for class_name, particle_slice, fixup in ( + ("box_to_qbx_source", qbx_user_source_slice, 0), + ("box_to_qbx_target", qbx_user_target_slice, -target_slice_start), + ("box_to_qbx_center", qbx_user_center_slice, -nsources), + ("box_to_qbx_panel", qbx_user_panel_slice, -panel_slice_start)): + flags.fill(0) + flags[particle_slice].fill(1) + flags.finish() + + from boxtree.tree import filter_target_lists_in_user_order + box_to_class = ( + filter_target_lists_in_user_order(queue, tree, flags) + .with_queue(queue)) + + if fixup: + box_to_class.target_lists += fixup + particle_classes[class_name + "_starts"] = box_to_class.target_starts + particle_classes[class_name + "_lists"] = box_to_class.target_lists + + del flags + del box_to_class + + # Compute panel => source relation + if use_fine_discr: + density_discr = lpot_source.fine_density_discr + else: + density_discr = lpot_source.density_discr + + qbx_panel_to_source_starts = cl.array.empty( + queue, npanels + 1, dtype=tree.particle_id_dtype) + el_offset = 0 + for group in density_discr.groups: + qbx_panel_to_source_starts[el_offset:el_offset + group.nelements] = \ + cl.array.arange(queue, group.node_nr_base, + group.node_nr_base + group.nnodes, + group.nunit_nodes, + dtype=tree.particle_id_dtype) + el_offset += group.nelements + qbx_panel_to_source_starts[-1] = nsources + + # Compute panel => center relation + qbx_panel_to_center_starts = ( + 2 * qbx_panel_to_source_starts + if not use_fine_discr + else None) + + # Transfer all tree attributes. + tree_attrs = {} + for attr_name in tree.__class__.fields: + try: + tree_attrs[attr_name] = getattr(tree, attr_name) + except AttributeError: + pass + + tree_attrs.update(particle_classes) + + logger.info("done building tree with qbx metadata") + + return TreeWithQBXMetadata( + qbx_panel_to_source_starts=qbx_panel_to_source_starts, + qbx_panel_to_center_starts=qbx_panel_to_center_starts, + qbx_user_source_slice=qbx_user_source_slice, + qbx_user_panel_slice=qbx_user_panel_slice, + qbx_user_center_slice=qbx_user_center_slice, + qbx_user_target_slice=qbx_user_target_slice, + nqbxpanels=npanels, + nqbxsources=nsources, + nqbxcenters=ncenters, + nqbxtargets=ntargets, + **tree_attrs).with_queue(None) # }}} diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index 34e878ab..c7a27c92 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -108,7 +108,9 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory + factory = InterpolatoryQuadratureSimplexGroupFactory(order) + fine_factory = InterpolatoryQuadratureSimplexGroupFactory(4 * order) discr = Discretization(cl_ctx, mesh, factory) @@ -120,39 +122,39 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): lpot_source, conn = refine_for_global_qbx( lpot_source, RefinerCodeContainer(cl_ctx), - factory, helmholtz_k) + factory, fine_factory, helmholtz_k) discr_nodes = lpot_source.density_discr.nodes().get(queue) + fine_discr_nodes = lpot_source.fine_density_discr.nodes().get(queue) int_centers = lpot_source.centers(-1) int_centers = np.array([axis.get(queue) for axis in int_centers]) ext_centers = lpot_source.centers(+1) ext_centers = np.array([axis.get(queue) for axis in ext_centers]) panel_sizes = lpot_source.panel_sizes("npanels").get(queue) + fine_panel_sizes = lpot_source.fine_panel_sizes("npanels").get(queue) # }}} # {{{ check if satisfying criteria - def check_panel(panel): - # Check wavenumber to panel size ratio. - assert panel_sizes[panel.element_nr] * helmholtz_k <= 5 + def check_disk_undisturbed_by_sources(centers_panel, sources_panel): + h = panel_sizes[centers_panel.element_nr] - def check_panel_pair(panel_1, panel_2): - h_1 = panel_sizes[panel_1.element_nr] - h_2 = panel_sizes[panel_2.element_nr] - - if panel_1.element_nr == panel_2.element_nr: + if centers_panel.element_nr == sources_panel.element_nr: # Same panel return - panel_1_centers = int_centers[:, panel_1.discr_slice] - panel_2_nodes = discr_nodes[:, panel_2.discr_slice] + my_int_centers = int_centers[:, centers_panel.discr_slice] + my_ext_centers = ext_centers[:, centers_panel.discr_slice] + all_centers = np.append(my_int_centers, my_ext_centers, axis=-1) + + nodes = discr_nodes[:, sources_panel.discr_slice] - # =distance(centers of panel 1, panel 2) + # =distance(interior centers of panel 1, panel 2) dist = ( la.norm(( - panel_1_centers[..., np.newaxis] - - panel_2_nodes[:, np.newaxis, ...]).T, + all_centers[..., np.newaxis] - + nodes[:, np.newaxis, ...]).T, axis=-1) .min()) @@ -160,12 +162,43 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): # A center cannot be closer to another panel than to its originating # panel. - assert dist >= h_1 / 2, (dist, h_1, panel_1.element_nr, panel_2.element_nr) + assert dist >= h / 2, \ + (dist, h, centers_panel.element_nr, sources_panel.element_nr) + + def check_sufficient_quadrature_resolution(centers_panel, sources_panel): + h = fine_panel_sizes[sources_panel.element_nr] + + my_int_centers = int_centers[:, centers_panel.discr_slice] + my_ext_centers = ext_centers[:, centers_panel.discr_slice] + all_centers = np.append(my_int_centers, my_ext_centers, axis=-1) + + nodes = fine_discr_nodes[:, sources_panel.discr_slice] + + # =distance(interior centers of panel 1, panel 2) + dist = ( + la.norm(( + all_centers[..., np.newaxis] - + nodes[:, np.newaxis, ...]).T, + axis=-1) + .min()) + + # Criterion 2: + # A center cannot be closer to another panel than to its originating + # panel. + + assert dist >= h / 4, \ + (dist, h, centers_panel.element_nr, sources_panel.element_nr) + + def check_panel_size_to_helmholtz_k_ratio(panel): + # Check wavenumber to panel size ratio. + assert panel_sizes[panel.element_nr] * helmholtz_k <= 5 for panel_1 in iter_elements(lpot_source.density_discr): - check_panel(panel_1) for panel_2 in iter_elements(lpot_source.density_discr): - check_panel_pair(panel_1, panel_2) + check_disk_undisturbed_by_sources(panel_1, panel_2) + for panel_2 in iter_elements(lpot_source.fine_density_discr): + check_sufficient_quadrature_resolution(panel_1, panel_2) + check_panel_size_to_helmholtz_k_ratio(panel_1) # }}} -- GitLab From cd8e453c8411a9d8ec73cc2454ff81baa4514850 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 15 May 2017 15:51:52 -0500 Subject: [PATCH 03/22] [ci skip] Fix redundant title in documentation. --- pytential/qbx/refinement.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 8a939afb..3b3b1981 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -46,9 +46,6 @@ MAX_LEVELS_INCREMENT = 10 __doc__ = """ -Refinement -^^^^^^^^^^ - The refiner takes a layer potential source and refines it until it satisfies three global QBX refinement criteria: -- GitLab From 52e494641529c7a1736627412c7f4e63457041f4 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 15 May 2017 16:05:07 -0500 Subject: [PATCH 04/22] Refiner: Only oversample at the final level. --- pytential/qbx/__init__.py | 9 +++++---- pytential/qbx/refinement.py | 40 ++++++++++++++++++++++++++++--------- 2 files changed, 36 insertions(+), 13 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index ea96e87c..a9839be3 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -235,11 +235,11 @@ class QBXLayerPotentialSource(LayerPotentialSource): return self._fine_density_discr from meshmode.discretization.poly_element import ( - InterpolatoryQuadratureSimplexGroupFactory) + QuadratureSimplexGroupFactory) return Discretization( self.density_discr.cl_context, self.density_discr.mesh, - InterpolatoryQuadratureSimplexGroupFactory(self.fine_order), + QuadratureSimplexGroupFactory(self.fine_order), self.real_dtype) @property @@ -285,7 +285,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): from pytential.qbx.refinement import refine_for_global_qbx from meshmode.discretization.poly_element import ( - InterpolatoryQuadratureSimplexGroupFactory) + InterpolatoryQuadratureSimplexGroupFactory, + QuadratureSimplexGroupFactory) if target_order is None: target_order = self.density_discr.groups[0].order @@ -294,7 +295,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): self, self.refiner_code_container, InterpolatoryQuadratureSimplexGroupFactory(target_order), - InterpolatoryQuadratureSimplexGroupFactory(self.fine_order), + QuadratureSimplexGroupFactory(self.fine_order), maxiter=maxiter) return lpot, connection diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 3b3b1981..b46da151 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -482,7 +482,7 @@ def refine_for_global_qbx(lpot_source, code_container, :arg fine_group_factory: An instance of :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for - discretizing the fine refined mesh. + oversample the refined mesh at the finest level. :arg helmholtz_k: The Helmholtz parameter, or `None` if not applicable. @@ -501,21 +501,24 @@ def refine_for_global_qbx(lpot_source, code_container, # Algorithm: # - # 1. Do initial refinement, if requested. + # 1. Do an initial refinement, if requested. # # 2. While not converged: # Refine until each center is closest to its own # source and panel sizes bounded by wavelength. + # This becomes the coarse discretization. # # 3. While not converged: - # Refine fine density discretization until + # Refine density discretization until # sufficient quadrature resolution from all panels achieved. + # Oversample the result to get the fine density discretization. # # TODO: Stop doing redundant checks by avoiding panels which no longer need # refinement. from meshmode.mesh.refinement import Refiner - from meshmode.discretization.connection import ChainedDiscretizationConnection + from meshmode.discretization.connection import ( + ChainedDiscretizationConnection, make_same_mesh_connection) with cl.CommandQueue(lpot_source.cl_context) as queue: wrangler = code_container.get_wrangler(queue) @@ -578,10 +581,15 @@ def refine_for_global_qbx(lpot_source, code_container, del refiner + density_discr = lpot_source.density_discr + lpot_source = lpot_source.copy( + fine_density_discr=density_discr, + resampler=make_same_mesh_connection(density_discr, density_discr)) + must_refine = True niter = 0 - fine_refiner = Refiner(lpot_source.fine_density_discr.mesh) - fine_connections = [lpot_source.resampler] + fine_refiner = Refiner(lpot_source.density_discr.mesh) + fine_connections = [] while must_refine: must_refine = False @@ -606,7 +614,7 @@ def refine_for_global_qbx(lpot_source, code_container, if must_refine: conn = wrangler.refine( lpot_source.fine_density_discr, - fine_refiner, refine_flags, fine_group_factory, debug) + fine_refiner, refine_flags, group_factory, debug) fine_connections.append(conn) lpot_source = lpot_source.copy( fine_density_discr=conn.to_discr, @@ -618,10 +626,24 @@ def refine_for_global_qbx(lpot_source, code_container, # }}} - lpot_source = lpot_source.copy(debug=debug, refined_for_global_qbx=True) + # Oversample the fine mesh. + from meshmode.discretization import Discretization + ovsmp_fine_density_discr = Discretization( + lpot_source.fine_density_discr.cl_context, + lpot_source.fine_density_discr.mesh, + fine_group_factory, + lpot_source.real_dtype) + + fine_connections.append(make_same_mesh_connection( + lpot_source.fine_density_discr, ovsmp_fine_density_discr)) + + lpot_source = lpot_source.copy( + fine_density_discr=ovsmp_fine_density_discr, + resampler=ChainedDiscretizationConnection(fine_connections), + debug=debug, + refined_for_global_qbx=True) if len(connections) == 0: - from meshmode.discretization.connection import make_same_mesh_connection connection = make_same_mesh_connection( lpot_source.density_discr, lpot_source.density_discr) -- GitLab From ffaf7980069a972b2ac4a079b18929738bebfca7 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 15 May 2017 19:21:19 -0500 Subject: [PATCH 05/22] More fixes: * Get test_global_qbx working. * Do oversampling in the layer potential source. * Close #37. --- pytential/qbx/__init__.py | 83 ++++++++++++++++++++----------- pytential/qbx/refinement.py | 97 ++++++++++++------------------------- pytential/qbx/utils.py | 17 ++++--- test/test_global_qbx.py | 82 +++++++++++++++---------------- 4 files changed, 134 insertions(+), 145 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index a9839be3..8101e1fe 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -140,8 +140,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_order=None, fmm_level_to_order=None, target_stick_out_factor=1e-10, - fine_density_discr=None, - resampler=None, + base_fine_density_discr=None, + base_resampler=None, # begin undocumented arguments # FIXME default debug=False once everything works @@ -150,21 +150,36 @@ class QBXLayerPotentialSource(LayerPotentialSource): performance_data_file=None): """ :arg fine_order: The total degree to which the (upsampled) - underlying quadrature is exact. - :arg fine_density_discr: A discretization or *None.* - If non-*None*, should also supply *resampler.* - :arg resampler: A connection used for resampling. - If non-*None*, should also supply *fine_density_discr*. + underlying quadrature is exact. + :arg base_fine_density_discr: A discretization or *None.* + Represents the (non-upsampled) fine density discretization. + If *None*, defaults to *density_discr*. + If non-*None*, should also supply *resampler*. + :arg base_resampler: A connection used for resampling from *density_discr* + to *base_fine_density_discr*. + If non-*None*, should also supply *base_fine_density_discr*. :arg fmm_order: `False` for direct calculation. ``None`` will set a reasonable(-ish?) default. """ - if fine_density_discr is None and resampler is not None: + if base_fine_density_discr is None and base_resampler is not None: raise ValueError( - "resampler is supplied; must also supply fine_density_discr") - if fine_density_discr is not None and resampler is None: + "base_resampler is supplied; must also supply " + "base_fine_density_discr") + + if base_fine_density_discr is not None and base_resampler is None: raise ValueError( - "fine_density_discr is supplied; must also supply resampler") + "base_fine_density_discr is supplied; must also supply " + "base_resampler") + + if base_resampler is not None and base_fine_density_discr is not None: + if base_resampler.from_discr is not density_discr: + raise ValueError("density_discr is not the same as " + "base_resampler.from_discr") + + if base_resampler.to_discr is not base_fine_density_discr: + raise ValueError("base_fine_density_discr is not the same as " + "base_resampler.to_discr") if fmm_level_to_order is None: if fmm_order is None and qbx_order is not None: @@ -187,8 +202,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): self.target_stick_out_factor = target_stick_out_factor # Default values are lazily provided if these are None - self._fine_density_discr = fine_density_discr - self._resampler = resampler + self._base_fine_density_discr = base_fine_density_discr + self._base_resampler = base_resampler self.debug = debug self.refined_for_global_qbx = refined_for_global_qbx @@ -201,8 +216,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): qbx_order=None, fmm_level_to_order=None, target_stick_out_factor=None, - fine_density_discr=None, - resampler=None, + base_fine_density_discr=None, + base_resampler=None, debug=None, refined_for_global_qbx=None, @@ -218,8 +233,9 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_level_to_order or self.fmm_level_to_order), target_stick_out_factor=( target_stick_out_factor or self.target_stick_out_factor), - fine_density_discr=fine_density_discr or self._fine_density_discr, - resampler=resampler or self._resampler, + base_fine_density_discr=( + base_fine_density_discr or self._base_fine_density_discr), + base_resampler=base_resampler or self._base_resampler, debug=( debug if debug is not None else self.debug), @@ -231,26 +247,30 @@ class QBXLayerPotentialSource(LayerPotentialSource): @property @memoize_method def fine_density_discr(self): - if self._fine_density_discr is not None: - return self._fine_density_discr - from meshmode.discretization.poly_element import ( QuadratureSimplexGroupFactory) + base_fine_density_discr = ( + self._base_fine_density_discr or self.density_discr) + return Discretization( - self.density_discr.cl_context, self.density_discr.mesh, + self.density_discr.cl_context, base_fine_density_discr.mesh, QuadratureSimplexGroupFactory(self.fine_order), self.real_dtype) @property @memoize_method def resampler(self): - if self._resampler is not None: - return self._resampler + from meshmode.discretization.connection import ( + make_same_mesh_connection, ChainedDiscretizationConnection) + + conn = make_same_mesh_connection( + self.fine_density_discr, self._base_fine_density_discr) - from meshmode.discretization.connection import make_same_mesh_connection - return make_same_mesh_connection( - self.fine_density_discr, self.density_discr) + if self._base_resampler is not None: + return ChainedDiscretizationConnection([self._base_resampler, conn]) + + return conn def el_view(self, discr, group_nr, global_array): """Return a view of *global_array* of shape @@ -363,7 +383,9 @@ class QBXLayerPotentialSource(LayerPotentialSource): @memoize_method def fine_panel_centers_of_mass(self): - return self._centers_of_mass_for_discr(self.fine_density_discr) + base_fine_density_discr = ( + self._base_fine_density_discr or self.density_discr) + return self._centers_of_mass_for_discr(base_fine_density_discr) def _panel_sizes_for_discr(self, discr, last_dim_length): assert last_dim_length in ("nsources", "ncenters", "npanels") @@ -426,7 +448,12 @@ class QBXLayerPotentialSource(LayerPotentialSource): @memoize_method def fine_panel_sizes(self, last_dim_length="nsources"): - return self._panel_sizes_for_discr(self.fine_density_discr, last_dim_length) + if last_dim_length != "npanels": + raise NotImplementedError() + + base_fine_density_discr = ( + self._base_fine_density_discr or self.density_discr) + return self._panel_sizes_for_discr(base_fine_density_discr, last_dim_length) @memoize_method def centers(self, sign): diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index b46da151..29fca3e1 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -72,8 +72,8 @@ three global QBX refinement criteria: EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( extra_args=r""" /* input */ - particle_id_t *box_to_panel_starts, - particle_id_t *box_to_panel_lists, + particle_id_t *box_to_source_starts, + particle_id_t *box_to_source_lists, particle_id_t *panel_to_source_starts, particle_id_t *panel_to_center_starts, particle_id_t source_offset, @@ -92,39 +92,35 @@ EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( %endfor """, ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r""" + /* Find the panel associated with this center. */ particle_id_t my_panel = bsearch(panel_to_center_starts, npanels + 1, i); ${load_particle("INDEX_FOR_CENTER_PARTICLE(i)", ball_center)} ${ball_radius} = panel_sizes[my_panel] / 2; """, leaf_found_op=QBX_TREE_MAKO_DEFS + r""" - for (particle_id_t panel_idx = box_to_panel_starts[${leaf_box_id}]; - panel_idx < box_to_panel_starts[${leaf_box_id} + 1]; - ++panel_idx) + /* Check that each source in the leaf box is sufficiently far from the + center; if not, mark the panel for refinement. */ + + for (particle_id_t source_idx = box_to_source_starts[${leaf_box_id}]; + source_idx < box_to_source_starts[${leaf_box_id} + 1]; + ++source_idx) { - particle_id_t panel = box_to_panel_lists[panel_idx]; + particle_id_t source = box_to_source_lists[source_idx]; + particle_id_t panel = bsearch( + panel_to_source_starts, npanels + 1, source); - // Skip self. if (my_panel == panel) { continue; } - bool is_close = false; - - for (particle_id_t source = panel_to_source_starts[panel]; - source < panel_to_source_starts[panel + 1]; - ++source) - { - coord_vec_t source_coords; + coord_vec_t source_coords; + ${load_particle("INDEX_FOR_SOURCE_PARTICLE(source)", "source_coords")} - ${load_particle( - "INDEX_FOR_SOURCE_PARTICLE(source)", "source_coords")} - - is_close |= ( - distance(${ball_center}, source_coords) - <= panel_sizes[my_panel] / 2); - } + bool is_close = ( + distance(${ball_center}, source_coords) + <= panel_sizes[my_panel] / 2); if (is_close) { @@ -301,8 +297,8 @@ class RefinerWrangler(object): evt = knl( *unwrap_args( tree, peer_lists, - tree.box_to_qbx_panel_starts, - tree.box_to_qbx_panel_lists, + tree.box_to_qbx_source_starts, + tree.box_to_qbx_source_lists, tree.qbx_panel_to_source_starts, tree.qbx_panel_to_center_starts, tree.qbx_user_source_slice.start, @@ -347,7 +343,7 @@ class RefinerWrangler(object): if debug: npanels_to_refine_prev = cl.array.sum(refine_flags).get() - logger.info("refiner: checking adequate quadrature resolution") + logger.info("refiner: checking sufficient quadrature resolution") found_panel_to_refine = cl.array.zeros(self.queue, 1, np.int32) found_panel_to_refine.finish() @@ -383,7 +379,7 @@ class RefinerWrangler(object): logger.debug("refiner: found {} panel(s) to refine".format( npanels_to_refine - npanels_to_refine_prev)) - logger.info("refiner: done checking adequate quadrature resolution") + logger.info("refiner: done checking sufficient quadrature resolution") return found_panel_to_refine.get()[0] == 1 @@ -499,20 +495,6 @@ def refine_for_global_qbx(lpot_source, code_container, going from the original mesh to the refined mesh. """ - # Algorithm: - # - # 1. Do an initial refinement, if requested. - # - # 2. While not converged: - # Refine until each center is closest to its own - # source and panel sizes bounded by wavelength. - # This becomes the coarse discretization. - # - # 3. While not converged: - # Refine density discretization until - # sufficient quadrature resolution from all panels achieved. - # Oversample the result to get the fine density discretization. - # # TODO: Stop doing redundant checks by avoiding panels which no longer need # refinement. @@ -579,18 +561,12 @@ def refine_for_global_qbx(lpot_source, code_container, # {{{ second stage refinement - del refiner - - density_discr = lpot_source.density_discr - lpot_source = lpot_source.copy( - fine_density_discr=density_discr, - resampler=make_same_mesh_connection(density_discr, density_discr)) - must_refine = True niter = 0 - fine_refiner = Refiner(lpot_source.density_discr.mesh) fine_connections = [] + base_fine_density_discr = lpot_source.density_discr + while must_refine: must_refine = False niter += 1 @@ -613,12 +589,14 @@ def refine_for_global_qbx(lpot_source, code_container, if must_refine: conn = wrangler.refine( - lpot_source.fine_density_discr, - fine_refiner, refine_flags, group_factory, debug) + base_fine_density_discr, + refiner, refine_flags, group_factory, debug) + base_fine_density_discr = conn.to_discr fine_connections.append(conn) lpot_source = lpot_source.copy( - fine_density_discr=conn.to_discr, - resampler=ChainedDiscretizationConnection(fine_connections)) + base_fine_density_discr=base_fine_density_discr, + base_resampler=ChainedDiscretizationConnection( + fine_connections)) del tree del refine_flags @@ -626,22 +604,7 @@ def refine_for_global_qbx(lpot_source, code_container, # }}} - # Oversample the fine mesh. - from meshmode.discretization import Discretization - ovsmp_fine_density_discr = Discretization( - lpot_source.fine_density_discr.cl_context, - lpot_source.fine_density_discr.mesh, - fine_group_factory, - lpot_source.real_dtype) - - fine_connections.append(make_same_mesh_connection( - lpot_source.fine_density_discr, ovsmp_fine_density_discr)) - - lpot_source = lpot_source.copy( - fine_density_discr=ovsmp_fine_density_discr, - resampler=ChainedDiscretizationConnection(fine_connections), - debug=debug, - refined_for_global_qbx=True) + lpot_source = lpot_source.copy(debug=debug, refined_for_global_qbx=True) if len(connections) == 0: connection = make_same_mesh_connection( diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index 076b4098..87b11eff 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -239,12 +239,15 @@ def build_tree_with_qbx_metadata( # Slices qbx_user_source_slice = slice(0, nsources) - panel_slice_start = nsources + ncenters - qbx_user_center_slice = slice(nsources, panel_slice_start) - target_slice_start = panel_slice_start + npanels + + center_slice_start = nsources + qbx_user_center_slice = slice(center_slice_start, center_slice_start + ncenters) + + panel_slice_start = center_slice_start + ncenters qbx_user_panel_slice = slice(panel_slice_start, panel_slice_start + npanels) - qbx_user_target_slice = slice(target_slice_start, - target_slice_start + ntargets) + + target_slice_start = panel_slice_start + npanels + qbx_user_target_slice = slice(target_slice_start, target_slice_start + ntargets) # Build tree with sources, centers, and centers of mass. Split boxes # only because of sources. @@ -265,7 +268,7 @@ def build_tree_with_qbx_metadata( for class_name, particle_slice, fixup in ( ("box_to_qbx_source", qbx_user_source_slice, 0), ("box_to_qbx_target", qbx_user_target_slice, -target_slice_start), - ("box_to_qbx_center", qbx_user_center_slice, -nsources), + ("box_to_qbx_center", qbx_user_center_slice, -center_slice_start), ("box_to_qbx_panel", qbx_user_panel_slice, -panel_slice_start)): flags.fill(0) flags[particle_slice].fill(1) @@ -291,7 +294,7 @@ def build_tree_with_qbx_metadata( density_discr = lpot_source.density_discr qbx_panel_to_source_starts = cl.array.empty( - queue, npanels + 1, dtype=tree.particle_id_dtype) + queue, npanels + 1, dtype=tree.particle_id_dtype) el_offset = 0 for group in density_discr.groups: qbx_panel_to_source_starts[el_offset:el_offset + group.nelements] = \ diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index c7a27c92..c290530c 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -39,7 +39,7 @@ from pytential.qbx import QBXLayerPotentialSource from functools import partial from meshmode.mesh.generation import ( # noqa ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, - make_curve_mesh) + make_curve_mesh, generate_icosphere, generate_torus) from extra_curve_data import horseshoe @@ -58,59 +58,38 @@ FAR_TARGET_DIST_FROM_SOURCE = 10 class ElementInfo(RecordWithoutPickling): """ .. attribute:: element_nr - .. attribute:: neighbors .. attribute:: discr_slice - .. attribute:: mesh_slice - .. attribute:: element_group - .. attribute:: mesh_element_group """ __slots__ = ["element_nr", - "neighbors", "discr_slice"] def iter_elements(discr): discr_nodes_idx = 0 element_nr = 0 - adjacency = discr.mesh.nodal_adjacency for discr_group in discr.groups: start = element_nr for element_nr in range(start, start + discr_group.nelements): yield ElementInfo( element_nr=element_nr, - neighbors=list(adjacency.neighbors[ - slice(*adjacency.neighbors_starts[ - element_nr:element_nr+2])]), discr_slice=slice(discr_nodes_idx, discr_nodes_idx + discr_group.nunit_nodes)) discr_nodes_idx += discr_group.nunit_nodes -# }}} - -@pytest.mark.parametrize(("curve_name", "curve_f", "nelements"), [ - ("20-to-1 ellipse", partial(ellipse, 20), 100), - ("horseshoe", horseshoe, 64), - ]) -def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): +def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): cl_ctx = ctx_getter() queue = cl.CommandQueue(cl_ctx) - # {{{ generate lpot source, run refiner - - order = 16 - helmholtz_k = 10 - - mesh = make_curve_mesh(curve_f, np.linspace(0, 1, nelements+1), order) - from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory + from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory, + QuadratureSimplexGroupFactory) factory = InterpolatoryQuadratureSimplexGroupFactory(order) - fine_factory = InterpolatoryQuadratureSimplexGroupFactory(4 * order) + fine_factory = QuadratureSimplexGroupFactory(4 * order) discr = Discretization(cl_ctx, mesh, factory) @@ -133,8 +112,6 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): panel_sizes = lpot_source.panel_sizes("npanels").get(queue) fine_panel_sizes = lpot_source.fine_panel_sizes("npanels").get(queue) - # }}} - # {{{ check if satisfying criteria def check_disk_undisturbed_by_sources(centers_panel, sources_panel): @@ -150,7 +127,7 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): nodes = discr_nodes[:, sources_panel.discr_slice] - # =distance(interior centers of panel 1, panel 2) + # =distance(centers of panel 1, panel 2) dist = ( la.norm(( all_centers[..., np.newaxis] - @@ -158,7 +135,7 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): axis=-1) .min()) - # Criterion 1: + # Criterion: # A center cannot be closer to another panel than to its originating # panel. @@ -182,10 +159,9 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): axis=-1) .min()) - # Criterion 2: - # A center cannot be closer to another panel than to its originating - # panel. - + # Criterion: + # The quadrature contribution from each panel is as accurate + # as from the center's own source panel. assert dist >= h / 4, \ (dist, h, centers_panel.element_nr, sources_panel.element_nr) @@ -193,15 +169,40 @@ def test_source_refinement(ctx_getter, curve_name, curve_f, nelements): # Check wavenumber to panel size ratio. assert panel_sizes[panel.element_nr] * helmholtz_k <= 5 - for panel_1 in iter_elements(lpot_source.density_discr): + for i, panel_1 in enumerate(iter_elements(lpot_source.density_discr)): for panel_2 in iter_elements(lpot_source.density_discr): check_disk_undisturbed_by_sources(panel_1, panel_2) for panel_2 in iter_elements(lpot_source.fine_density_discr): check_sufficient_quadrature_resolution(panel_1, panel_2) - check_panel_size_to_helmholtz_k_ratio(panel_1) + if helmholtz_k is not None: + check_panel_size_to_helmholtz_k_ratio(panel_1) # }}} +# }}} + + +@pytest.mark.parametrize(("curve_name", "curve_f", "nelements"), [ + ("20-to-1 ellipse", partial(ellipse, 20), 100), + ("horseshoe", horseshoe, 64), + ]) +def test_source_refinement_2d(ctx_getter, curve_name, curve_f, nelements): + # {{{ generate lpot source, run refiner + helmholtz_k = 10 + order = 8 + + mesh = make_curve_mesh(curve_f, np.linspace(0, 1, nelements+1), order) + run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k) + + +@pytest.mark.parametrize(("surface_name", "surface_f", "order"), [ + ("sphere", partial(generate_icosphere, 1), 4), + ("torus", partial(generate_torus, 3, 1, n_inner=10, n_outer=7), 6), + ]) +def test_source_refinement_3d(ctx_getter, surface_name, surface_f, order): + mesh = surface_f(order=order) + run_source_refinement_test(ctx_getter, mesh, order) + @pytest.mark.parametrize(("curve_name", "curve_f", "nelements"), [ ("20-to-1 ellipse", partial(ellipse, 20), 100), @@ -225,13 +226,8 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements): discr = Discretization(cl_ctx, mesh, factory) - from pytential.qbx.refinement import QBXLayerPotentialSourceRefiner - - lpot_source = QBXLayerPotentialSource(discr, order) + lpot_source, conn = QBXLayerPotentialSource(discr, order).with_refinement() del discr - refiner = QBXLayerPotentialSourceRefiner(cl_ctx) - - lpot_source, conn = refiner(lpot_source, factory) int_centers = lpot_source.centers(-1) ext_centers = lpot_source.centers(+1) -- GitLab From 943e1785a7bebe19656444d451eaa7270b9339db Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 15 May 2017 19:31:41 -0500 Subject: [PATCH 06/22] Fix from_discr. --- pytential/qbx/__init__.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 8101e1fe..1f9cd60d 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -264,8 +264,11 @@ class QBXLayerPotentialSource(LayerPotentialSource): from meshmode.discretization.connection import ( make_same_mesh_connection, ChainedDiscretizationConnection) + base_fine_density_discr = ( + self._base_fine_density_discr or self.density_discr) + conn = make_same_mesh_connection( - self.fine_density_discr, self._base_fine_density_discr) + self.fine_density_discr, base_fine_density_discr) if self._base_resampler is not None: return ChainedDiscretizationConnection([self._base_resampler, conn]) -- GitLab From eb4919c96b07733b747ceb5af97927162fee4a3d Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 18 May 2017 12:58:39 -0500 Subject: [PATCH 07/22] Refiner: Don't specify refinement conditions by their numbers. --- pytential/qbx/refinement.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 29fca3e1..3fbddf09 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -68,7 +68,6 @@ three global QBX refinement criteria: # {{{ kernels -# Refinement checker for Condition 1. EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( extra_args=r""" /* input */ @@ -134,7 +133,6 @@ EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( preamble=str(InlineBinarySearch("particle_id_t"))) -# Refinement checker for Condition 2. SUFFICIENT_SOURCE_QUADRATURE_RESOLUTION_CHECKER = AreaQueryElementwiseTemplate( extra_args=r""" /* input */ @@ -537,11 +535,9 @@ def refine_for_global_qbx(lpot_source, code_container, peer_lists = wrangler.find_peer_lists(tree) refine_flags = make_empty_refine_flags(queue, lpot_source) - # Check condition 1. must_refine |= wrangler.check_expansion_disks_undisturbed_by_sources( lpot_source, tree, peer_lists, refine_flags, debug) - # Check condition 3. if helmholtz_k is not None: must_refine |= wrangler.check_helmholtz_k_to_panel_size_ratio( lpot_source, helmholtz_k, refine_flags, debug) -- GitLab From 4f513216b4da484e196ea5e0068fa9c4d03f92ce Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 18 May 2017 13:00:37 -0500 Subject: [PATCH 08/22] QBXLayerPotentialSource: Only use base_resampler. --- pytential/qbx/__init__.py | 64 ++++++++++--------------------------- pytential/qbx/refinement.py | 1 - 2 files changed, 17 insertions(+), 48 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 1f9cd60d..f4954bb3 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -140,7 +140,6 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_order=None, fmm_level_to_order=None, target_stick_out_factor=1e-10, - base_fine_density_discr=None, base_resampler=None, # begin undocumented arguments @@ -151,36 +150,15 @@ class QBXLayerPotentialSource(LayerPotentialSource): """ :arg fine_order: The total degree to which the (upsampled) underlying quadrature is exact. - :arg base_fine_density_discr: A discretization or *None.* - Represents the (non-upsampled) fine density discretization. - If *None*, defaults to *density_discr*. - If non-*None*, should also supply *resampler*. - :arg base_resampler: A connection used for resampling from *density_discr* - to *base_fine_density_discr*. - If non-*None*, should also supply *base_fine_density_discr*. + :arg base_resampler: A connection used for resampling from + *density_discr* the fine density discretization. It is assumed + that the fine density discretization given by + *base_resampler.to_discr* is *not* already upsampled. May + be *None*. :arg fmm_order: `False` for direct calculation. ``None`` will set - a reasonable(-ish?) default. + a reasonable(-ish?) default. """ - if base_fine_density_discr is None and base_resampler is not None: - raise ValueError( - "base_resampler is supplied; must also supply " - "base_fine_density_discr") - - if base_fine_density_discr is not None and base_resampler is None: - raise ValueError( - "base_fine_density_discr is supplied; must also supply " - "base_resampler") - - if base_resampler is not None and base_fine_density_discr is not None: - if base_resampler.from_discr is not density_discr: - raise ValueError("density_discr is not the same as " - "base_resampler.from_discr") - - if base_resampler.to_discr is not base_fine_density_discr: - raise ValueError("base_fine_density_discr is not the same as " - "base_resampler.to_discr") - if fmm_level_to_order is None: if fmm_order is None and qbx_order is not None: fmm_order = qbx_order + 1 @@ -202,7 +180,6 @@ class QBXLayerPotentialSource(LayerPotentialSource): self.target_stick_out_factor = target_stick_out_factor # Default values are lazily provided if these are None - self._base_fine_density_discr = base_fine_density_discr self._base_resampler = base_resampler self.debug = debug @@ -216,7 +193,6 @@ class QBXLayerPotentialSource(LayerPotentialSource): qbx_order=None, fmm_level_to_order=None, target_stick_out_factor=None, - base_fine_density_discr=None, base_resampler=None, debug=None, @@ -233,8 +209,6 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_level_to_order or self.fmm_level_to_order), target_stick_out_factor=( target_stick_out_factor or self.target_stick_out_factor), - base_fine_density_discr=( - base_fine_density_discr or self._base_fine_density_discr), base_resampler=base_resampler or self._base_resampler, debug=( @@ -244,17 +218,20 @@ class QBXLayerPotentialSource(LayerPotentialSource): else self.refined_for_global_qbx), performance_data_file=self.performance_data_file) + @property + def _base_fine_density_discr(self): + return (self._base_resampler.to_discr + if self._base_resampler is not None + else self.density_discr) + @property @memoize_method def fine_density_discr(self): from meshmode.discretization.poly_element import ( QuadratureSimplexGroupFactory) - base_fine_density_discr = ( - self._base_fine_density_discr or self.density_discr) - return Discretization( - self.density_discr.cl_context, base_fine_density_discr.mesh, + self.density_discr.cl_context, self._base_fine_density_discr.mesh, QuadratureSimplexGroupFactory(self.fine_order), self.real_dtype) @@ -264,11 +241,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): from meshmode.discretization.connection import ( make_same_mesh_connection, ChainedDiscretizationConnection) - base_fine_density_discr = ( - self._base_fine_density_discr or self.density_discr) - conn = make_same_mesh_connection( - self.fine_density_discr, base_fine_density_discr) + self.fine_density_discr, self._base_fine_density_discr) if self._base_resampler is not None: return ChainedDiscretizationConnection([self._base_resampler, conn]) @@ -386,9 +360,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): @memoize_method def fine_panel_centers_of_mass(self): - base_fine_density_discr = ( - self._base_fine_density_discr or self.density_discr) - return self._centers_of_mass_for_discr(base_fine_density_discr) + return self._centers_of_mass_for_discr(self._base_fine_density_discr) def _panel_sizes_for_discr(self, discr, last_dim_length): assert last_dim_length in ("nsources", "ncenters", "npanels") @@ -453,10 +425,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): def fine_panel_sizes(self, last_dim_length="nsources"): if last_dim_length != "npanels": raise NotImplementedError() - - base_fine_density_discr = ( - self._base_fine_density_discr or self.density_discr) - return self._panel_sizes_for_discr(base_fine_density_discr, last_dim_length) + return self._panel_sizes_for_discr( + self._base_fine_density_discr, last_dim_length) @memoize_method def centers(self, sign): diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 3fbddf09..3bb295d3 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -590,7 +590,6 @@ def refine_for_global_qbx(lpot_source, code_container, base_fine_density_discr = conn.to_discr fine_connections.append(conn) lpot_source = lpot_source.copy( - base_fine_density_discr=base_fine_density_discr, base_resampler=ChainedDiscretizationConnection( fine_connections)) -- GitLab From 1f1db7ca51ab3147085fb5ed7f996b071baedc89 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 18 May 2017 13:10:40 -0500 Subject: [PATCH 09/22] Fix accidental maxiter increase. --- pytential/qbx/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index f4954bb3..a8302ca8 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -272,7 +272,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): return RefinerCodeContainer(self.cl_context) @memoize_method - def with_refinement(self, target_order=None, maxiter=10): + def with_refinement(self, target_order=None, maxiter=3): """ :returns: a tuple ``(lpot_src, cnx)``, where ``lpot_src`` is a :class:`QBXLayerPotentialSource` and ``cnx`` is a -- GitLab From 5c28185fe979bf5c23c55d8157cfb203dec95eb2 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 18 May 2017 13:47:50 -0500 Subject: [PATCH 10/22] Revert "Refiner: Don't specify refinement conditions by their numbers." This reverts commit eb4919c96b07733b747ceb5af97927162fee4a3d. --- pytential/qbx/refinement.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 3bb295d3..3db3b753 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -68,6 +68,7 @@ three global QBX refinement criteria: # {{{ kernels +# Refinement checker for Condition 1. EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( extra_args=r""" /* input */ @@ -133,6 +134,7 @@ EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( preamble=str(InlineBinarySearch("particle_id_t"))) +# Refinement checker for Condition 2. SUFFICIENT_SOURCE_QUADRATURE_RESOLUTION_CHECKER = AreaQueryElementwiseTemplate( extra_args=r""" /* input */ @@ -535,9 +537,11 @@ def refine_for_global_qbx(lpot_source, code_container, peer_lists = wrangler.find_peer_lists(tree) refine_flags = make_empty_refine_flags(queue, lpot_source) + # Check condition 1. must_refine |= wrangler.check_expansion_disks_undisturbed_by_sources( lpot_source, tree, peer_lists, refine_flags, debug) + # Check condition 3. if helmholtz_k is not None: must_refine |= wrangler.check_helmholtz_k_to_panel_size_ratio( lpot_source, helmholtz_k, refine_flags, debug) -- GitLab From 222aa25b772f655460d19de43aa199893dade1fe Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 18 May 2017 14:53:02 -0500 Subject: [PATCH 11/22] Document TreeWithQBXMetadata. --- pytential/qbx/utils.py | 53 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 52 insertions(+), 1 deletion(-) diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index 87b11eff..08eb1a1d 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -154,27 +154,78 @@ def plot_discr(lpot_source, outfilename="discr.pdf"): class TreeWithQBXMetadata(Tree): - """ + """A subclass of :class:`boxtree.tree.Tree`. Has all of that class's + attributes, along with the following: + .. attribute:: nqbxpanels .. attribuet:: nqbxsources .. attribute:: nqbxcenters .. attribute:: nqbxtargets + .. ------------------------------------------------------------------------ + .. rubric:: Box properties + .. ------------------------------------------------------------------------ + + Box to QBX panels + ---- + .. attribute:: box_to_qbx_panel_starts + + ``box_id_t [nboxes + 1]`` + .. attribute:: box_to_qbx_panel_lists + ``particle_id_t [*]`` + + Box to QBX sources + ---- + .. attribute:: box_to_qbx_source_starts + + ``box_id_t [nboxes + 1]`` + .. attribute:: box_to_qbx_source_lists + ``particle_id_t [*]`` + + Box to QBX centers + ---- + .. attribute:: box_to_qbx_center_starts + + ``box_id_t [nboxes + 1]`` + .. attribute:: box_to_qbx_center_lists + ``particle_id_t [*]`` + + Box to QBX targets + ---- + .. attribute:: box_to_qbx_target_starts + + ``box_id_t [nboxes + 1]`` + .. attribute:: box_to_qbx_target_lists + ``particle_id_t [*]`` + + .. ------------------------------------------------------------------------ + .. rubric:: Panel properties + .. ------------------------------------------------------------------------ + .. attribute:: qbx_panel_to_source_starts + + ``particle_id_t [nqbxpanels + 1]`` + .. attribute:: qbx_panel_to_center_starts + ``particle_id_t [nqbxpanels + 1]`` + + .. ------------------------------------------------------------------------ + .. rubric:: Particle order indices + .. ------------------------------------------------------------------------ + .. attribute:: qbx_user_source_slice .. attribute:: qbx_user_center_slice .. attribute:: qbx_user_panel_slice -- GitLab From 976b88d207199158e42cdfd981c5f9e4d2251d9d Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 18 May 2017 16:42:16 -0500 Subject: [PATCH 12/22] TreeWithQBXMetadata docs: Change title to rubric. --- pytential/qbx/utils.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index 08eb1a1d..b1ecd4ad 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -166,8 +166,7 @@ class TreeWithQBXMetadata(Tree): .. rubric:: Box properties .. ------------------------------------------------------------------------ - Box to QBX panels - ---- + .. rubric:: Box to QBX panels .. attribute:: box_to_qbx_panel_starts @@ -177,8 +176,7 @@ class TreeWithQBXMetadata(Tree): ``particle_id_t [*]`` - Box to QBX sources - ---- + .. rubric:: Box to QBX sources .. attribute:: box_to_qbx_source_starts @@ -188,8 +186,7 @@ class TreeWithQBXMetadata(Tree): ``particle_id_t [*]`` - Box to QBX centers - ---- + .. rubric:: Box to QBX centers .. attribute:: box_to_qbx_center_starts @@ -199,8 +196,7 @@ class TreeWithQBXMetadata(Tree): ``particle_id_t [*]`` - Box to QBX targets - ---- + .. rubric:: Box to QBX targets .. attribute:: box_to_qbx_target_starts -- GitLab From 804166c8d428b94569cab336a4144b15506398fe Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 21 May 2017 16:48:49 -0500 Subject: [PATCH 13/22] Refiner: Make failure to converge a silenceable error. --- pytential/__init__.py | 7 +++++++ pytential/qbx/refinement.py | 21 +++++++++++++++------ 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/pytential/__init__.py b/pytential/__init__.py index 98e62c71..7a4e1f4e 100644 --- a/pytential/__init__.py +++ b/pytential/__init__.py @@ -48,7 +48,14 @@ def _set_up_logging_from_environment(): set_up_logging(pytential_log_var.split(":"), level=level) +def _set_up_errors(): + import warnings + from pytential.qbx.refinement import RefinerNotConvergedWarning + warnings.filterwarnings("error", category=RefinerNotConvergedWarning) + + _set_up_logging_from_environment() +_set_up_errors() @memoize_on_first_arg diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 3db3b753..4a419baa 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -61,6 +61,7 @@ three global QBX refinement criteria: .. autoclass:: RefinerCodeContainer .. autoclass:: RefinerWrangler +.. autoclass:: RefinerNotConvergedWarning .. autofunction:: make_empty_refine_flags .. autofunction:: refine_for_global_qbx @@ -444,6 +445,10 @@ class RefinerWrangler(object): # }}} +class RefinerNotConvergedWarning(UserWarning): + pass + + def make_empty_refine_flags(queue, lpot_source, use_fine_discr=False): """Return an array on the device suitable for use as element refine flags. @@ -526,9 +531,11 @@ def refine_for_global_qbx(lpot_source, code_container, niter += 1 if niter > maxiter: - logger.warning( - "Max iteration count reached in QBX layer potential source" - " refiner (first stage).") + from warnings import warn + warn( + "Max iteration count reached in QBX layer potential source" + " refiner.", + RefinerNotConvergedWarning) break # Build tree and auxiliary data. @@ -572,9 +579,11 @@ def refine_for_global_qbx(lpot_source, code_container, niter += 1 if niter > maxiter: - logger.warning( - "Max iteration count reached in QBX layer potential source" - " refiner (second stage).") + from warnings import warn + warn( + "Max iteration count reached in QBX layer potential source" + " refiner.", + RefinerNotConvergedWarning) break # Build tree and auxiliary data. -- GitLab From 1ec4c18dd879f96f1dbb75026e0ad1fe097b71ed Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 21 May 2017 17:02:16 -0500 Subject: [PATCH 14/22] Refiner: Change helmholtz_k to kernel_length_scale. Also, add and start using a kernel_length_scale parameter to with_refinement(). (closes #42). --- pytential/qbx/__init__.py | 4 +++- pytential/qbx/refinement.py | 35 +++++++++++++++++++---------------- test/test_global_qbx.py | 2 +- test/test_layer_pot.py | 20 ++++++++++++++++---- 4 files changed, 39 insertions(+), 22 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index a8302ca8..46921700 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -272,7 +272,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): return RefinerCodeContainer(self.cl_context) @memoize_method - def with_refinement(self, target_order=None, maxiter=3): + def with_refinement(self, target_order=None, kernel_length_scale=None, + maxiter=3): """ :returns: a tuple ``(lpot_src, cnx)``, where ``lpot_src`` is a :class:`QBXLayerPotentialSource` and ``cnx`` is a @@ -293,6 +294,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): self.refiner_code_container, InterpolatoryQuadratureSimplexGroupFactory(target_order), QuadratureSimplexGroupFactory(self.fine_order), + kernel_length_scale=kernel_length_scale, maxiter=maxiter) return lpot, connection diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 4a419baa..50ac9f41 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -56,8 +56,9 @@ three global QBX refinement criteria: The quadrature contribution from each panel is as accurate as from the center's own source panel. - * *Condition 3* (Panel size bounded based on wavelength) - (Helmholtz only) The panel size is bounded with respect to the wavelength. + * *Condition 3* (Panel size bounded based on kernel length scale) + The panel size is bounded by a kernel length scale. This + applies only to Helmholtz kernels. .. autoclass:: RefinerCodeContainer .. autoclass:: RefinerWrangler @@ -224,12 +225,12 @@ class RefinerCodeContainer(object): extra_type_aliases=(("particle_id_t", particle_id_dtype),)) @memoize_method - def helmholtz_k_to_panel_size_ratio_checker(self): + def kernel_length_scale_to_panel_size_ratio_checker(self): knl = lp.make_kernel( "{[panel]: 0<=panel oversize = panel_sizes[panel] * helmholtz_k > 5 + <> oversize = panel_sizes[panel] > kernel_length_scale if oversize refine_flags[panel] = 1 refine_flags_updated = 1 {id=write_refine_flags_updated} @@ -238,7 +239,7 @@ class RefinerCodeContainer(object): """, options="return_dict", silenced_warnings="write_race(write_refine_flags_updated)", - name="refine_helmholtz_k_to_panel_size_ratio") + name="refine_kernel_length_scale_to_panel_size_ratio") knl = lp.split_iname(knl, "panel", 128, inner_tag="l.0", outer_tag="g.0") return knl @@ -384,11 +385,11 @@ class RefinerWrangler(object): return found_panel_to_refine.get()[0] == 1 - def check_helmholtz_k_to_panel_size_ratio(self, lpot_source, - helmholtz_k, refine_flags, debug, wait_for=None): - knl = self.code_container.helmholtz_k_to_panel_size_ratio_checker() + def check_kernel_length_scale_to_panel_size_ratio(self, lpot_source, + kernel_length_scale, refine_flags, debug, wait_for=None): + knl = self.code_container.kernel_length_scale_to_panel_size_ratio_checker() - logger.info("refiner: checking helmholtz k to panel size ratio") + logger.info("refiner: checking kernel length scale to panel size ratio") if debug: npanels_to_refine_prev = cl.array.sum(refine_flags).get() @@ -397,7 +398,7 @@ class RefinerWrangler(object): panel_sizes=lpot_source.panel_sizes("npanels"), refine_flags=refine_flags, refine_flags_updated=np.array(0), - helmholtz_k=np.array(helmholtz_k), + kernel_length_scale=np.array(kernel_length_scale), wait_for=wait_for) cl.wait_for_events([evt]) @@ -408,7 +409,7 @@ class RefinerWrangler(object): logger.debug("refiner: found {} panel(s) to refine".format( npanels_to_refine - npanels_to_refine_prev)) - logger.info("refiner: done checking helmholtz k to panel size ratio") + logger.info("refiner: done checking kernel length scale to panel size ratio") return (out["refine_flags_updated"].get() == 1).all() @@ -467,7 +468,7 @@ def make_empty_refine_flags(queue, lpot_source, use_fine_discr=False): def refine_for_global_qbx(lpot_source, code_container, - group_factory, fine_group_factory, helmholtz_k=None, + group_factory, fine_group_factory, kernel_length_scale=None, # FIXME: Set debug=False once everything works. refine_flags=None, debug=True, maxiter=50): """ @@ -485,7 +486,8 @@ def refine_for_global_qbx(lpot_source, code_container, :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for oversample the refined mesh at the finest level. - :arg helmholtz_k: The Helmholtz parameter, or `None` if not applicable. + :arg kernel_length_scale: The kernel length scale, or *None* if not + applicable. All panels are refined to below this size. :arg refine_flags: A :class:`pyopencl.array.Array` indicating which panels should get refined initially, or `None` if no initial @@ -549,9 +551,10 @@ def refine_for_global_qbx(lpot_source, code_container, lpot_source, tree, peer_lists, refine_flags, debug) # Check condition 3. - if helmholtz_k is not None: - must_refine |= wrangler.check_helmholtz_k_to_panel_size_ratio( - lpot_source, helmholtz_k, refine_flags, debug) + if kernel_length_scale is not None: + must_refine |= ( + wrangler.check_kernel_length_scale_to_panel_size_ratio( + lpot_source, kernel_length_scale, refine_flags, debug)) if must_refine: conn = wrangler.refine( diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index c290530c..d7b4fb92 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -101,7 +101,7 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): lpot_source, conn = refine_for_global_qbx( lpot_source, RefinerCodeContainer(cl_ctx), - factory, fine_factory, helmholtz_k) + factory, fine_factory, kernel_length_scale=5*helmholtz_k) discr_nodes = lpot_source.density_discr.nodes().get(queue) fine_discr_nodes = lpot_source.fine_density_discr.nodes().get(queue) diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 1c704e52..15713a5e 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -309,10 +309,15 @@ def run_int_eq_test( if source_order is None: source_order = 4*target_order + refiner_extra_kwargs = {} + + if k != 0: + refiner_extra_kwargs["kernel_length_scale"] = 5*k + qbx, _ = QBXLayerPotentialSource( pre_density_discr, fine_order=source_order, qbx_order=qbx_order, # Don't use FMM for now - fmm_order=False).with_refinement() + fmm_order=False).with_refinement(**refiner_extra_kwargs) density_discr = qbx.density_discr @@ -798,9 +803,16 @@ def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, elif d == 3: order_bump = 6 - qbx, _ = QBXLayerPotentialSource( - pre_density_discr, 4*target_order, - qbx_order, fmm_order=qbx_order + order_bump).with_refinement() + refiner_extra_kwargs = {} + + if k != 0: + refiner_extra_kwargs["kernel_length_scale"] = 5*k + + qbx, _ = ( + QBXLayerPotentialSource( + pre_density_discr, 4*target_order, + qbx_order, fmm_order=qbx_order + order_bump) + .with_refinement(**refiner_extra_kwargs)) density_discr = qbx.density_discr # {{{ compute values of a solution to the PDE -- GitLab From 7ae32a2a8313ea59ec5b4c056174ea9cc5d63cde Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 21 May 2017 17:09:34 -0500 Subject: [PATCH 15/22] LayerPotentialSource.with_refinement(): Increase default maxiter to 10. --- pytential/qbx/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 46921700..936fc232 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -273,7 +273,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): @memoize_method def with_refinement(self, target_order=None, kernel_length_scale=None, - maxiter=3): + maxiter=10): """ :returns: a tuple ``(lpot_src, cnx)``, where ``lpot_src`` is a :class:`QBXLayerPotentialSource` and ``cnx`` is a -- GitLab From b9a22b6df5b65a6c17fb705cebefa42d28bd312c Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 21 May 2017 17:13:21 -0500 Subject: [PATCH 16/22] Refiner: Improve variable name and usage for center/source danger zone radii. --- pytential/qbx/refinement.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 50ac9f41..94f59206 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -81,7 +81,7 @@ EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( particle_id_t source_offset, particle_id_t center_offset, particle_id_t *sorted_target_ids, - coord_t *panel_sizes, + coord_t *center_danger_zone_radii_by_panel, int npanels, /* output */ @@ -98,7 +98,7 @@ EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( particle_id_t my_panel = bsearch(panel_to_center_starts, npanels + 1, i); ${load_particle("INDEX_FOR_CENTER_PARTICLE(i)", ball_center)} - ${ball_radius} = panel_sizes[my_panel] / 2; + ${ball_radius} = center_danger_zone_radii_by_panel[my_panel]; """, leaf_found_op=QBX_TREE_MAKO_DEFS + r""" /* Check that each source in the leaf box is sufficiently far from the @@ -122,7 +122,7 @@ EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( bool is_close = ( distance(${ball_center}, source_coords) - <= panel_sizes[my_panel] / 2); + <= center_danger_zone_radii_by_panel[my_panel]); if (is_close) { @@ -146,7 +146,7 @@ SUFFICIENT_SOURCE_QUADRATURE_RESOLUTION_CHECKER = AreaQueryElementwiseTemplate( particle_id_t source_offset, particle_id_t center_offset, particle_id_t *sorted_target_ids, - coord_t *ball_radii_by_panel, + coord_t *source_danger_zone_radii_by_panel, int npanels, /* output */ @@ -163,7 +163,7 @@ SUFFICIENT_SOURCE_QUADRATURE_RESOLUTION_CHECKER = AreaQueryElementwiseTemplate( particle_id_t my_panel = bsearch(panel_to_source_starts, npanels + 1, i); ${load_particle("INDEX_FOR_SOURCE_PARTICLE(i)", ball_center)} - ${ball_radius} = ball_radii_by_panel[my_panel]; + ${ball_radius} = source_danger_zone_radii_by_panel[my_panel]; """, leaf_found_op=QBX_TREE_MAKO_DEFS + r""" /* Check that each center in the leaf box is sufficiently far from the @@ -181,7 +181,7 @@ SUFFICIENT_SOURCE_QUADRATURE_RESOLUTION_CHECKER = AreaQueryElementwiseTemplate( bool is_close = ( distance(${ball_center}, center_coords) - <= ball_radii_by_panel[my_panel]); + <= source_danger_zone_radii_by_panel[my_panel]); if (is_close) { @@ -296,6 +296,10 @@ class RefinerWrangler(object): found_panel_to_refine.finish() unwrap_args = AreaQueryElementwiseTemplate.unwrap_args + center_danger_zone_radii_by_panel = ( + lpot_source.panel_sizes("npanels") + .with_queue(self.queue) / 2) + evt = knl( *unwrap_args( tree, peer_lists, @@ -306,7 +310,7 @@ class RefinerWrangler(object): tree.qbx_user_source_slice.start, tree.qbx_user_center_slice.start, tree.sorted_target_ids, - lpot_source.panel_sizes("npanels"), + center_danger_zone_radii_by_panel, tree.nqbxpanels, refine_flags, found_panel_to_refine, @@ -350,7 +354,7 @@ class RefinerWrangler(object): found_panel_to_refine = cl.array.zeros(self.queue, 1, np.int32) found_panel_to_refine.finish() - ball_radii_by_panel = ( + source_danger_zone_radii_by_panel = ( lpot_source.fine_panel_sizes("npanels") .with_queue(self.queue) / 4) unwrap_args = AreaQueryElementwiseTemplate.unwrap_args @@ -364,7 +368,7 @@ class RefinerWrangler(object): tree.qbx_user_source_slice.start, tree.qbx_user_center_slice.start, tree.sorted_target_ids, - ball_radii_by_panel, + source_danger_zone_radii_by_panel, tree.nqbxpanels, refine_flags, found_panel_to_refine, -- GitLab From b84edfebf82801c395dc5e19c5065ea1240563b7 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 21 May 2017 17:42:18 -0500 Subject: [PATCH 17/22] Fix handling of helmholtz_k parameter in refiner tests. --- test/test_global_qbx.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index d7b4fb92..fc09dc83 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -99,9 +99,13 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): lpot_source = QBXLayerPotentialSource(discr, order) del discr + refiner_extra_kwargs = {} + if helmholtz_k is not None: + refiner_extra_kwargs["kernel_length_scale"] = 5*helmholtz_k + lpot_source, conn = refine_for_global_qbx( lpot_source, RefinerCodeContainer(cl_ctx), - factory, fine_factory, kernel_length_scale=5*helmholtz_k) + factory, fine_factory, **refiner_extra_kwargs) discr_nodes = lpot_source.density_discr.nodes().get(queue) fine_discr_nodes = lpot_source.fine_density_discr.nodes().get(queue) -- GitLab From dbd7c1eb8f53d394363068e0065b38a6025b43c8 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 21 May 2017 18:10:52 -0500 Subject: [PATCH 18/22] Fix Helmholtz k max panel size: 5*k -> 5/k --- test/test_global_qbx.py | 2 +- test/test_layer_pot.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index fc09dc83..5071ce84 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -101,7 +101,7 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): refiner_extra_kwargs = {} if helmholtz_k is not None: - refiner_extra_kwargs["kernel_length_scale"] = 5*helmholtz_k + refiner_extra_kwargs["kernel_length_scale"] = 5/helmholtz_k lpot_source, conn = refine_for_global_qbx( lpot_source, RefinerCodeContainer(cl_ctx), diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 15713a5e..870cc6f0 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -312,7 +312,7 @@ def run_int_eq_test( refiner_extra_kwargs = {} if k != 0: - refiner_extra_kwargs["kernel_length_scale"] = 5*k + refiner_extra_kwargs["kernel_length_scale"] = 5/k qbx, _ = QBXLayerPotentialSource( pre_density_discr, fine_order=source_order, qbx_order=qbx_order, @@ -806,7 +806,7 @@ def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, refiner_extra_kwargs = {} if k != 0: - refiner_extra_kwargs["kernel_length_scale"] = 5*k + refiner_extra_kwargs["kernel_length_scale"] = 5/k qbx, _ = ( QBXLayerPotentialSource( -- GitLab From 6db0166c7a10cbfba4b9a7a6f778a30b89ccae59 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 25 May 2017 13:03:57 -0500 Subject: [PATCH 19/22] Refinement: Refine over the (non-upsampled) base_fine_density_discr, instead of using the upsampled version. --- pytential/qbx/__init__.py | 12 +++++++----- pytential/qbx/refinement.py | 15 +++++++-------- pytential/qbx/utils.py | 28 +++++++++++++++++----------- 3 files changed, 31 insertions(+), 24 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 936fc232..aec6b28f 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -219,7 +219,9 @@ class QBXLayerPotentialSource(LayerPotentialSource): performance_data_file=self.performance_data_file) @property - def _base_fine_density_discr(self): + def base_fine_density_discr(self): + """The fine density discretization before upsampling is applied. + """ return (self._base_resampler.to_discr if self._base_resampler is not None else self.density_discr) @@ -231,7 +233,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): QuadratureSimplexGroupFactory) return Discretization( - self.density_discr.cl_context, self._base_fine_density_discr.mesh, + self.density_discr.cl_context, self.base_fine_density_discr.mesh, QuadratureSimplexGroupFactory(self.fine_order), self.real_dtype) @@ -242,7 +244,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): make_same_mesh_connection, ChainedDiscretizationConnection) conn = make_same_mesh_connection( - self.fine_density_discr, self._base_fine_density_discr) + self.fine_density_discr, self.base_fine_density_discr) if self._base_resampler is not None: return ChainedDiscretizationConnection([self._base_resampler, conn]) @@ -362,7 +364,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): @memoize_method def fine_panel_centers_of_mass(self): - return self._centers_of_mass_for_discr(self._base_fine_density_discr) + return self._centers_of_mass_for_discr(self.base_fine_density_discr) def _panel_sizes_for_discr(self, discr, last_dim_length): assert last_dim_length in ("nsources", "ncenters", "npanels") @@ -428,7 +430,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): if last_dim_length != "npanels": raise NotImplementedError() return self._panel_sizes_for_discr( - self._base_fine_density_discr, last_dim_length) + self.base_fine_density_discr, last_dim_length) @memoize_method def centers(self, sign): diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 94f59206..9ae54bb6 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -256,7 +256,6 @@ class RefinerCodeContainer(object): def get_wrangler(self, queue): """ :arg queue: - :arg refiner: """ return RefinerWrangler(self, queue) @@ -417,11 +416,11 @@ class RefinerWrangler(object): return (out["refine_flags_updated"].get() == 1).all() - def build_tree(self, lpot_source, use_fine_discr=False): + def build_tree(self, lpot_source, use_base_fine_discr=False): tb = self.code_container.tree_builder() from pytential.qbx.utils import build_tree_with_qbx_metadata return build_tree_with_qbx_metadata( - self.queue, tb, lpot_source, use_fine_discr=use_fine_discr) + self.queue, tb, lpot_source, use_base_fine_discr=use_base_fine_discr) def find_peer_lists(self, tree): plf = self.code_container.peer_list_finder() @@ -454,7 +453,7 @@ class RefinerNotConvergedWarning(UserWarning): pass -def make_empty_refine_flags(queue, lpot_source, use_fine_discr=False): +def make_empty_refine_flags(queue, lpot_source, use_base_fine_discr=False): """Return an array on the device suitable for use as element refine flags. :arg queue: An instance of :class:`pyopencl.CommandQueue`. @@ -463,8 +462,8 @@ def make_empty_refine_flags(queue, lpot_source, use_fine_discr=False): :returns: A :class:`pyopencl.array.Array` suitable for use as refine flags, initialized to zero. """ - discr = (lpot_source.fine_density_discr - if use_fine_discr + discr = (lpot_source.base_fine_density_discr + if use_base_fine_discr else lpot_source.density_discr) result = cl.array.zeros(queue, discr.mesh.nelements, np.int32) result.finish() @@ -595,10 +594,10 @@ def refine_for_global_qbx(lpot_source, code_container, # Build tree and auxiliary data. # FIXME: The tree should not have to be rebuilt at each iteration. - tree = wrangler.build_tree(lpot_source, use_fine_discr=True) + tree = wrangler.build_tree(lpot_source, use_base_fine_discr=True) peer_lists = wrangler.find_peer_lists(tree) refine_flags = make_empty_refine_flags( - queue, lpot_source, use_fine_discr=True) + queue, lpot_source, use_base_fine_discr=True) must_refine |= wrangler.check_sufficient_source_quadrature_resolution( lpot_source, tree, peer_lists, refine_flags, debug) diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index b1ecd4ad..d39dc292 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -235,10 +235,14 @@ MAX_REFINE_WEIGHT = 64 def build_tree_with_qbx_metadata( queue, tree_builder, lpot_source, targets_list=(), - use_fine_discr=False): - """ - Return a :class:`TreeWithQBXMetadata` built from the given layer - potential source. + use_base_fine_discr=False): + """Return a :class:`TreeWithQBXMetadata` built from the given layer + potential source. This contains particles of four different types: + + * source particles and panel centers of mass either from + ``lpot_source.density_discr`` or ``lpot_source.base_fine_density_discr`` + * centers from ``lpot_source.centers()`` + * targets from ``targets_list``. :arg queue: An instance of :class:`pyopencl.CommandQueue` @@ -247,9 +251,11 @@ def build_tree_with_qbx_metadata( :arg targets_list: A list of :class:`pytential.target.TargetBase` - :arg use_fine_discr: If *True*, + :arg use_base_fine_discr: If *True*, builds a tree with sources/centers of + mass from ``lpot_source.base_fine_density_discr``. If *False* (default), + they are from ``lpot_source.density_discr``. """ - # The ordering of particles is as follows: + # the ordering of particles is as follows: # - sources go first # - then centers # - then panels (=centers of mass) @@ -259,14 +265,14 @@ def build_tree_with_qbx_metadata( sources = ( lpot_source.density_discr.nodes() - if not use_fine_discr + if not use_base_fine_discr else lpot_source.fine_density_discr.nodes()) centers = get_interleaved_centers(queue, lpot_source) centers_of_mass = ( lpot_source.panel_centers_of_mass() - if not use_fine_discr + if not use_base_fine_discr else lpot_source.fine_panel_centers_of_mass()) targets = (tgt.nodes() for tgt in targets_list) @@ -281,7 +287,7 @@ def build_tree_with_qbx_metadata( nsources = len(sources[0]) ncenters = len(centers[0]) # Each source gets an interior / exterior center. - assert 2 * nsources == ncenters or use_fine_discr + assert 2 * nsources == ncenters or use_base_fine_discr ntargets = sum(tgt.nnodes for tgt in targets_list) # Slices @@ -335,7 +341,7 @@ def build_tree_with_qbx_metadata( del box_to_class # Compute panel => source relation - if use_fine_discr: + if use_base_fine_discr: density_discr = lpot_source.fine_density_discr else: density_discr = lpot_source.density_discr @@ -355,7 +361,7 @@ def build_tree_with_qbx_metadata( # Compute panel => center relation qbx_panel_to_center_starts = ( 2 * qbx_panel_to_source_starts - if not use_fine_discr + if not use_base_fine_discr else None) # Transfer all tree attributes. -- GitLab From 7b64d69e0b26639c538d039a74ad5cad5a5c473c Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 25 May 2017 13:08:55 -0500 Subject: [PATCH 20/22] [ci skip] Fixes to comments. --- pytential/qbx/utils.py | 2 +- test/test_global_qbx.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py index d39dc292..8b363250 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -255,7 +255,7 @@ def build_tree_with_qbx_metadata( mass from ``lpot_source.base_fine_density_discr``. If *False* (default), they are from ``lpot_source.density_discr``. """ - # the ordering of particles is as follows: + # The ordering of particles is as follows: # - sources go first # - then centers # - then panels (=centers of mass) diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index 5071ce84..27bdea91 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -155,7 +155,7 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): nodes = fine_discr_nodes[:, sources_panel.discr_slice] - # =distance(interior centers of panel 1, panel 2) + # =distance(centers of panel 1, panel 2) dist = ( la.norm(( all_centers[..., np.newaxis] - -- GitLab From e105bab5e8e7797c0f6353502efe11ee4a144132 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 25 May 2017 13:11:04 -0500 Subject: [PATCH 21/22] [ci skip] More fixes to comments. --- test/test_global_qbx.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index 27bdea91..4013b8f3 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -53,7 +53,7 @@ RNG_SEED = 10 FAR_TARGET_DIST_FROM_SOURCE = 10 -# {{{ utilities for iterating over panels +# {{{ source refinement checker class ElementInfo(RecordWithoutPickling): """ @@ -191,7 +191,6 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): ("horseshoe", horseshoe, 64), ]) def test_source_refinement_2d(ctx_getter, curve_name, curve_f, nelements): - # {{{ generate lpot source, run refiner helmholtz_k = 10 order = 8 -- GitLab From aacfa01e22d53d357b2355446d81b40a1113f24f Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 25 May 2017 19:11:53 -0500 Subject: [PATCH 22/22] Change default argument of panel_sizes() to 'npanels'. --- pytential/qbx/__init__.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index aec6b28f..9c9b4ea3 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -422,11 +422,11 @@ class QBXLayerPotentialSource(LayerPotentialSource): return panel_sizes.with_queue(None) @memoize_method - def panel_sizes(self, last_dim_length="nsources"): + def panel_sizes(self, last_dim_length="npanels"): return self._panel_sizes_for_discr(self.density_discr, last_dim_length) @memoize_method - def fine_panel_sizes(self, last_dim_length="nsources"): + def fine_panel_sizes(self, last_dim_length="npanels"): if last_dim_length != "npanels": raise NotImplementedError() return self._panel_sizes_for_discr( @@ -441,7 +441,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): with cl.CommandQueue(self.cl_context) as queue: nodes = bind(self.density_discr, sym.nodes(adim))(queue) normals = bind(self.density_discr, sym.normal(adim, dim=dim))(queue) - panel_sizes = self.panel_sizes().with_queue(queue) + panel_sizes = self.panel_sizes("nsources").with_queue(queue) return (nodes + normals * sign * panel_sizes / 2).as_vector(np.object) @memoize_method -- GitLab