diff --git a/experiments/layerpot-coverage.py b/experiments/layerpot-coverage.py new file mode 100644 index 0000000000000000000000000000000000000000..6b7596249cef39c278e61a2a636c62509d7d2014 --- /dev/null +++ b/experiments/layerpot-coverage.py @@ -0,0 +1,143 @@ +from sumpy import set_caching_enabled +set_caching_enabled(False) + +import logging +if 0: + logging.basicConfig(level='INFO') + +import pyopencl as cl + +from functools import partial +import numpy as np +from meshmode.mesh.generation import make_curve_mesh + +def get_ellipse_mesh(resolution, aspect_ratio, mesh_order): + from meshmode.mesh.generation import ellipse + curve_func = partial(ellipse, aspect_ratio) + return make_curve_mesh(curve_func, + np.linspace(0, 1, resolution+1), + mesh_order) + +def run_test(cl_ctx, queue): + + q_order = 5 + qbx_order = q_order + fmm_backend = "sumpy" + mesh = get_ellipse_mesh(20, 40, mesh_order=5) + a = 1 + b = 1 / 40 + + if 0: + from meshmode.mesh.visualization import draw_curve + import matplotlib.pyplot as plt + draw_curve(mesh) + plt.axes().set_aspect('equal') + plt.show() + + from pytential.qbx import QBXLayerPotentialSource + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + pre_density_discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(q_order)) + + refiner_extra_kwargs = { + # "_expansion_disturbance_tolerance": 0.05, + "_scaled_max_curvature_threshold": 1, + "maxiter": 10, + } + + qbx, _ = QBXLayerPotentialSource( + pre_density_discr, + fine_order=4 * q_order, + qbx_order=qbx_order, + fmm_backend=fmm_backend, + fmm_order=qbx_order + 5, + ).with_refinement(**refiner_extra_kwargs) + + if 1: + print("%d stage-1 elements after refinement" + % qbx.density_discr.mesh.nelements) + print("%d stage-2 elements after refinement" + % qbx.stage2_density_discr.mesh.nelements) + print("quad stage-2 elements have %d nodes" + % qbx.quad_stage2_density_discr.groups[0].nunit_nodes) + + def reference_solu(rvec): + # a harmonic function + x, y = rvec + return 2.1 * x * y + (x**2 - y**2) * 0.5 + x + + bvals = reference_solu(qbx.density_discr.nodes().with_queue(queue)) + + from pytential.symbolic.pde.scalar import DirichletOperator + from sumpy.kernel import LaplaceKernel + from pytential import sym, bind + op = DirichletOperator(LaplaceKernel(2), -1) + + bound_op = bind( + qbx.copy(target_association_tolerance=0.5), + op.operator(sym.var('sigma'))) + rhs = bind(qbx.density_discr, op.prepare_rhs(sym.var("bc")))(queue, bc=bvals) + + from pytential.solve import gmres + gmres_result = gmres( + bound_op.scipy_op(queue, "sigma", dtype=np.float64), + rhs, + tol=1e-12, + progress=True, + hard_failure=True, + stall_iterations=50, no_progress_factor=1.05) + + + from sumpy.visualization import FieldPlotter + from pytential.target import PointsTarget + pltsize = b * 1.5 + fplot = FieldPlotter(np.array([-1 + pltsize * 0.5, 0]), extent=pltsize * 1.05, npoints=500) + plt_targets = cl.array.to_device(queue, fplot.points) + + interior_pts = (fplot.points[0]**2 / a**2 + fplot.points[1]**2 / b**2) < 0.99 + + exact_vals = reference_solu(fplot.points) + out_errs = [] + + for assotol in [0.05]: + + qbx_stick_out = qbx.copy(target_association_tolerance=0.05) + + vol_solution = bind((qbx_stick_out, PointsTarget(plt_targets)), + op.representation(sym.var('sigma')))( + queue, sigma=gmres_result.solution).get() + + + interior_error_linf = ( + np.linalg.norm( + np.abs(vol_solution - exact_vals)[interior_pts], ord=np.inf) + / + np.linalg.norm(exact_vals[interior_pts], ord=np.inf) + ) + + interior_error_l2 = ( + np.linalg.norm( + np.abs(vol_solution - exact_vals)[interior_pts], ord=2) + / + np.linalg.norm(exact_vals[interior_pts], ord=2) + ) + + print("\nassotol = %f" % assotol) + print("L_inf Error = %e " % interior_error_linf) + print("L_2 Error = %e " % interior_error_l2) + + out_errs.append(("error-%f" % assotol, np.abs(vol_solution - exact_vals))) + + if 1: + fplot.write_vtk_file("results.vts", out_errs) + + +if __name__ == '__main__': + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + run_test(cl_ctx, queue) + diff --git a/pytential/qbx/target_assoc.py b/pytential/qbx/target_assoc.py index 06fe1472608113ebf63f306757606648b4a50e60..d8f32f54fd9e72e8d3fdadd6612a2087ef839ecf 100644 --- a/pytential/qbx/target_assoc.py +++ b/pytential/qbx/target_assoc.py @@ -233,7 +233,12 @@ QBX_CENTER_FINDER = AreaQueryElementwiseTemplate( /* output */ int *target_to_center, - coord_t *min_dist_to_center, + int *target_to_center_plus, + int *target_to_center_minus, + coord_t *min_dist_to_center_plus, + coord_t *min_dist_to_center_minus, + coord_t *min_rel_dist_to_center_plus, + coord_t *min_rel_dist_to_center_minus, /* input, dim-dependent size */ %for ax in AXIS_NAMES[:dimensions]: @@ -272,17 +277,74 @@ QBX_CENTER_FINDER = AreaQueryElementwiseTemplate( coord_vec_t center_coords; ${load_particle( "INDEX_FOR_CENTER_PARTICLE(center)", "center_coords")} + coord_t my_dist_to_center = distance(tgt_coords, center_coords); - if (my_dist_to_center - <= expansion_radii_by_center_with_tolerance[center] - && my_dist_to_center < min_dist_to_center[i]) + // Relative distance is weighted by disk radius, or, + // equivalently, by the expanded radius (constant times disk radius). + coord_t my_rel_dist_to_center = my_dist_to_center / ( + expansion_radii_by_center_with_tolerance[center]); + + if (my_rel_dist_to_center > 1) + { + continue; + } + + // The idea is to use relative distance to determine the closest + // disk center on each side, and then pick the one with smaller + // absolute distance. + // + // Specifically, the following code does two things: + // + // 1. Find the center with minimal relative distance on either side + // 2. Find the side with minimal Euclidean distance from one of the + // choosen centers in step 1. + // + // min_dist_to_center_plus and min_dist_to_center_minus + // hold the absolute distances from minimum + // relative distance centers. + // + // Refer to: + // - https://gitlab.tiker.net/inducer/pytential/issues/132 + // - https://gitlab.tiker.net/inducer/pytential/merge_requests/181 + if (center_side > 0) + { + if (my_rel_dist_to_center < min_rel_dist_to_center_plus[i]) + { + target_status[i] = MARKED_QBX_CENTER_FOUND; + min_rel_dist_to_center_plus[i] = my_rel_dist_to_center; + min_dist_to_center_plus[i] = my_dist_to_center; + target_to_center_plus[i] = center; + } + } + else { - target_status[i] = MARKED_QBX_CENTER_FOUND; - min_dist_to_center[i] = my_dist_to_center; - target_to_center[i] = center; + if (my_rel_dist_to_center < min_rel_dist_to_center_minus[i]) + { + target_status[i] = MARKED_QBX_CENTER_FOUND; + min_rel_dist_to_center_minus[i] = my_rel_dist_to_center; + min_dist_to_center_minus[i] = my_dist_to_center; + target_to_center_minus[i] = center; + } } } + + // If the "winner" on either side is updated, + // bind to the newly found center if it is closer than the + // currently found one. + + if (min_dist_to_center_plus[i] < min_dist_to_center_minus[i]) + { + target_to_center[i] = target_to_center_plus[i]; + } + else + { + // This also includes the case where both distances are INFINITY, + // thus target_to_center_plus and target_to_center_minus are + // required to be initialized with the same default center + // -1 as make_default_target_association + target_to_center[i] = target_to_center_minus[i]; + } } """, name="find_centers", @@ -561,7 +623,7 @@ class TargetAssociationWrangler(TreeWranglerBase): # # (1) Tag leaf boxes around centers with max distance to usable center. # (2) Area query from targets with those radii to find closest eligible - # center. + # center in terms of relative distance. box_to_search_dist, evt = self.code_container.space_invader_query()( self.queue, @@ -572,11 +634,21 @@ class TargetAssociationWrangler(TreeWranglerBase): wait_for=wait_for) wait_for = [evt] - min_dist_to_center = cl.array.empty( - self.queue, tree.nqbxtargets, tree.coord_dtype) - min_dist_to_center.fill(np.inf) + def make_target_field(fill_val, dtype=tree.coord_dtype): + arr = cl.array.empty( + self.queue, tree.nqbxtargets, dtype) + arr.fill(fill_val) + wait_for.extend(arr.events) + return arr + + target_to_center_plus = make_target_field(-1, np.int32) + target_to_center_minus = make_target_field(-1, np.int32) + + min_dist_to_center_plus = make_target_field(np.inf) + min_dist_to_center_minus = make_target_field(np.inf) - wait_for.extend(min_dist_to_center.events) + min_rel_dist_to_center_plus = make_target_field(np.inf) + min_rel_dist_to_center_minus = make_target_field(np.inf) evt = knl( *unwrap_args( @@ -591,7 +663,12 @@ class TargetAssociationWrangler(TreeWranglerBase): target_flags, target_status, target_assoc.target_to_center, - min_dist_to_center, + target_to_center_plus, + target_to_center_minus, + min_dist_to_center_plus, + min_dist_to_center_minus, + min_rel_dist_to_center_plus, + min_rel_dist_to_center_minus, *tree.sources), range=slice(tree.nqbxtargets), queue=self.queue,