From 7c9eb77e84e9a8cfaa65975fa0d6e9eaeb6fc4ed Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 21 Sep 2017 20:53:11 -0500 Subject: [PATCH 1/6] Refiner: Add expansion disturbance tolerance, use per-center radii --- pytential/qbx/__init__.py | 5 ++-- pytential/qbx/refinement.py | 49 +++++++++++++++++++++++-------------- 2 files changed, 33 insertions(+), 21 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 8a366ae4..57f67b2f 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=0): """ :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 e716fe46..c7d21ec4 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,27 @@ 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) - { + /* Find the panel associated with this center. */ + particle_id_t center_panel = bsearch(panel_to_center_starts, npanels + 1, + icenter); + + if (center_panel == source_panel) continue; - } 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 +290,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 +307,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 +317,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 +329,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 +486,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 +521,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 + # TODO: Stop doing redundant checks by avoiding panels which no longer need # refinement. @@ -587,7 +596,9 @@ 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) + lpot_source, tree, peer_lists, + expansion_disturbance_tolerance, + refine_flags, debug) must_refine = must_refine or has_disturbed_expansions if has_disturbed_expansions: visualize_refinement(niter, "disturbed-expansions", refine_flags) -- GitLab From 9255d9ccd3563325f91a626ec1878edf29e79296 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 21 Sep 2017 20:57:33 -0500 Subject: [PATCH 2/6] Refiner: generate more useful error messages --- pytential/qbx/refinement.py | 63 +++++++++++++++++++++++-------------- 1 file changed, 40 insertions(+), 23 deletions(-) diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index c7d21ec4..5dcd4de7 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -544,9 +544,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 @@ -575,16 +572,36 @@ 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_not_converging(): + 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 tiggering refinement in each iteration " + "were: %s. " % ( + len(violated_criteria), + ", ".join(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_not_converging() break # Build tree and auxiliary data. @@ -599,8 +616,8 @@ def refine_for_global_qbx(lpot_source, wrangler, lpot_source, tree, peer_lists, expansion_disturbance_tolerance, refine_flags, debug) - must_refine = must_refine or has_disturbed_expansions if has_disturbed_expansions: + iter_violated_criteria.append("disturbed expansions") visualize_refinement(niter, "disturbed-expansions", refine_flags) # Check condition 3. @@ -609,12 +626,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) @@ -629,22 +648,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_not_converging() break # Build tree and auxiliary data. @@ -657,11 +672,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) -- GitLab From 84ac7a259f0e7aef2a19f7025015d9bdc76a2066 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 22 Sep 2017 10:24:12 -0500 Subject: [PATCH 3/6] Refiner expansion disturbance check: Non-zero tolerance by default, stop giving the self panel a free pass on that criterion --- pytential/qbx/__init__.py | 2 +- pytential/qbx/refinement.py | 5 +---- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 57f67b2f..3a737a2b 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, _expansion_disturbance_tolerance=0): + maxiter=None, visualize=False, _expansion_disturbance_tolerance=0.025): """ :returns: a tuple ``(lpot_src, cnx)``, where ``lpot_src`` is a :class:`QBXLayerPotentialSource` and ``cnx`` is a diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 5dcd4de7..99b3a36a 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -130,9 +130,6 @@ EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate( particle_id_t center_panel = bsearch(panel_to_center_starts, npanels + 1, icenter); - if (center_panel == source_panel) - continue; - coord_vec_t source_coords; ${load_particle("INDEX_FOR_SOURCE_PARTICLE(source)", "source_coords")} @@ -522,7 +519,7 @@ def refine_for_global_qbx(lpot_source, wrangler, debug = True if expansion_disturbance_tolerance is None: - expansion_disturbance_tolerance = 0 + raise ValueError("must specify expansion_disturbance_tolerance") # TODO: Stop doing redundant checks by avoiding panels which no longer need # refinement. -- GitLab From e75d406c2a5c7ed05c27cefe6f644f67e8eec0ab Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 22 Sep 2017 10:27:19 -0500 Subject: [PATCH 4/6] Improve refiner max-iter message --- pytential/qbx/refinement.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 99b3a36a..9e1a57c1 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -569,7 +569,7 @@ def refine_for_global_qbx(lpot_source, wrangler, vis.write_vtk_file("refinement-%03d-%s.vtu" % (niter, stage), vis_data) - def warn_not_converging(): + def warn_max_iterations(): from warnings import warn warn( "QBX layer potential source refiner did not terminate " @@ -582,10 +582,12 @@ def refine_for_global_qbx(lpot_source, wrangler, "As a last resort, " "you may use Python's warning filtering mechanism to " "not treat this warning as an error. " - "The criteria tiggering refinement in each iteration " + "The criteria triggering refinement in each iteration " "were: %s. " % ( len(violated_criteria), - ", ".join(violated_criteria)), + ", ".join( + "%d: %s" % (i+1, vc_text) + for i, vc_text in enumerate(violated_criteria))), RefinerNotConvergedWarning) violated_criteria = [] @@ -598,7 +600,7 @@ def refine_for_global_qbx(lpot_source, wrangler, niter += 1 if niter > maxiter: - warn_not_converging() + warn_max_iterations() break # Build tree and auxiliary data. @@ -656,7 +658,7 @@ def refine_for_global_qbx(lpot_source, wrangler, niter += 1 if niter > maxiter: - warn_not_converging() + warn_max_iterations() break # Build tree and auxiliary data. -- GitLab From 3e6f21b77e78fbc81471b81f660f9af842c0e0d5 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 22 Sep 2017 13:13:17 -0500 Subject: [PATCH 5/6] Refiner: push expansion disturbance tolerance default down to refine_for_global_qbx entrypoint --- pytential/qbx/__init__.py | 2 +- pytential/qbx/refinement.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 3a737a2b..60f439ef 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, _expansion_disturbance_tolerance=0.025): + 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 diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 9e1a57c1..ce183e58 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -519,7 +519,7 @@ def refine_for_global_qbx(lpot_source, wrangler, debug = True if expansion_disturbance_tolerance is None: - raise ValueError("must specify expansion_disturbance_tolerance") + expansion_disturbance_tolerance = 0.025 # TODO: Stop doing redundant checks by avoiding panels which no longer need # refinement. -- GitLab From e4c1b9e4c5fa0e949981b3b6a49273e7934dd732 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 23 Sep 2017 10:28:21 -0500 Subject: [PATCH 6/6] Tweak global QBX refiner test to take into account expansion_disturbance_tolerance --- test/test_global_qbx.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index e3a0575b..fe45746e 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): -- GitLab