diff --git a/examples/curve-pot.py b/examples/curve-pot.py index c1c113662..c7a4879ab 100644 --- a/examples/curve-pot.py +++ b/examples/curve-pot.py @@ -36,10 +36,6 @@ def process_kernel(knl, what_operator): elif what_operator == "D": from sumpy.kernel import DirectionalSourceDerivative source_knl = DirectionalSourceDerivative(knl) - # DirectionalTargetDerivative (temporarily?) removed - # elif what_operator == "S'": - # from sumpy.kernel import DirectionalTargetDerivative - # knl = DirectionalTargetDerivative(knl) else: raise RuntimeError(f"unrecognized operator '{what_operator}'") diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 184e0cadd..4697efb89 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -149,9 +149,12 @@ def coefficients_from_source(self, kernel, avec, bvec, rscale, sac=None): def evaluate(self, tgt_kernel, coeffs, bvec, rscale, sac=None): # no point in heeding rscale here--just ignore it - return sym.Add(*( - coeffs[self.get_storage_index(i)] / math.factorial(i) - for i in self.get_coefficient_identifiers())) + + # NOTE: We can't meaningfully apply target derivatives here. + # Instead, this is handled in LayerPotentialBase._evaluate. + return sym.Add(*( + coeffs[self.get_storage_index(i)] / math.factorial(i) + for i in self.get_coefficient_identifiers())) def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, dvec, tgt_rscale, sac=None, m2l_translation_classes_dependent_data=None): diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index ef440cf01..753a98819 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -100,10 +100,9 @@ def evaluate(self, kernel, coeffs, bvec, rscale, sac=None): rscale = 1 base_taker = kernel.get_derivative_taker(bvec, rscale, sac) - # Following is a no-op, but AxisTargetDerivative.postprocess_at_target and - # DirectionalTargetDerivative.postprocess_at_target only handles - # DifferentiatedExprDerivativeTaker and sympy expressions, so we need to - # make the taker a DifferentitatedExprDerivativeTaker instance. + # Following is a no-op, but AxisTargetDerivative.postprocess_at_target + # only handles DifferentiatedExprDerivativeTaker and sympy expressions, + # so we need to make the taker a DifferentitatedExprDerivativeTaker instance. base_taker = DifferentiatedExprDerivativeTaker(base_taker, {tuple([0]*self.dim): 1}) taker = kernel.postprocess_at_target(base_taker, bvec) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 81556bf92..d9dab1bca 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -79,7 +79,6 @@ .. autoclass:: AxisTargetDerivative .. autoclass:: AxisSourceDerivative .. autoclass:: DirectionalSourceDerivative -.. autoclass:: DirectionalTargetDerivative Transforming kernels -------------------- @@ -1143,74 +1142,6 @@ def __repr__(self): return f"{type(self).__name__}({self.inner_kernel!r}, {self.dir_vec_name})" -class DirectionalTargetDerivative(DirectionalDerivative): - directional_kind = "tgt" - target_array_name = "targets" - - def get_code_transformer(self): - from sumpy.codegen import VectorComponentRewriter - vcr = VectorComponentRewriter([self.dir_vec_name]) - from pymbolic.primitives import Variable - via = _VectorIndexAdder(self.dir_vec_name, (Variable("itgt"),)) - - inner_transform = self.inner_kernel.get_code_transformer() - - def transform(expr): - return via(vcr(inner_transform(expr))) - - return transform - - def postprocess_at_target(self, expr, bvec): - from sumpy.derivative_taker import ( - DifferentiatedExprDerivativeTaker, - diff_derivative_coeff_dict, - ) - from sumpy.symbolic import make_sym_vector as make_sympy_vector - dir_vec = make_sympy_vector(self.dir_vec_name, self.dim) - target_vec = make_sympy_vector(self.target_array_name, self.dim) - - expr = self.inner_kernel.postprocess_at_target(expr, bvec) - - # bvec = tgt - center - if not isinstance(expr, DifferentiatedExprDerivativeTaker): - result = 0 - for axis in range(self.dim): - # Since `bvec` and `tgt` are two different symbolic variables - # need to differentiate by both to get the correct answer - result += (expr.diff(bvec[axis]) + expr.diff(target_vec[axis])) \ - * dir_vec[axis] - return result - - new_transformation = defaultdict(lambda: 0) - for axis in range(self.dim): - axis_transformation = diff_derivative_coeff_dict( - expr.derivative_coeff_dict, axis, target_vec) - for mi, coeff in axis_transformation.items(): - new_transformation[mi] += coeff * dir_vec[axis] - - return DifferentiatedExprDerivativeTaker(expr.taker, - dict(new_transformation)) - - def get_args(self): - return [ - KernelArgument( - loopy_arg=lp.GlobalArg( - self.dir_vec_name, - None, - shape=(self.dim, "ntargets"), - offset=lp.auto - ), - ), - *self.inner_kernel.get_args() - ] - - def prepare_loopy_kernel(self, loopy_knl): - loopy_knl = self.inner_kernel.prepare_loopy_kernel(loopy_knl) - return lp.tag_array_axes(loopy_knl, self.dir_vec_name, "sep,C") - - mapper_method = "map_directional_target_derivative" - - class DirectionalSourceDerivative(DirectionalDerivative): directional_kind = "src" diff --git a/sumpy/qbx.py b/sumpy/qbx.py index ea3cded78..b3016f3a1 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -105,12 +105,12 @@ def _expand(self, sac, avec, bvec, rscale, isrc): def _evaluate(self, sac, avec, bvec, rscale, expansion_nr, coefficients): from sumpy.expansion.local import LineTaylorLocalExpansion tgt_knl = self.target_kernels[expansion_nr] - if isinstance(tgt_knl, LineTaylorLocalExpansion): + if isinstance(self.expansion, LineTaylorLocalExpansion): # In LineTaylorLocalExpansion.evaluate, we can't run # postprocess_at_target because the coefficients are assigned # symbols and postprocess with a derivative will make them zero. # Instead run postprocess here before the coefficients are assigned. - coefficients = [tgt_knl.postprocess_at_target(coeff, bvec) for + coefficients = [tgt_knl.postprocess_at_target(coeff, avec) for coeff in coefficients] assigned_coeffs = [ @@ -512,7 +512,6 @@ def find_jump_term(kernel, arg_provider): AxisTargetDerivative, DerivativeBase, DirectionalSourceDerivative, - DirectionalTargetDerivative, ) tgt_derivatives = [] @@ -522,9 +521,6 @@ def find_jump_term(kernel, arg_provider): if isinstance(kernel, AxisTargetDerivative): tgt_derivatives.append(kernel.axis) kernel = kernel.kernel - elif isinstance(kernel, DirectionalTargetDerivative): - tgt_derivatives.append(kernel.dir_vec_name) - kernel = kernel.kernel elif isinstance(kernel, AxisSourceDerivative): src_derivatives.append(kernel.axis) kernel = kernel.kernel diff --git a/sumpy/tools.py b/sumpy/tools.py index 8251c8bb6..14f76261a 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -298,7 +298,7 @@ def __init__(self, ctx: Any, device: Any | None = None) -> None: """ :arg target_kernels: list of :class:`~sumpy.kernel.Kernel` instances, - with :class:`sumpy.kernel.DirectionalTargetDerivative` as + with :class:`sumpy.kernel.AxisTargetDerivative` as the outermost kernel wrappers, if present. :arg source_kernels: list of :class:`~sumpy.kernel.Kernel` instances with :class:`~sumpy.kernel.DirectionalSourceDerivative` as the diff --git a/test/test_directionaltargetderivative.py b/test/test_directionaltargetderivative.py new file mode 100644 index 000000000..c484da6c1 --- /dev/null +++ b/test/test_directionaltargetderivative.py @@ -0,0 +1,172 @@ +from __future__ import annotations + +import numpy as np +import pytest +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.discretization import Discretization +from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory, +) +from meshmode.mesh.generation import make_curve_mesh, starfish +from pytential import GeometryCollection, bind, sym +from pytential.qbx import QBXLayerPotentialSource + +import pyopencl as cl +from arraycontext import flatten +from pytools.convergence import EOCRecorder + +from sumpy.expansion.local import LineTaylorLocalExpansion +from sumpy.kernel import AxisTargetDerivative, LaplaceKernel + + +def starfish_parametrization(t, n_arms=5, amplitude=0.25): + """ + Parametrization: + (x(θ), y(θ)) = r(θ)(cos(θ), sin(θ)), r(θ) = 1 + amplitude * sin(n_arms * θ). + It is used to compute normal vectors and expansion radius at different + refinement levels for arbitrary boundary targets. + """ + + theta = 2 * np.pi * t + + r = 1 + amplitude * np.sin(n_arms * theta) + dr_dt = amplitude * n_arms * 2 * np.pi * np.cos(n_arms * theta) + + x = r * np.cos(theta) + y = r * np.sin(theta) + + dx_dt = dr_dt * np.cos(theta) - r * np.sin(theta) * 2 * np.pi + dy_dt = dr_dt * np.sin(theta) + r * np.cos(theta) * 2 * np.pi + + jacobian_norm = np.sqrt(dx_dt**2 + dy_dt**2) + + tangent_x = dx_dt / jacobian_norm + tangent_y = dy_dt / jacobian_norm + + normal_x = tangent_y + normal_y = -tangent_x + + coords = np.vstack([x, y]) + tangents = np.vstack([tangent_x, tangent_y]) + normals = np.vstack([normal_x, normal_y]) + + return coords, tangents, normals, jacobian_norm + + +@pytest.mark.parametrize("kernel_type", ["laplace"]) +def test_lpot_dx_jump_relation_convergence(kernel_type): + """Test convergence of jump relations for single layer potential derivatives.""" + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + if kernel_type == "laplace": + knl = LaplaceKernel(2) + else: + raise ValueError(f"Unknown kernel type: {kernel_type}") + + qbx_order = 5 + nelements = [100, 150, 200] + target_order = 5 + upsampling_factor = 5 + + ntargets = 20 + rng = np.random.default_rng(42) + t = rng.uniform(0, 1, ntargets) + targets_h, _, targets_normals_h, jac = starfish_parametrization( + t, n_arms=5, amplitude=0.25 + ) + targets = actx.from_numpy(targets_h) + + from sumpy.qbx import LayerPotential + expansion = LineTaylorLocalExpansion(knl, qbx_order) + lplot_dx = LayerPotential( + actx.context, + expansion=expansion, + target_kernels=(AxisTargetDerivative(0, knl),), + source_kernels=(knl,) + ) + lplot_dy = LayerPotential( + actx.context, + expansion=expansion, + target_kernels=(AxisTargetDerivative(1, knl),), + source_kernels=(knl,) + ) + eocrec = EOCRecorder() + + for nelement in nelements: + mesh = make_curve_mesh(starfish, np.linspace(0, 1, nelement + 1), target_order) + pre_density_discr = Discretization( + actx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order) + ) + + qbx = QBXLayerPotentialSource( + pre_density_discr, + upsampling_factor * target_order, + qbx_order, + fmm_order=False + ) + places = GeometryCollection({"qbx": qbx}, auto_where=("qbx")) + + source_discr = places.get_discretization("qbx", sym.QBX_SOURCE_QUAD_STAGE2) + sources_h = actx.to_numpy(flatten(source_discr.nodes(), actx)).reshape(2, -1) + sources = actx.from_numpy(sources_h) + + dofdesc = sym.DOFDescriptor("qbx", sym.QBX_SOURCE_QUAD_STAGE2) + weights_nodes = bind( + places, + sym.weights_and_area_elements(ambient_dim=2, dim=1, dofdesc=dofdesc) + )(actx) + weights_nodes_h = actx.to_numpy(flatten(weights_nodes, actx)) + strengths = (actx.from_numpy(weights_nodes_h),) + + expansion_radii_h = jac / (2 * nelement) + centers_in = actx.from_numpy(targets_h - targets_normals_h * expansion_radii_h) + centers_out = actx.from_numpy(targets_h + targets_normals_h * expansion_radii_h) + + _, (eval_in_dx,) = lplot_dx( + actx.queue, + targets, sources, centers_in, strengths, + expansion_radii=expansion_radii_h + ) + + _, (eval_in_dy,) = lplot_dy( + actx.queue, + targets, sources, centers_in, strengths, + expansion_radii=expansion_radii_h + ) + + _, (eval_out_dx,) = lplot_dx( + actx.queue, + targets, sources, centers_out, strengths, + expansion_radii=expansion_radii_h + ) + + _, (eval_out_dy,) = lplot_dy( + actx.queue, + targets, sources, centers_out, strengths, + expansion_radii=expansion_radii_h + ) + + eval_in_dx = actx.to_numpy(eval_in_dx) + eval_in_dy = actx.to_numpy(eval_in_dy) + eval_out_dx = actx.to_numpy(eval_out_dx) + eval_out_dy = actx.to_numpy(eval_out_dy) + + eval_in = eval_in_dx * targets_normals_h[0] + \ + eval_in_dy * targets_normals_h[1] + eval_out = eval_out_dx * targets_normals_h[0] + \ + eval_out_dy * targets_normals_h[1] + + # check jump relation: S'_int - S'_ext = sigma (=1 for constant density) + jump_error = np.abs(eval_in - eval_out - 1) + + h_max = actx.to_numpy(bind(places, sym.h_max(places.ambient_dim))(actx)) + eocrec.add_data_point(h_max, np.max(jump_error)) + + assert eocrec.order_estimate() > qbx_order - 1 + + +if __name__ == "__main__": + pytest.main([__file__])