diff --git a/doc/qbx.rst b/doc/qbx.rst index d297db5eeb540327dc159c0b02e0316f86f5939b..46244da1086be469757fed4251ebdd946ed83af3 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/__init__.py b/pytential/__init__.py index 98e62c7164105a4c8cb0d1225706ad666ce233d1..7a4e1f4ec91fed3ae83e72396a9e703271a5908b 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/__init__.py b/pytential/qbx/__init__.py index a18196f401ce355cd5c42f2df0416cf5bc08eb88..9c9b4ea399f56895c10b80003752c56f6fe95e98 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,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_order=None, fmm_level_to_order=None, target_stick_out_factor=1e-10, + base_resampler=None, # begin undocumented arguments # FIXME default debug=False once everything works @@ -149,9 +149,14 @@ class QBXLayerPotentialSource(LayerPotentialSource): performance_data_file=None): """ :arg fine_order: The total degree to which the (upsampled) - underlying quadrature is exact. + underlying quadrature is exact. + :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 fmm_level_to_order is None: @@ -174,6 +179,9 @@ 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._base_resampler = base_resampler + self.debug = debug self.refined_for_global_qbx = refined_for_global_qbx self.performance_data_file = performance_data_file @@ -185,6 +193,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): qbx_order=None, fmm_level_to_order=None, target_stick_out_factor=None, + base_resampler=None, debug=None, refined_for_global_qbx=None, @@ -200,6 +209,7 @@ 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_resampler=base_resampler or self._base_resampler, debug=( debug if debug is not None else self.debug), @@ -208,19 +218,38 @@ class QBXLayerPotentialSource(LayerPotentialSource): else self.refined_for_global_qbx), performance_data_file=self.performance_data_file) + @property + 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) + @property @memoize_method def fine_density_discr(self): + from meshmode.discretization.poly_element import ( + QuadratureSimplexGroupFactory) + return Discretization( - self.density_discr.cl_context, self.density_discr.mesh, - QuadratureSimplexGroupFactory(self.fine_order), self.real_dtype) + self.density_discr.cl_context, self.base_fine_density_discr.mesh, + QuadratureSimplexGroupFactory(self.fine_order), + self.real_dtype) @property @memoize_method def resampler(self): - from meshmode.discretization.connection import make_same_mesh_connection - return make_same_mesh_connection( - self.fine_density_discr, self.density_discr) + from meshmode.discretization.connection import ( + make_same_mesh_connection, ChainedDiscretizationConnection) + + conn = make_same_mesh_connection( + self.fine_density_discr, self.base_fine_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 @@ -238,23 +267,38 @@ class QBXLayerPotentialSource(LayerPotentialSource): global_array.shape[:-1] + (group.nelements,)) + @property @memoize_method - def with_refinement(self, target_order=None, maxiter=3): + 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, kernel_length_scale=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, + QuadratureSimplexGroupFactory) + 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), + QuadratureSimplexGroupFactory(self.fine_order), + kernel_length_scale=kernel_length_scale, maxiter=maxiter) + return lpot, connection @property @@ -284,8 +328,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 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 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} @@ -389,21 +239,47 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin): """, 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 - # }}} + @memoize_method + def tree_builder(self): + from boxtree.tree_build import TreeBuilder + return TreeBuilder(self.cl_context) - # {{{ refinement triggering + @memoize_method + def peer_list_finder(self): + from boxtree.area_query import PeerListFinder + return PeerListFinder(self.cl_context) + + def get_wrangler(self, queue): + """ + :arg queue: + """ + return RefinerWrangler(self, queue) + +# }}} + + +# {{{ wrangler + +class RefinerWrangler(object): + + def __init__(self, code_container, queue): + self.code_container = code_container + self.queue = queue + + def check_expansion_disks_undisturbed_by_sources(self, + lpot_source, tree, peer_lists, refine_flags, + debug, wait_for=None): - def refinement_check_center_is_closest_to_orig_panel(self, queue, tree, - lpot_source, peer_lists, 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_closest_to_orig_panel_refiner( + knl = self.code_container.expansion_disk_undisturbed_by_sources_checker( tree.dimensions, tree.coord_dtype, tree.box_id_dtype, peer_lists.peer_list_starts.dtype, @@ -415,26 +291,32 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin): 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() + 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, - 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, 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, *tree.sources), range=slice(tree.nqbxcenters), - queue=queue) + queue=self.queue, + wait_for=wait_for) cl.wait_for_events([evt]) @@ -448,51 +330,51 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin): return found_panel_to_refine.get()[0] == 1 - def refinement_check_center_is_far_from_nonneighbor_panels(self, queue, - tree, lpot_source, peer_lists, tq_dists, refine_flags, debug, - wait_for=None): + def check_sufficient_source_quadrature_resolution( + self, lpot_source, tree, peer_lists, 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.sufficient_source_quadrature_resolution_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") - if debug: npanels_to_refine_prev = cl.array.sum(refine_flags).get() - found_panel_to_refine = cl.array.zeros(queue, 1, np.int32) + logger.info("refiner: checking sufficient quadrature resolution") + + 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) + source_danger_zone_radii_by_panel = ( + lpot_source.fine_panel_sizes("npanels") + .with_queue(self.queue) / 4) + unwrap_args = AreaQueryElementwiseTemplate.unwrap_args evt = knl( *unwrap_args( tree, peer_lists, - tree.box_to_qbx_panel_starts, - tree.box_to_qbx_panel_lists, + tree.box_to_qbx_center_starts, + tree.box_to_qbx_center_lists, tree.qbx_panel_to_source_starts, - 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, + source_danger_zone_radii_by_panel, tree.nqbxpanels, - tq_dists, refine_flags, found_panel_to_refine, *tree.sources), - range=slice(tree.nqbxcenters), - queue=queue) + range=slice(tree.nqbxsources), + queue=self.queue, + wait_for=wait_for) cl.wait_for_events([evt]) @@ -502,24 +384,24 @@ 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 sufficient quadrature resolution") return found_panel_to_refine.get()[0] == 1 - def refinement_check_helmholtz_k_to_panel_size_ratio(self, queue, lpot_source, - helmholtz_k, refine_flags, debug, wait_for=None): - knl = self.get_helmholtz_k_to_panel_ratio_refiner() + 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() - 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), - helmholtz_k=np.array(helmholtz_k), + kernel_length_scale=np.array(kernel_length_scale), wait_for=wait_for) cl.wait_for_events([evt]) @@ -530,246 +412,222 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin): 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() - 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 + 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_base_fine_discr=use_base_fine_discr) - 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 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 get_tunnel_query_dists(self, queue, tree, lpot_source): + def refine(self, density_discr, refiner, refine_flags, factory, debug): """ - Compute radii for the tubular neighborhood around each panel center of mass. + Refine the underlying mesh and discretization. """ - 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)) + if isinstance(refine_flags, cl.array.Array): + refine_flags = refine_flags.get(self.queue) + refine_flags = refine_flags.astype(np.bool) - cl.wait_for_events([evt]) + logger.info("refiner: calling meshmode") - if tree.coord_dtype != tq_dists.dtype: - tq_dists = tq_dists.astype(tree.coord_dtype) + refiner.refine(refine_flags) + from meshmode.discretization.connection import make_refinement_connection + conn = make_refinement_connection(refiner, density_discr, factory) - return tq_dists, evt + logger.info("refiner: done calling meshmode") - 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 conn - 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`. +class RefinerNotConvergedWarning(UserWarning): + pass - :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 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. - def refine(self, queue, lpot_source, refine_flags, refiner, 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.astype(np.bool) + :arg queue: An instance of :class:`pyopencl.CommandQueue`. + :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`. - logger.info("refiner: calling meshmode") + :returns: A :class:`pyopencl.array.Array` suitable for use as refine flags, + initialized to zero. + """ + 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() + return result + + +def refine_for_global_qbx(lpot_source, code_container, + group_factory, fine_group_factory, kernel_length_scale=None, + # FIXME: Set debug=False once everything works. + refine_flags=None, debug=True, maxiter=50): + """ + Entry point for calling the refiner. - refiner.refine(refine_flags) - from meshmode.discretization.connection import make_refinement_connection + :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`. - conn = make_refinement_connection( - refiner, lpot_source.density_discr, - factory) + :arg code_container: An instance of :class:`RefinerCodeContainer`. - logger.info("refiner: done calling meshmode") + :arg group_factory: An instance of + :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for + discretizing the coarse refined mesh. - new_density_discr = conn.to_discr + :arg fine_group_factory: An instance of + :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for + oversample the refined mesh at the finest level. - new_lpot_source = lpot_source.copy( - density_discr=new_density_discr, - refined_for_global_qbx=True) + :arg kernel_length_scale: The kernel length scale, or *None* if not + applicable. All panels are refined to below this size. - return new_lpot_source, conn + :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()`. - 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 maxiter: The maximum number of refiner iterations. - :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`. + :returns: A tuple ``(lpot_source, *conn*)`` where ``lpot_source`` is the + refined layer potential source, and ``conn`` is a + :class:`meshmode.discretization.connection.DiscretizationConnection` + going from the original mesh to the refined mesh. + """ - :arg group_factory: An instance of - :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for - discretizing the refined mesh. + # TODO: Stop doing redundant checks by avoiding panels which no longer need + # refinement. - :arg helmholtz_k: The Helmholtz parameter, or `None` if not applicable. + from meshmode.mesh.refinement import Refiner + from meshmode.discretization.connection import ( + ChainedDiscretizationConnection, make_same_mesh_connection) - :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()`. + with cl.CommandQueue(lpot_source.cl_context) as queue: + wrangler = code_container.get_wrangler(queue) - :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) + # Do initial refinement. + if refine_flags is not None: + 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) - 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) + # {{{ first stage refinement - done_refining = False - niter = 0 + must_refine = True + niter = 0 - while not done_refining: - niter += 1 + while must_refine: + must_refine = False + niter += 1 - if niter > maxiter: - logger.warning( + if niter > maxiter: + from warnings import warn + warn( "Max iteration count reached in QBX layer potential source" - " refiner.") - break + " refiner.", + RefinerNotConvergedWarning) + break - # 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 = [] + # 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) + + # Check condition 1. + must_refine |= wrangler.check_expansion_disks_undisturbed_by_sources( + lpot_source, tree, peer_lists, refine_flags, debug) + + # Check condition 3. + 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( + lpot_source.density_discr, refiner, refine_flags, + group_factory, debug) + connections.append(conn) + lpot_source = lpot_source.copy(density_discr=conn.to_discr) - peer_lists, evt = self.peer_list_finder(queue, tree, wait_for) - wait_for = [evt] + del tree + del refine_flags + del peer_lists - refine_flags, evt = self.get_refine_flags(queue, lpot_source) - wait_for.append(evt) + # }}} - tq_dists, evt = self.get_tunnel_query_dists(queue, tree, lpot_source) - wait_for.append(evt) + # {{{ second stage refinement - # Run refinement checkers. - must_refine = False + must_refine = True + niter = 0 + fine_connections = [] - must_refine |= \ - self.refinement_check_center_is_closest_to_orig_panel( - queue, tree, lpot_source, peer_lists, refine_flags, - debug, wait_for) + base_fine_density_discr = lpot_source.density_discr - must_refine |= \ - self.refinement_check_center_is_far_from_nonneighbor_panels( - queue, tree, lpot_source, peer_lists, tq_dists, - refine_flags, debug, wait_for) + while must_refine: + must_refine = False + niter += 1 - must_refine |= \ - self.refinement_check_2_to_1_panel_size_ratio( - queue, lpot_source, refine_flags, debug, wait_for) + if niter > maxiter: + from warnings import warn + warn( + "Max iteration count reached in QBX layer potential source" + " refiner.", + RefinerNotConvergedWarning) + break - if helmholtz_k: - must_refine |= \ - self.refinement_check_helmholtz_k_to_panel_size_ratio( - queue, lpot_source, helmholtz_k, refine_flags, debug, - wait_for) + # Build tree and auxiliary data. + # FIXME: The tree should not have to be rebuilt at each iteration. + 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_base_fine_discr=True) - if must_refine: - lpot_source, conn = self.refine( - queue, lpot_source, refine_flags, refiner, discr_factory, - debug) - connections.append(conn) + must_refine |= wrangler.check_sufficient_source_quadrature_resolution( + lpot_source, tree, peer_lists, refine_flags, debug) - del tree - del peer_lists - del tq_dists - del refine_flags - done_refining = not must_refine + if must_refine: + conn = wrangler.refine( + 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( + base_resampler=ChainedDiscretizationConnection( + fine_connections)) - return lpot_source, connections + del tree + del refine_flags + del peer_lists + + # }}} + + lpot_source = lpot_source.copy(debug=debug, refined_for_global_qbx=True) + + if len(connections) == 0: + 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 8890d992106e09509ceee424d73f8fa13ec9fa4f..b5f6b153e2e340f8e4f2315be2543017a469fe2f 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 60da4b5b2a4e914696ba0ea1083c0cc1a6f825d0..8b363250d644bfaebb65d738208ebab057c677aa 100644 --- a/pytential/qbx/utils.py +++ b/pytential/qbx/utils.py @@ -118,203 +118,276 @@ 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): + """A subclass of :class:`boxtree.tree.Tree`. Has all of that class's + attributes, along with the following: - class TreeWithQBXMetadata(Tree): - """ - .. attribute:: nqbxpanels - .. attribuet:: nqbxsources - .. attribute:: nqbxcenters - .. attribute:: nqbxtargets + .. attribute:: nqbxpanels + .. attribuet:: nqbxsources + .. attribute:: nqbxcenters + .. attribute:: nqbxtargets - .. attribute:: box_to_qbx_panel_starts - .. attribute:: box_to_qbx_panel_lists + .. ------------------------------------------------------------------------ + .. rubric:: Box properties + .. ------------------------------------------------------------------------ - .. attribute:: box_to_qbx_source_starts - .. attribute:: box_to_qbx_source_lists + .. rubric:: Box to QBX panels - .. attribute:: box_to_qbx_center_starts - .. attribute:: box_to_qbx_center_lists + .. attribute:: box_to_qbx_panel_starts - .. attribute:: box_to_qbx_target_starts - .. attribute:: box_to_qbx_target_lists + ``box_id_t [nboxes + 1]`` - .. attribute:: qbx_panel_to_source_starts - .. attribute:: qbx_panel_to_center_starts + .. attribute:: box_to_qbx_panel_lists - .. attribute:: qbx_user_source_slice - .. attribute:: qbx_user_center_slice - .. attribute:: qbx_user_panel_slice - .. attribute:: qbx_user_target_slice - """ - pass + ``particle_id_t [*]`` - def __init__(self, context): - self.context = context - from boxtree.tree_build import TreeBuilder - self.tree_builder = TreeBuilder(self.context) + .. rubric:: Box to QBX sources - def __call__(self, queue, lpot_source, targets_list=()): - """ - Return a :class:`TreeWithQBXMetadata` built from the given layer - potential source. + .. attribute:: box_to_qbx_source_starts - :arg queue: An instance of :class:`pyopencl.CommandQueue` + ``box_id_t [nboxes + 1]`` - :arg lpot_source: An instance of - :class:`pytential.qbx.NewQBXLayerPotentialSource`. + .. attribute:: box_to_qbx_source_lists - :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( + ``particle_id_t [*]`` + + .. rubric:: Box to QBX centers + + .. attribute:: box_to_qbx_center_starts + + ``box_id_t [nboxes + 1]`` + + .. attribute:: box_to_qbx_center_lists + + ``particle_id_t [*]`` + + .. rubric:: 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 + .. attribute:: qbx_user_target_slice + """ + pass + + +MAX_REFINE_WEIGHT = 64 + + +def build_tree_with_qbx_metadata( + queue, tree_builder, lpot_source, targets_list=(), + 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` + + :arg lpot_source: An instance of + :class:`pytential.qbx.NewQBXLayerPotentialSource`. + + :arg targets_list: A list of :class:`pytential.target.TargetBase` + + :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: + # - 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_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_base_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_base_fine_discr + ntargets = sum(tgt.nnodes for tgt in targets_list) + + # Slices + qbx_user_source_slice = slice(0, nsources) + + 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) + + 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. + 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, -center_slice_start), + ("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_base_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 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) + 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_base_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 d7e514510d320602496a8b13a562a9da99733eba..4013b8f31081d6e0d4deda9a846fd4bf0e0af552 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 @@ -53,134 +53,159 @@ RNG_SEED = 10 FAR_TARGET_DIST_FROM_SOURCE = 10 -# {{{ utilities for iterating over panels +# {{{ source refinement checker 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 = QuadratureSimplexGroupFactory(4 * order) 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) + 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, **refiner_extra_kwargs) 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 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]) + def check_disk_undisturbed_by_sources(centers_panel, sources_panel): + h = panel_sizes[centers_panel.element_nr] - # Check wavenumber to panel size ratio. - assert panel_sizes[panel.element_nr] * helmholtz_k <= 5 - - 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) 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()) - # Criterion 1: + # Criterion: # 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) - # Criterion 2: - # A center cannot be closer to another panel than that panel's - # centers - unless the panels are adjacent, to allow for refinement. + def check_sufficient_quadrature_resolution(centers_panel, sources_panel): + h = fine_panel_sizes[sources_panel.element_nr] - if panel_2.element_nr in panel_1.neighbors: - return + 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) - assert dist >= h_2 / 2, (dist, h_2, panel_1.element_nr, panel_2.element_nr) + nodes = fine_discr_nodes[:, sources_panel.discr_slice] - for panel_1 in iter_elements(lpot_source.density_discr): - check_panel(panel_1) + # =distance(centers of panel 1, panel 2) + dist = ( + la.norm(( + all_centers[..., np.newaxis] - + nodes[:, np.newaxis, ...]).T, + axis=-1) + .min()) + + # 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) + + 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 i, panel_1 in enumerate(iter_elements(lpot_source.density_discr)): 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) + 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): + 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), @@ -204,13 +229,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) diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 1c704e5226058b8b19b454413ebb8d809d5b6057..870cc6f0438502e767dc23426d839243619a3023 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