diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py index f7cd6559f06878bbfc7811fec1f6c5d814afc193..8393f744caae96c798dd293be299b16127df1dc4 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 c1f4719c66ca60a60f97413a8c8aa455b565bb1b..b591c5bee25db99c69fb7d69c1dc243d352692cf 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 5e6bbcf43452b0354cc90bfce0d2e8f9214c4b3e..111265cb358cb7eb4b3b0671518c2954afb42236 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 b3f016a8109bb3eb9070d77c9515ea85d5df1f48..ed8c6106fcde56e7919344fb44b9500f7013aed0 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 cf8cf166fa96c6457a59cc0c93c4bb21b6620d2f..183fc915cb286ccb7731a76e9c1ed6cc3869efd5 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 4f6e7bee4158fc39d0a9f00ad9f276c22604e905..977ae8329991e5f5acabdb1bbec3024fa55ac39e 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 0c4bfdf366d9adb96eea1cdfe298e9345afac07e..9ca037df724cad76ee1d20705f541cc325c4f847 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,21 +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. @@ -84,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 @@ -103,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 @@ -119,9 +154,10 @@ 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 def copy( @@ -130,12 +166,40 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): fine_order=None, qbx_order=None, fmm_level_to_order=None, - target_stick_out_factor=None, base_resampler=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( @@ -145,20 +209,29 @@ 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), - performance_data_file=self.performance_data_file, + # 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 _not_provided + else self._expansion_disk_stick_out_factor), + performance_data_file=( + performance_data_file or self.performance_data_file), fmm_backend=self.fmm_backend) # }}} @@ -303,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) # }}} @@ -355,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/fmm.py b/pytential/qbx/fmm.py index c7d6c029448ce9c5c47ebedfc92a91a509a35157..6afa047b692e02a536d93def31f5c08097bcc587 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) - start, end = traversal.sep_smaller_starts[itgt_box:itgt_box+2] + outf.write("mp_eval = {cost}\n" + .format(cost=nmp_eval * p_fmm)) - nmp_eval += ntargets * (end-start) - - 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,71 +653,86 @@ 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)) + + process_form_qbxl() # }}} # {{{ 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] + + start, stop = ( + ssn.starts[icontaining_tgt_box], + ssn.starts[icontaining_tgt_box+1]) - for itgt_center, tgt_icenter in enumerate(global_qbx_centers): - itgt_box = qbx_center_to_target_box[tgt_icenter] + nqbx_m2l += stop-start - start, end = traversal.sep_smaller_starts[itgt_box:itgt_box+2] - nqbx_m2l += end - 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() # }}} diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index f2669274bcd018cbc33f5bda5dbd338b173c0446..53f4de33eeab5c8ad618651bb6bd066b4f676ea3 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() @@ -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: @@ -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 8efd9a8800f614f9d54a5d30285f108c42072ba1..aac3f0ebf5f9a39e9cbaec212f377e68c91f5ad1 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 cbf2a26c62c4936cf729b50ee96c088212d1f38f..22b2351ba82ce692cca49b5b73e1d1730604924b 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_global_qbx.py b/test/test_global_qbx.py index c5f77476b908ab7e5681784e6de73fa647ef428c..642d5478e6caff1abd717804f1706acaeb9cfb7d 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) # }}} diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 994c96735834eb6c76f36956731a1e84644c1f12..34d11d3f2924be38c109467a68d87fbf4617dbe0 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) @@ -564,13 +567,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: @@ -1057,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 @@ -1224,13 +1229,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, + _expansion_disks_in_tree_have_extent=True, + 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 02a71b1bc09b8a727c81a2c8f3a4e8b424afce7e..e605f99a5c0737027608adcd53a8163ccdfbd1fd 100644 --- a/test/test_stokes.py +++ b/test/test_stokes.py @@ -69,10 +69,12 @@ 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, + _expansion_disks_in_tree_have_extent=True, ).with_refinement() density_discr = qbx.density_discr