From bb129448c4e40dea67554016c516f9dda4550a53 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 9 Jun 2017 04:29:34 -0500 Subject: [PATCH 1/9] Make performance model work again and make it possible for the user to tweak the stick out. --- pytential/qbx/__init__.py | 12 +++++++++++- pytential/qbx/fmm.py | 11 ++++++----- pytential/qbx/geometry.py | 4 ++-- 3 files changed, 19 insertions(+), 8 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 5292f7bb..3d9a52e8 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -78,6 +78,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): debug=True, refined_for_global_qbx=False, expansion_disks_in_tree_have_extent=False, + _expansion_disk_stick_out_factor=0, performance_data_file=None, fmm_backend="sumpy"): """ @@ -130,6 +131,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): self.refined_for_global_qbx = refined_for_global_qbx self.expansion_disks_in_tree_have_extent = \ expansion_disks_in_tree_have_extent + self._expansion_disk_stick_out_factor = _expansion_disk_stick_out_factor self.performance_data_file = performance_data_file def copy( @@ -140,6 +142,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_level_to_order=None, target_stick_out_factor=None, base_resampler=None, + _expansion_disk_stick_out_factor=None, + performance_data_file=None, debug=None, refined_for_global_qbx=None, @@ -166,7 +170,13 @@ class QBXLayerPotentialSource(LayerPotentialSource): else self.refined_for_global_qbx), expansion_disks_in_tree_have_extent=( self.expansion_disks_in_tree_have_extent), - performance_data_file=self.performance_data_file, + _expansion_disk_stick_out_factor=( + # 0 is a valid value here + _expansion_disk_stick_out_factor + if _expansion_disk_stick_out_factor is not None + else self._expansion_disk_stick_out_factor), + performance_data_file=( + performance_data_file or self.performance_data_file), fmm_backend=self.fmm_backend) # }}} diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index c7d6c029..49026ef5 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -600,9 +600,9 @@ def write_performance_model(outf, geo_data): for itgt_box, tgt_ibox in enumerate(traversal.target_boxes): ntargets = box_target_counts_nonchild[tgt_ibox] - start, end = traversal.sep_smaller_starts[itgt_box:itgt_box+2] - - nmp_eval += ntargets * (end-start) + for sep_smaller_list in traversal.sep_smaller_by_level: + start, end = sep_smaller_list.starts[itgt_box:itgt_box+2] + nmp_eval += ntargets * (end-start) outf.write("mp_eval = {cost}\n" .format(cost=nmp_eval * p_fmm)) @@ -680,8 +680,9 @@ def write_performance_model(outf, geo_data): for itgt_center, tgt_icenter in enumerate(global_qbx_centers): itgt_box = qbx_center_to_target_box[tgt_icenter] - start, end = traversal.sep_smaller_starts[itgt_box:itgt_box+2] - nqbx_m2l += end - start + for sep_smaller_list in traversal.sep_smaller_by_level: + start, end = sep_smaller_list.starts[tgt_ibox:tgt_ibox+2] + nqbx_m2l += ntargets * (end-start) outf.write("qbx_m2l = {cost}\n" .format(cost=nqbx_m2l * p_fmm * p_qbx)) diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index f2669274..2c681b89 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -517,10 +517,10 @@ class QBXFMMGeometryData(object): particles=lpot_src.fine_density_discr.nodes(), targets=target_info.targets, target_radii=target_radii, - max_leaf_refine_weight=256, + max_leaf_refine_weight=32, refine_weights=refine_weights, debug=self.debug, - stick_out_factor=0, + stick_out_factor=lpot_src._expansion_disk_stick_out_factor, kind="adaptive") if self.debug: -- GitLab From 497e73cb9de4576d34ce0447ff68682f189ab1ea Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 9 Jun 2017 18:04:46 -0400 Subject: [PATCH 2/9] Revamp performance model to prevent identifier leakage [ci skip] --- pytential/qbx/fmm.py | 153 +++++++++++++++++++++++++------------------ 1 file changed, 88 insertions(+), 65 deletions(-) diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index 49026ef5..c3ff6433 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -565,64 +565,75 @@ def write_performance_model(outf, geo_data): # {{{ direct evaluation from neighbor source boxes ("list 1") - npart_direct = 0 - for itgt_box, tgt_ibox in enumerate(traversal.target_boxes): - ntargets = box_target_counts_nonchild[tgt_ibox] + def process_list1(): + npart_direct = 0 + for itgt_box, tgt_ibox in enumerate(traversal.target_boxes): + ntargets = box_target_counts_nonchild[tgt_ibox] - start, end = traversal.neighbor_source_boxes_starts[itgt_box:itgt_box+2] - for src_ibox in traversal.neighbor_source_boxes_lists[start:end]: - nsources = tree.box_source_counts_nonchild[src_ibox] + start, end = traversal.neighbor_source_boxes_starts[itgt_box:itgt_box+2] + for src_ibox in traversal.neighbor_source_boxes_lists[start:end]: + nsources = tree.box_source_counts_nonchild[src_ibox] - npart_direct += ntargets * nsources + npart_direct += ntargets * nsources - outf.write("part_direct = {cost}\n" - .format(cost=npart_direct)) + outf.write("part_direct = {cost}\n" + .format(cost=npart_direct)) + + process_list1() # }}} # {{{ translate separated siblings' ("list 2") mpoles to local - nm2l = 0 - for itgt_box, tgt_ibox in enumerate(traversal.target_or_target_parent_boxes): - start, end = traversal.sep_siblings_starts[itgt_box:itgt_box+2] + def process_list2(): + nm2l = 0 + for itgt_box, tgt_ibox in enumerate(traversal.target_or_target_parent_boxes): + start, end = traversal.sep_siblings_starts[itgt_box:itgt_box+2] + + nm2l += (end-start) - nm2l += (end-start) + outf.write("m2l = {cost}\n" + .format(cost=nm2l * p_fmm**2)) - outf.write("m2l = {cost}\n" - .format(cost=nm2l * p_fmm**2)) + process_list2() # }}} # {{{ evaluate sep. smaller mpoles ("list 3") at particles - nmp_eval = 0 + def process_list3(): + nmp_eval = 0 + for itgt_box, tgt_ibox in enumerate(traversal.target_boxes): + ntargets = box_target_counts_nonchild[tgt_ibox] - for itgt_box, tgt_ibox in enumerate(traversal.target_boxes): - ntargets = box_target_counts_nonchild[tgt_ibox] + for sep_smaller_list in traversal.sep_smaller_by_level: + start, end = sep_smaller_list.starts[itgt_box:itgt_box+2] + nmp_eval += ntargets * (end-start) - for sep_smaller_list in traversal.sep_smaller_by_level: - start, end = sep_smaller_list.starts[itgt_box:itgt_box+2] - nmp_eval += ntargets * (end-start) + outf.write("mp_eval = {cost}\n" + .format(cost=nmp_eval * p_fmm)) - outf.write("mp_eval = {cost}\n" - .format(cost=nmp_eval * p_fmm)) + process_list3() # }}} # {{{ form locals for separated bigger mpoles ("list 4") - nform_local = 0 + def process_list4(): + nform_local = 0 + + for itgt_box, tgt_ibox in enumerate(traversal.target_or_target_parent_boxes): + start, end = traversal.sep_bigger_starts[itgt_box:itgt_box+2] - for itgt_box, tgt_ibox in enumerate(traversal.target_or_target_parent_boxes): - start, end = traversal.sep_bigger_starts[itgt_box:itgt_box+2] + for src_ibox in traversal.sep_bigger_lists[start:end]: + nsources = tree.box_source_counts_nonchild[src_ibox] - for src_ibox in traversal.sep_bigger_lists[start:end]: - nsources = tree.box_source_counts_nonchild[src_ibox] + nform_local += nsources - nform_local += nsources + outf.write("form_local = {cost}\n" + .format(cost=nform_local * p_fmm)) - outf.write("form_local = {cost}\n" - .format(cost=nform_local * p_fmm)) + process_list4() # }}} @@ -642,72 +653,84 @@ def write_performance_model(outf, geo_data): # {{{ form global qbx locals - nqbxl_direct = 0 - global_qbx_centers = geo_data.global_qbx_centers() qbx_center_to_target_box = geo_data.qbx_center_to_target_box() center_to_targets_starts = geo_data.center_to_tree_targets().starts - ncenters = geo_data.ncenters + with cl.CommandQueue(geo_data.cl_context) as queue: + global_qbx_centers = global_qbx_centers.get( + queue=queue) + qbx_center_to_target_box = qbx_center_to_target_box.get( + queue=queue) + center_to_targets_starts = center_to_targets_starts.get( + queue=queue) - outf.write("ncenters = {cost}\n" - .format(cost=ncenters)) + def process_form_qbxl(): + nqbxl_direct = 0 - with cl.CommandQueue(geo_data.cl_context) as queue: - global_qbx_centers = global_qbx_centers.get(queue=queue) - qbx_center_to_target_box = qbx_center_to_target_box.get(queue=queue) - center_to_targets_starts = center_to_targets_starts.get(queue=queue) + ncenters = geo_data.ncenters - for itgt_center, tgt_icenter in enumerate(global_qbx_centers): - itgt_box = qbx_center_to_target_box[tgt_icenter] - tgt_ibox = traversal.target_or_target_parent_boxes[itgt_box] + outf.write("ncenters = {cost}\n" + .format(cost=ncenters)) - start, end = traversal.neighbor_source_boxes_starts[itgt_box:itgt_box+2] - for src_ibox in traversal.neighbor_source_boxes_lists[start:end]: - nsources = tree.box_source_counts_nonchild[src_ibox] + for itgt_center, tgt_icenter in enumerate(global_qbx_centers): + itgt_box = qbx_center_to_target_box[tgt_icenter] - nqbxl_direct += nsources + start, end = traversal.neighbor_source_boxes_starts[itgt_box:itgt_box+2] + for src_ibox in traversal.neighbor_source_boxes_lists[start:end]: + nsources = tree.box_source_counts_nonchild[src_ibox] - outf.write("qbxl_direct = {cost}\n" - .format(cost=nqbxl_direct * p_qbx)) + nqbxl_direct += nsources + + outf.write("qbxl_direct = {cost}\n" + .format(cost=nqbxl_direct * p_qbx)) # }}} # {{{ translate from list 3 multipoles to qbx local expansions - nqbx_m2l = 0 + def process_m2qbxl(): + nqbx_m2l = 0 + + for isrc_level, ssn in enumerate(traversal.sep_smaller_by_level): + for itgt_center, tgt_icenter in enumerate(global_qbx_centers): + icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter] - for itgt_center, tgt_icenter in enumerate(global_qbx_centers): - itgt_box = qbx_center_to_target_box[tgt_icenter] + start, stop = ( + ssn.starts[icontaining_tgt_box], + ssn.starts[icontaining_tgt_box+1]) - for sep_smaller_list in traversal.sep_smaller_by_level: - start, end = sep_smaller_list.starts[tgt_ibox:tgt_ibox+2] - nqbx_m2l += ntargets * (end-start) + nqbx_m2l += stop-start - outf.write("qbx_m2l = {cost}\n" - .format(cost=nqbx_m2l * p_fmm * p_qbx)) + outf.write("qbx_m2l = {cost}\n" + .format(cost=nqbx_m2l * p_fmm * p_qbx)) + + process_m2qbxl() # }}} # {{{ translate from box local expansions to qbx local expansions outf.write("qbx_l2l = {cost}\n" - .format(cost=ncenters * p_fmm * p_qbx)) + .format(cost=geo_data.ncenters * p_fmm * p_qbx)) # }}} # {{{ evaluate qbx local expansions - nqbx_eval = 0 + def process_eval_qbxl(): + nqbx_eval = 0 + + for iglobal_center in range(geo_data.ncenters): + src_icenter = global_qbx_centers[iglobal_center] - for iglobal_center in range(ncenters): - src_icenter = global_qbx_centers[iglobal_center] + start, end = center_to_targets_starts[src_icenter:src_icenter+2] + nqbx_eval += end-start - start, end = center_to_targets_starts[src_icenter:src_icenter+2] - nqbx_eval += end-start + outf.write("qbx_eval = {cost}\n" + .format(cost=nqbx_eval * p_qbx)) - outf.write("qbx_eval = {cost}\n" - .format(cost=nqbx_eval * p_qbx)) + process_eval_qbxl() # }}} -- GitLab From a934c93b6c44aae60285f7cb1f1d86ec5405e903 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 10 Jun 2017 12:00:04 -0400 Subject: [PATCH 3/9] Perf estimate fix: Add missing process_form_qbxl [ci skip] --- pytential/qbx/fmm.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index c3ff6433..6afa047b 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -685,6 +685,8 @@ def write_performance_model(outf, geo_data): outf.write("qbxl_direct = {cost}\n" .format(cost=nqbxl_direct * p_qbx)) + process_form_qbxl() + # }}} # {{{ translate from list 3 multipoles to qbx local expansions -- GitLab From c71db26e239bae901a73757a5fb20fe9361eb7bc Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 15 Jun 2017 16:09:08 -0500 Subject: [PATCH 4/9] Rework argument passing for the QBX LP source, rename target stick-out to target assoc tol (Closes #48 on gitlab) --- examples/cahn-hilliard.py | 2 +- examples/helmholtz-dirichlet.py | 2 +- examples/layerpot.py | 2 +- examples/maxwell.py | 4 +- examples/scaling-study.py | 4 +- examples/stokes-2d-interior.py | 33 ++++----- pytential/qbx/__init__.py | 115 ++++++++++++++++++++++++-------- pytential/qbx/geometry.py | 11 +-- pytential/qbx/refinement.py | 2 +- pytential/qbx/target_assoc.py | 22 +++--- test/test_layer_pot.py | 8 +-- test/test_stokes.py | 5 +- 12 files changed, 138 insertions(+), 72 deletions(-) diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index f7cd6559..8393f744 100644 --- a/examples/cahn-hilliard.py +++ b/examples/cahn-hilliard.py @@ -91,7 +91,7 @@ def main(): targets = cl.array.to_device(queue, fplot.points) - qbx_stick_out = qbx.copy(target_stick_out_factor=0.05) + qbx_stick_out = qbx.copy(target_association_tolerance=0.05) indicator_qbx = qbx_stick_out.copy(qbx_order=2) diff --git a/examples/helmholtz-dirichlet.py b/examples/helmholtz-dirichlet.py index c1f4719c..b591c5be 100644 --- a/examples/helmholtz-dirichlet.py +++ b/examples/helmholtz-dirichlet.py @@ -143,7 +143,7 @@ def main(): u_incoming = u_incoming_func(targets) - qbx_stick_out = qbx.copy(target_stick_out_factor=0.05) + qbx_stick_out = qbx.copy(target_association_tolerance=0.05) indicator_qbx = qbx_stick_out.copy( fmm_level_to_order=lambda lev: 7, qbx_order=2) diff --git a/examples/layerpot.py b/examples/layerpot.py index 5e6bbcf4..111265cb 100644 --- a/examples/layerpot.py +++ b/examples/layerpot.py @@ -54,7 +54,7 @@ pre_density_discr = Discretization( qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4*target_order, qbx_order, fmm_order=qbx_order+3, - target_stick_out_factor=0.005).with_refinement() + target_association_tolerance=0.005).with_refinement() density_discr = qbx.density_discr diff --git a/examples/maxwell.py b/examples/maxwell.py index b3f016a8..ed8c6106 100644 --- a/examples/maxwell.py +++ b/examples/maxwell.py @@ -99,12 +99,12 @@ def main(): fplot = FieldPlotter(bbox_center, extent=2*bbox_size, npoints=(150, 150, 1)) - qbx_stick_out = qbx.copy(target_stick_out_factor=0.1) + qbx_tgt_tol = qbx.copy(target_association_tolerance=0.1) from pytential.target import PointsTarget from pytential.qbx import QBXTargetAssociationFailedException try: fld_in_vol = bind( - (qbx_stick_out, PointsTarget(fplot.points)), + (qbx_tgt_tol, PointsTarget(fplot.points)), sym.make_obj_array([ sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), qbx_forced_limit=None), diff --git a/examples/scaling-study.py b/examples/scaling-study.py index cf8cf166..183fc915 100644 --- a/examples/scaling-study.py +++ b/examples/scaling-study.py @@ -142,9 +142,9 @@ def timing_run(nx, ny): targets = cl.array.to_device(queue, fplot.points) - qbx_stick_out = qbx.copy(target_stick_out_factor=0.05) + qbx_tgt_tol = qbx.copy(target_association_tolerance=0.05) - indicator_qbx = qbx_stick_out.copy( + indicator_qbx = qbx_tgt_tol.copy( fmm_level_to_order=lambda lev: 7, qbx_order=2) ones_density = density_discr.zeros(queue) diff --git a/examples/stokes-2d-interior.py b/examples/stokes-2d-interior.py index 4f6e7bee..977ae832 100644 --- a/examples/stokes-2d-interior.py +++ b/examples/stokes-2d-interior.py @@ -23,7 +23,7 @@ fmm_order = 7 mu = 3 # Test solution type -- either 'fundamental' or 'couette' (default is couette) -soln_type = 'couette' +soln_type = 'couette' # }}} @@ -53,19 +53,20 @@ def main(nelements): InterpolatoryQuadratureSimplexGroupFactory(target_order)) from pytential.qbx import QBXLayerPotentialSource - stick_out = 0.05 + target_association_tolerance = 0.05 qbx, _ = QBXLayerPotentialSource( coarse_density_discr, fine_order=ovsmp_target_order, qbx_order=qbx_order, - fmm_order=fmm_order, target_stick_out_factor=stick_out + fmm_order=fmm_order, + target_association_tolerance=target_association_tolerance, ).with_refinement() density_discr = qbx.density_discr nodes = density_discr.nodes().with_queue(queue) - # Get normal vectors for the density discretization -- used in integration with stresslet + # Get normal vectors for the density discretization -- used in integration with stresslet mv_normal = bind(density_discr, sym.normal(2))(queue) normal = mv_normal.as_vector(np.object) - + # {{{ describe bvp @@ -88,7 +89,7 @@ def main(nelements): # Create stresslet object stresslet_obj = StressletWrapper(dim=2) - # Describe boundary operator + # Describe boundary operator bdry_op_sym = loc_sign * 0.5 * sigma_sym + sqrt_w * stresslet_obj.apply(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit='avg') # Bind to the qbx discretization @@ -106,14 +107,14 @@ def main(nelements): xcomp = (-cl.clmath.log(r) + (x - loc[0])**2/r**2) * scaling ycomp = ((x - loc[0])*(y - loc[1])/r**2) * scaling return [ xcomp, ycomp ] - + def couette_soln(x, y, dp, h): scaling = 1./(2*mu) xcomp = scaling * dp * ((y+(h/2.))**2 - h * (y+(h/2.))) ycomp = scaling * 0*y return [xcomp, ycomp] - - + + if soln_type == 'fundamental': pt_loc = np.array([2.0, 0.0]) @@ -122,8 +123,8 @@ def main(nelements): dp = -10. h = 2.5 bc = couette_soln(nodes[0], nodes[1], dp, h) - - # Get rhs vector + + # Get rhs vector bvp_rhs = bind(qbx, sqrt_w*sym.make_sym_vector("bc",dim))(queue, bc=bc) from pytential.solve import gmres @@ -137,7 +138,7 @@ def main(nelements): # {{{ postprocess/visualize sigma = gmres_result.solution - + # Describe representation of solution for evaluation in domain representation_sym = stresslet_obj.apply(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit=-2) @@ -163,7 +164,7 @@ def main(nelements): return np.array([ row[mask] for row in test_points]) - + def stride_hack(arr): from numpy.lib.stride_tricks import as_strided return np.array(as_strided(arr, strides=(8 * len(arr[0]), 8))) @@ -180,13 +181,13 @@ def main(nelements): if soln_type == 'fundamental': - exact_soln = fund_soln(eval_points_dev[0], eval_points_dev[1], pt_loc) + exact_soln = fund_soln(eval_points_dev[0], eval_points_dev[1], pt_loc) else: exact_soln = couette_soln(eval_points_dev[0], eval_points_dev[1], dp, h) err = vel - get_obj_array(exact_soln) print("@@@@@@@@") - + print("L2 error estimate: ", np.sqrt((2./(nsamp-1))**2*np.sum(err[0]*err[0]) + (2./(nsamp-1))**2*np.sum(err[1]*err[1]))) max_error_loc = [abs(err[0]).argmax(), abs(err[1]).argmax()] print("max error at sampled points: ", max(abs(err[0])), max(abs(err[1]))) @@ -194,7 +195,7 @@ def main(nelements): from pytential.symbolic.mappers import DerivativeTaker rep_pressure = stresslet_obj.apply_pressure(inv_sqrt_w_sigma, nvec_sym, mu_sym, qbx_forced_limit=-2) - pressure = bind((qbx, PointsTarget(eval_points_dev)), + pressure = bind((qbx, PointsTarget(eval_points_dev)), rep_pressure)(queue, sigma=sigma, mu=mu, normal=normal) pressure = pressure.get() print "pressure = ", pressure diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 9d114f19..9ca037df 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -46,6 +46,10 @@ __doc__ = """ # {{{ QBX layer potential source +class _not_provided: # noqa: N801 + pass + + class QBXLayerPotentialSource(LayerPotentialSourceBase): """A source discretization for a QBX layer potential. @@ -57,22 +61,25 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): # {{{ constructor / copy - def __init__(self, density_discr, fine_order, + def __init__(self, + density_discr, + fine_order, qbx_order=None, fmm_order=None, fmm_level_to_order=None, - target_stick_out_factor=1e-10, base_resampler=None, expansion_factory=None, + target_association_tolerance=_not_provided, # begin undocumented arguments # FIXME default debug=False once everything works debug=True, - refined_for_global_qbx=False, - expansion_disks_in_tree_have_extent=False, + _refined_for_global_qbx=False, + _expansion_disks_in_tree_have_extent=False, _expansion_disk_stick_out_factor=0, performance_data_file=None, - fmm_backend="sumpy"): + fmm_backend="sumpy", + target_stick_out_factor=_not_provided): """ :arg fine_order: The total degree to which the (upsampled) underlying quadrature is exact. @@ -85,6 +92,28 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): a reasonable(-ish?) default. """ + # {{{ argument processing + + if target_stick_out_factor is not _not_provided: + from warnings import warn + warn("target_stick_out_factor has been renamed to " + "target_association_tolerance. " + "Using target_stick_out_factor is deprecated " + "and will stop working in 2018.", + DeprecationWarning, stacklevel=2) + + if target_association_tolerance is not _not_provided: + raise TypeError("May not pass both target_association_tolerance and " + "target_stick_out_factor.") + + target_association_tolerance = target_stick_out_factor + + del target_stick_out_factor + + if target_association_tolerance is _not_provided: + target_association_tolerance = float( + np.finfo(density_discr.real_dtype).eps) * 1e3 + if fmm_level_to_order is None: if fmm_order is None and qbx_order is not None: fmm_order = qbx_order + 1 @@ -104,11 +133,16 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): def fmm_level_to_order(level): return fmm_order + # }}} + self.fine_order = fine_order self.qbx_order = qbx_order self.density_discr = density_discr self.fmm_level_to_order = fmm_level_to_order - self.target_stick_out_factor = target_stick_out_factor + + assert target_association_tolerance is not None + + self.target_association_tolerance = target_association_tolerance self.fmm_backend = fmm_backend # Default values are lazily provided if these are None @@ -120,9 +154,9 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): self.expansion_factory = expansion_factory self.debug = debug - self.refined_for_global_qbx = refined_for_global_qbx - self.expansion_disks_in_tree_have_extent = \ - expansion_disks_in_tree_have_extent + self._refined_for_global_qbx = _refined_for_global_qbx + self._expansion_disks_in_tree_have_extent = \ + _expansion_disks_in_tree_have_extent self._expansion_disk_stick_out_factor = _expansion_disk_stick_out_factor self.performance_data_file = performance_data_file @@ -132,14 +166,40 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): fine_order=None, qbx_order=None, fmm_level_to_order=None, - target_stick_out_factor=None, base_resampler=None, - _expansion_disk_stick_out_factor=None, + target_association_tolerance=_not_provided, + _expansion_disks_in_tree_have_extent=_not_provided, + _expansion_disk_stick_out_factor=_not_provided, performance_data_file=None, - debug=None, - refined_for_global_qbx=None, + debug=_not_provided, + _refined_for_global_qbx=None, + target_stick_out_factor=_not_provided, ): + + # {{{ argument processing + + if target_stick_out_factor is not _not_provided: + from warnings import warn + warn("target_stick_out_factor has been renamed to " + "target_association_tolerance. " + "Using target_stick_out_factor is deprecated " + "and will stop working in 2018.", + DeprecationWarning, stacklevel=2) + + if target_association_tolerance is not _not_provided: + raise TypeError("May not pass both target_association_tolerance and " + "target_stick_out_factor.") + + target_association_tolerance = target_stick_out_factor + + elif target_association_tolerance is _not_provided: + target_association_tolerance = self.target_association_tolerance + + del target_stick_out_factor + + # }}} + # FIXME Could/should share wrangler and geometry kernels # if no relevant changes have been made. return QBXLayerPotentialSource( @@ -149,23 +209,26 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): qbx_order=qbx_order if qbx_order is not None else self.qbx_order, fmm_level_to_order=( fmm_level_to_order or self.fmm_level_to_order), - target_stick_out_factor=( - target_stick_out_factor - if target_stick_out_factor is not None - else self.target_stick_out_factor), + target_association_tolerance=target_association_tolerance, base_resampler=base_resampler or self._base_resampler, debug=( - debug if debug is not None else self.debug), - refined_for_global_qbx=( - refined_for_global_qbx if refined_for_global_qbx is not None - else self.refined_for_global_qbx), - expansion_disks_in_tree_have_extent=( - self.expansion_disks_in_tree_have_extent), + # False is a valid value here + debug if debug is not _not_provided else self.debug), + _refined_for_global_qbx=( + # False is a valid value here + _refined_for_global_qbx + if _refined_for_global_qbx is not _not_provided + else self._refined_for_global_qbx), + _expansion_disks_in_tree_have_extent=( + # False is a valid value here + _expansion_disks_in_tree_have_extent + if _expansion_disks_in_tree_have_extent is not _not_provided + else self._expansion_disks_in_tree_have_extent), _expansion_disk_stick_out_factor=( # 0 is a valid value here _expansion_disk_stick_out_factor - if _expansion_disk_stick_out_factor is not None + if _expansion_disk_stick_out_factor is not _not_provided else self._expansion_disk_stick_out_factor), performance_data_file=( performance_data_file or self.performance_data_file), @@ -313,7 +376,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): return QBXFMMGeometryData(self.qbx_fmm_code_getter, self, target_discrs_and_qbx_sides, - target_stick_out_factor=self.target_stick_out_factor, + target_association_tolerance=self.target_association_tolerance, debug=self.debug) # }}} @@ -365,7 +428,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): from functools import partial oversample = partial(self.resampler, queue) - if not self.refined_for_global_qbx: + if not self._refined_for_global_qbx: from warnings import warn warn( "Executing global QBX without refinement. " diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 2c681b89..53f4de33 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -353,7 +353,7 @@ class QBXFMMGeometryData(object): def __init__(self, code_getter, lpot_source, target_discrs_and_qbx_sides, - target_stick_out_factor, debug): + target_association_tolerance, debug): """ .. rubric:: Constructor arguments @@ -368,7 +368,7 @@ class QBXFMMGeometryData(object): self.lpot_source = lpot_source self.target_discrs_and_qbx_sides = \ target_discrs_and_qbx_sides - self.target_stick_out_factor = target_stick_out_factor + self.target_association_tolerance = target_association_tolerance self.debug = debug @property @@ -493,7 +493,7 @@ class QBXFMMGeometryData(object): nparticles = nsources + target_info.ntargets target_radii = None - if self.lpot_source.expansion_disks_in_tree_have_extent: + if self.lpot_source._expansion_disks_in_tree_have_extent: target_radii = cl.array.zeros(queue, target_info.ntargets, self.coord_dtype) target_radii[:self.ncenters] = self.expansion_radii() @@ -545,7 +545,7 @@ class QBXFMMGeometryData(object): trav, _ = self.code_getter.build_traversal(queue, self.tree(), debug=self.debug) - if self.lpot_source.expansion_disks_in_tree_have_extent: + if self.lpot_source._expansion_disks_in_tree_have_extent: trav = trav.merge_close_lists(queue) return trav @@ -685,7 +685,8 @@ class QBXFMMGeometryData(object): # FIXME: try block... tgt_assoc_result = tgt_assoc(self.lpot_source, target_discrs_and_qbx_sides, - stick_out_factor=self.target_stick_out_factor) + target_association_tolerance=( + self.target_association_tolerance)) tree = self.tree() diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index 8efd9a88..aac3f0eb 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -624,7 +624,7 @@ def refine_for_global_qbx(lpot_source, code_container, # }}} - lpot_source = lpot_source.copy(debug=debug, refined_for_global_qbx=True) + lpot_source = lpot_source.copy(debug=debug, _refined_for_global_qbx=True) if len(connections) == 0: # FIXME: This is inefficient diff --git a/pytential/qbx/target_assoc.py b/pytential/qbx/target_assoc.py index cbf2a26c..22b2351b 100644 --- a/pytential/qbx/target_assoc.py +++ b/pytential/qbx/target_assoc.py @@ -178,7 +178,7 @@ QBX_CENTER_FINDER = AreaQueryElementwiseTemplate( particle_id_t center_offset, particle_id_t target_offset, particle_id_t *sorted_target_ids, - coord_t *expansion_radii_by_center_with_stick_out, + coord_t *expansion_radii_by_center_with_tolerance, coord_t *box_to_search_dist, int *target_flags, @@ -229,7 +229,7 @@ QBX_CENTER_FINDER = AreaQueryElementwiseTemplate( coord_t my_dist_to_center = distance(tgt_coords, center_coords); if (my_dist_to_center - <= expansion_radii_by_center_with_stick_out[center] + <= expansion_radii_by_center_with_tolerance[center] && my_dist_to_center < min_dist_to_center[i]) { target_status[i] = MARKED_QBX_CENTER_FOUND; @@ -480,7 +480,7 @@ class QBXTargetAssociator(object): def try_find_centers(self, queue, tree, peer_lists, lpot_source, target_status, target_flags, target_assoc, - stick_out_factor, debug, wait_for=None): + target_association_tolerance, debug, wait_for=None): # Round up level count--this gets included in the kernel as # a stack bound. Rounding avoids too many kernel versions. from pytools import div_ceil @@ -503,8 +503,8 @@ class QBXTargetAssociator(object): centers = [axis.with_queue(queue)[center_slice] for axis in tree.sources] expansion_radii_by_center = \ lpot_source._expansion_radii("ncenters").with_queue(queue) - expansion_radii_by_center_with_stick_out = \ - expansion_radii_by_center * (1 + stick_out_factor) + expansion_radii_by_center_with_tolerance = \ + expansion_radii_by_center * (1 + target_association_tolerance) # Idea: # @@ -516,7 +516,7 @@ class QBXTargetAssociator(object): queue, tree, centers, - expansion_radii_by_center_with_stick_out, + expansion_radii_by_center_with_tolerance, peer_lists, wait_for=wait_for) wait_for = [evt] @@ -537,7 +537,7 @@ class QBXTargetAssociator(object): tree.qbx_user_center_slice.start, tree.qbx_user_target_slice.start, tree.sorted_target_ids, - expansion_radii_by_center_with_stick_out, + expansion_radii_by_center_with_tolerance, box_to_search_dist, target_flags, target_status, @@ -655,7 +655,7 @@ class QBXTargetAssociator(object): return QBXTargetAssociation(target_to_center=target_to_center) def __call__(self, lpot_source, target_discrs_and_qbx_sides, - stick_out_factor=1e-10, debug=True, wait_for=None): + target_association_tolerance, debug=True, wait_for=None): """ Entry point for calling the target associator. @@ -732,7 +732,7 @@ class QBXTargetAssociator(object): self.try_find_centers(queue, tree, peer_lists, lpot_source, target_status, target_flags, target_assoc, - stick_out_factor, debug) + target_association_tolerance, debug) center_not_found = ( target_status == target_status_enum.MARKED_QBX_CENTER_PENDING) @@ -745,8 +745,8 @@ class QBXTargetAssociator(object): if (center_not_found & surface_target).any().get(): logger.warning("An on-surface target was not " "assigned a center. As a remedy you can try increasing " - "the \"stick_out_factor\" parameter, but this could " - "also cause an invalid center assignment.") + "the \"target_association_tolerance\" parameter, but " + "this could also cause an invalid center assignment.") refine_flags = cl.array.zeros(queue, tree.nqbxpanels, dtype=np.int32) have_panel_to_refine = self.mark_panels_for_refinement(queue, diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 994c9673..93e59c14 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -564,13 +564,13 @@ def run_int_eq_test(cl_ctx, queue, case, resolution): fplot = FieldPlotter( bbox_center, extent=2*2*bbox_size, npoints=(150, 150, 1)) - qbx_stick_out = qbx.copy(target_stick_out_factor=0.15) + qbx_tgt_tol = qbx.copy(target_association_tolerance=0.15) from pytential.target import PointsTarget from pytential.qbx import QBXTargetAssociationFailedException try: solved_pot = bind( - (qbx_stick_out, PointsTarget(fplot.points)), + (qbx_tgt_tol, PointsTarget(fplot.points)), op.representation(sym.var("u")) )(queue, u=u, k=case.k) except QBXTargetAssociationFailedException as e: @@ -1224,13 +1224,13 @@ def test_off_surface_eval_vs_direct(ctx_getter, do_plot=False): direct_qbx, _ = QBXLayerPotentialSource( pre_density_discr, 4*target_order, qbx_order, fmm_order=False, - target_stick_out_factor=0.05, + target_association_tolerance=0.05, ).with_refinement() fmm_qbx, _ = QBXLayerPotentialSource( pre_density_discr, 4*target_order, qbx_order, fmm_order=qbx_order + 3, expansion_disks_in_tree_have_extent=True, - target_stick_out_factor=0.05, + target_association_tolerance=0.05, ).with_refinement() fplot = FieldPlotter(np.zeros(2), extent=5, npoints=1000) diff --git a/test/test_stokes.py b/test/test_stokes.py index 02a71b1b..091d58fc 100644 --- a/test/test_stokes.py +++ b/test/test_stokes.py @@ -69,10 +69,11 @@ def run_exterior_stokes_2d(ctx_factory, nelements, InterpolatoryQuadratureSimplexGroupFactory(target_order)) from pytential.qbx import QBXLayerPotentialSource - stick_out = 0.05 + target_association_tolerance = 0.05 qbx, _ = QBXLayerPotentialSource( coarse_density_discr, fine_order=ovsmp_target_order, qbx_order=qbx_order, - fmm_order=fmm_order, target_stick_out_factor=stick_out + fmm_order=fmm_order, + target_association_tolerance=target_association_tolerance, ).with_refinement() density_discr = qbx.density_discr -- GitLab From e8acc41d8bc1e01338d9aee46fec87ae0c7fcd90 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 15 Jun 2017 17:12:49 -0500 Subject: [PATCH 5/9] Fix global QBX tests for new internal target assoc interface --- test/test_global_qbx.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py index c5f77476..642d5478 100644 --- a/test/test_global_qbx.py +++ b/test/test_global_qbx.py @@ -291,7 +291,8 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements): from pytential.qbx.target_assoc import QBXTargetAssociator target_assoc = ( - QBXTargetAssociator(cl_ctx)(lpot_source, target_discrs) + QBXTargetAssociator(cl_ctx)(lpot_source, target_discrs, + target_association_tolerance=1e-10) .get(queue=queue)) expansion_radii = lpot_source._expansion_radii("ncenters").get(queue) @@ -392,7 +393,8 @@ def test_target_association_failure(ctx_getter): QBXTargetAssociator, QBXTargetAssociationFailedException) with pytest.raises(QBXTargetAssociationFailedException): - QBXTargetAssociator(cl_ctx)(lpot_source, targets) + QBXTargetAssociator(cl_ctx)(lpot_source, targets, + target_association_tolerance=1e-10) # }}} -- GitLab From b20ed002752beb6d44b110e2b4d53958dafbbc0c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 15 Jun 2017 17:24:02 -0500 Subject: [PATCH 6/9] Switch ellipse eigenvalue tests to use zero-stick-out to recover accuracy issues from max_leaf_refine_weight change in bb129448c4e --- test/test_layer_pot.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 93e59c14..158d854a 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -153,8 +153,11 @@ def test_ellipse_eigenvalues(ctx_getter, ellipse_aspect, mode_nr, qbx_order): pre_density_discr = Discretization( cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) - qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4*target_order, - qbx_order, fmm_order=fmm_order).with_refinement() + qbx, _ = QBXLayerPotentialSource( + pre_density_discr, 4*target_order, + qbx_order, fmm_order=fmm_order, + _expansion_disks_in_tree_have_extent=True, + ).with_refinement() density_discr = qbx.density_discr nodes = density_discr.nodes().with_queue(queue) -- GitLab From 9cf9a5e6b6b6e61edd6dec1d13ea53910478bcef Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 15 Jun 2017 18:02:08 -0500 Subject: [PATCH 7/9] Switch LP identity tests to use zero-stick-out to recover accuracy issues from max_leaf_refine_weight change in bb129448c4e --- test/test_layer_pot.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 158d854a..8fd4172e 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -1060,7 +1060,9 @@ def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, qbx, _ = QBXLayerPotentialSource( pre_density_discr, 4*target_order, - qbx_order, fmm_order=qbx_order + order_bump + qbx_order, + fmm_order=qbx_order + order_bump, + _expansion_disks_in_tree_have_extent=True, ).with_refinement(**refiner_extra_kwargs) density_discr = qbx.density_discr -- GitLab From 65770cd6be92ba7f4544b8dee8b00713211a55c8 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 15 Jun 2017 20:25:00 -0500 Subject: [PATCH 8/9] Fix zero-stick-out parameter for off-surface eval test --- test/test_layer_pot.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 8fd4172e..34d11d3f 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -1234,7 +1234,7 @@ def test_off_surface_eval_vs_direct(ctx_getter, do_plot=False): fmm_qbx, _ = QBXLayerPotentialSource( pre_density_discr, 4*target_order, qbx_order, fmm_order=qbx_order + 3, - expansion_disks_in_tree_have_extent=True, + _expansion_disks_in_tree_have_extent=True, target_association_tolerance=0.05, ).with_refinement() -- GitLab From 1eb5e58afbe3a43f598843036370543dbf833d69 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 15 Jun 2017 22:57:21 -0500 Subject: [PATCH 9/9] Switch Stokes tests to use zero-stick-out to recover accuracy issues from max_leaf_refine_weight change in bb129448c4e --- test/test_stokes.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_stokes.py b/test/test_stokes.py index 091d58fc..e605f99a 100644 --- a/test/test_stokes.py +++ b/test/test_stokes.py @@ -74,6 +74,7 @@ def run_exterior_stokes_2d(ctx_factory, nelements, coarse_density_discr, fine_order=ovsmp_target_order, qbx_order=qbx_order, fmm_order=fmm_order, target_association_tolerance=target_association_tolerance, + _expansion_disks_in_tree_have_extent=True, ).with_refinement() density_discr = qbx.density_discr -- GitLab