Add a test for target derivative

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}'")
 

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):

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)

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"
 

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

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

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__])