diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 8a366ae4ac945acd6cb26a8bca852003ead97719..60f439ef09a3a37b3b997f4f2afa933ec2f04eb1 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -331,7 +331,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): @memoize_method def with_refinement(self, target_order=None, kernel_length_scale=None, - maxiter=None, visualize=False): + maxiter=None, visualize=False, _expansion_disturbance_tolerance=None): """ :returns: a tuple ``(lpot_src, cnx)``, where ``lpot_src`` is a :class:`QBXLayerPotentialSource` and ``cnx`` is a @@ -352,7 +352,8 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): self.refiner_code_container.get_wrangler(queue), InterpolatoryQuadratureSimplexGroupFactory(target_order), kernel_length_scale=kernel_length_scale, - maxiter=maxiter, visualize=visualize) + maxiter=maxiter, visualize=visualize, + expansion_disturbance_tolerance=_expansion_disturbance_tolerance) return lpot, connection diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index e716fe46e359dddb4ed038bd51de20220eba80a6..ce183e582cbb45e018e3e0d7ef90c6313b2dd000 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -94,7 +94,8 @@ 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 *center_danger_zone_radii_by_panel, + coord_t *center_danger_zone_radii, + coord_t expansion_disturbance_tolerance, int npanels, /* output */ @@ -107,11 +108,11 @@ 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); + particle_id_t icenter = i; - ${load_particle("INDEX_FOR_CENTER_PARTICLE(i)", ball_center)} - ${ball_radius} = center_danger_zone_radii_by_panel[my_panel]; + ${load_particle("INDEX_FOR_CENTER_PARTICLE(icenter)", ball_center)} + ${ball_radius} = (1-expansion_disturbance_tolerance) + * center_danger_zone_radii[icenter]; """, leaf_found_op=QBX_TREE_MAKO_DEFS + r""" /* Check that each source in the leaf box is sufficiently far from the @@ -122,24 +123,24 @@ EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( ++source_idx) { particle_id_t source = box_to_source_lists[source_idx]; - particle_id_t panel = bsearch( + particle_id_t source_panel = bsearch( panel_to_source_starts, npanels + 1, source); - if (my_panel == panel) - { - continue; - } + /* Find the panel associated with this center. */ + particle_id_t center_panel = bsearch(panel_to_center_starts, npanels + 1, + icenter); coord_vec_t source_coords; ${load_particle("INDEX_FOR_SOURCE_PARTICLE(source)", "source_coords")} bool is_close = ( distance(${ball_center}, source_coords) - <= center_danger_zone_radii_by_panel[my_panel]); + <= (1-expansion_disturbance_tolerance) + * center_danger_zone_radii[icenter]); if (is_close) { - panel_refine_flags[my_panel] = 1; + panel_refine_flags[center_panel] = 1; *found_panel_to_refine = 1; break; } @@ -286,7 +287,9 @@ class RefinerWrangler(TreeWranglerBase): # {{{ check subroutines for conditions 1-3 def check_expansion_disks_undisturbed_by_sources(self, - lpot_source, tree, peer_lists, refine_flags, + lpot_source, tree, peer_lists, + expansion_disturbance_tolerance, + refine_flags, debug, wait_for=None): # Avoid generating too many kernels. @@ -301,7 +304,8 @@ class RefinerWrangler(TreeWranglerBase): tree.particle_id_dtype, max_levels) - logger.info("refiner: checking center is closest to orig panel") + logger.info("refiner: checking that expansion disk is " + "undisturbed by sources") if debug: npanels_to_refine_prev = cl.array.sum(refine_flags).get() @@ -310,9 +314,7 @@ class RefinerWrangler(TreeWranglerBase): found_panel_to_refine.finish() unwrap_args = AreaQueryElementwiseTemplate.unwrap_args - # FIXME: This shouldn't be by panel - center_danger_zone_radii_by_panel = \ - lpot_source._expansion_radii("npanels") + center_danger_zone_radii = lpot_source._expansion_radii("ncenters") evt = knl( *unwrap_args( @@ -324,7 +326,8 @@ class RefinerWrangler(TreeWranglerBase): tree.qbx_user_source_slice.start, tree.qbx_user_center_slice.start, tree.sorted_target_ids, - center_danger_zone_radii_by_panel, + center_danger_zone_radii, + expansion_disturbance_tolerance, tree.nqbxpanels, refine_flags, found_panel_to_refine, @@ -480,7 +483,7 @@ def make_empty_refine_flags(queue, lpot_source, use_stage2_discr=False): def refine_for_global_qbx(lpot_source, wrangler, group_factory, kernel_length_scale=None, refine_flags=None, debug=None, maxiter=None, - visualize=None): + visualize=None, expansion_disturbance_tolerance=None): """ Entry point for calling the refiner. @@ -515,6 +518,9 @@ def refine_for_global_qbx(lpot_source, wrangler, # FIXME: Set debug=False by default once everything works. debug = True + if expansion_disturbance_tolerance is None: + expansion_disturbance_tolerance = 0.025 + # TODO: Stop doing redundant checks by avoiding panels which no longer need # refinement. @@ -535,9 +541,6 @@ def refine_for_global_qbx(lpot_source, wrangler, # {{{ first stage refinement - must_refine = True - niter = 0 - def visualize_refinement(niter, stage, flags): if not visualize: return @@ -566,16 +569,38 @@ def refine_for_global_qbx(lpot_source, wrangler, vis.write_vtk_file("refinement-%03d-%s.vtu" % (niter, stage), vis_data) - while must_refine: - must_refine = False + def warn_max_iterations(): + from warnings import warn + warn( + "QBX layer potential source refiner did not terminate " + "after %d iterations (the maximum). " + "You may pass 'visualize=True' to with_refinement() " + "to see what area of the geometry is causing trouble. " + "If the issue is disturbance of expansion disks, you may " + "pass a slightly increased value (say 1e-2) for " + "_expansion_disturbance_tolerance in with_refinement(). " + "As a last resort, " + "you may use Python's warning filtering mechanism to " + "not treat this warning as an error. " + "The criteria triggering refinement in each iteration " + "were: %s. " % ( + len(violated_criteria), + ", ".join( + "%d: %s" % (i+1, vc_text) + for i, vc_text in enumerate(violated_criteria))), + RefinerNotConvergedWarning) + + violated_criteria = [] + iter_violated_criteria = ["start"] + + niter = 0 + + while iter_violated_criteria: + iter_violated_criteria = [] niter += 1 if niter > maxiter: - from warnings import warn - warn( - "Max iteration count reached in QBX layer potential source" - " refiner.", - RefinerNotConvergedWarning) + warn_max_iterations() break # Build tree and auxiliary data. @@ -587,9 +612,11 @@ def refine_for_global_qbx(lpot_source, wrangler, # Check condition 1. has_disturbed_expansions = \ wrangler.check_expansion_disks_undisturbed_by_sources( - lpot_source, tree, peer_lists, refine_flags, debug) - must_refine = must_refine or has_disturbed_expansions + lpot_source, tree, peer_lists, + expansion_disturbance_tolerance, + refine_flags, debug) if has_disturbed_expansions: + iter_violated_criteria.append("disturbed expansions") visualize_refinement(niter, "disturbed-expansions", refine_flags) # Check condition 3. @@ -598,12 +625,14 @@ def refine_for_global_qbx(lpot_source, wrangler, violates_kernel_length_scale = \ wrangler.check_kernel_length_scale_to_panel_size_ratio( lpot_source, kernel_length_scale, refine_flags, debug) - must_refine = must_refine or violates_kernel_length_scale if violates_kernel_length_scale: + iter_violated_criteria.append("kernel length scale") visualize_refinement(niter, "kernel-length-scale", refine_flags) - if must_refine: + if iter_violated_criteria: + violated_criteria.append(" and ".join(iter_violated_criteria)) + conn = wrangler.refine( lpot_source.density_discr, refiner, refine_flags, group_factory, debug) @@ -618,22 +647,18 @@ def refine_for_global_qbx(lpot_source, wrangler, # {{{ second stage refinement - must_refine = True + iter_violated_criteria = ["start"] niter = 0 fine_connections = [] stage2_density_discr = lpot_source.density_discr - while must_refine: - must_refine = False + while iter_violated_criteria: + iter_violated_criteria = [] niter += 1 if niter > maxiter: - from warnings import warn - warn( - "Max iteration count reached in QBX layer potential source" - " refiner.", - RefinerNotConvergedWarning) + warn_max_iterations() break # Build tree and auxiliary data. @@ -646,11 +671,13 @@ def refine_for_global_qbx(lpot_source, wrangler, has_insufficient_quad_res = \ wrangler.check_sufficient_source_quadrature_resolution( lpot_source, tree, peer_lists, refine_flags, debug) - must_refine = must_refine or has_insufficient_quad_res if has_insufficient_quad_res: + iter_violated_criteria.append("insufficient quadrature resolution") visualize_refinement(niter, "quad-resolution", refine_flags) - if must_refine: + if iter_violated_criteria: + violated_criteria.append(" and ".join(iter_violated_criteria)) + conn = wrangler.refine( stage2_density_discr, refiner, refine_flags, group_factory, debug) diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index e3a0575b8fea84a468ac580c5fb0b51d38ede377..fe45746e79e6fe1e03738f278601946cf4260dca 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -97,7 +97,10 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): lpot_source = QBXLayerPotentialSource(discr, order) del discr - refiner_extra_kwargs = {} + expansion_disturbance_tolerance = 0.025 + refiner_extra_kwargs = { + "expansion_disturbance_tolerance": expansion_disturbance_tolerance, + } if helmholtz_k is not None: refiner_extra_kwargs["kernel_length_scale"] = 5/helmholtz_k @@ -144,7 +147,7 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None): # panel. rad = expansion_radii[centers_panel.element_nr] - assert dist >= rad, \ + assert dist >= rad * (1-expansion_disturbance_tolerance), \ (dist, rad, centers_panel.element_nr, sources_panel.element_nr) def check_sufficient_quadrature_resolution(centers_panel, sources_panel):