Add Gauss transform kernel

Author: isuruf

sumpy/kernel.py

@@ -933,6 +933,49 @@ def get_pde_as_diff_op(self):
         return laplacian(w)
 
 
+class GaussTransformKernel(ExpressionKernel):
+    init_arg_names = ("dim",)
+
+    def __init__(self, dim, gauss_delta_name="delta"):
+        r = pymbolic_real_norm_2(make_sym_vector("d", dim))
+        delta = SpatialConstant(gauss_delta_name)
+        expr = var("exp")(-r**2/delta)
+        scaling = 1
+
+        super().__init__(
+                dim,
+                expression=expr,
+                global_scaling_const=scaling,
+                is_complex_valued=False)
+
+        self.gauss_delta_name = gauss_delta_name
+
+    def __getinitargs__(self):
+        return (self.dim, self.gauss_delta_name)
+
+    def update_persistent_hash(self, key_hash, key_builder):
+        key_hash.update(type(self).__name__.encode("utf8"))
+        key_builder.rec(key_hash, (self.dim, self.gauss_delta_name))
+
+    def __repr__(self):
+        return f"GaussKnl{self.dim}D"
+
+    def get_args(self):
+        return [
+                KernelArgument(
+                    loopy_arg=lp.ValueArg(self.gauss_delta_name, np.float64),
+                    )]
+
+    mapper_method = "map_gauss_transform_kernel"
+
+    def get_derivative_taker(self, dvec, rscale, sac):
+        """Return a :class:`sumpy.derivative_taker.ExprDerivativeTaker` instance
+        that supports taking derivatives of the base kernel with respect to dvec.
+        """
+        from sumpy.derivative_taker import RadialDerivativeTaker
+        return RadialDerivativeTaker(self.get_expression(dvec), dvec, rscale,
+                sac)
+
 # }}}
 
 
@@ -1350,6 +1393,7 @@ def map_expression_kernel(self, kernel):
     map_elasticity_kernel = map_expression_kernel
     map_line_of_compression_kernel = map_expression_kernel
     map_stresslet_kernel = map_expression_kernel
+    map_gauss_transform_kernel = map_expression_kernel
 
     def map_axis_target_derivative(self, kernel):
         return type(kernel)(kernel.axis, self.rec(kernel.inner_kernel))
@@ -1406,6 +1450,7 @@ def map_expression_kernel(self, kernel):
     map_yukawa_kernel = map_expression_kernel
     map_line_of_compression_kernel = map_expression_kernel
     map_stresslet_kernel = map_expression_kernel
+    map_gauss_transform_kernel = map_expression_kernel
 
     def map_axis_target_derivative(self, kernel):
         return 1 + self.rec(kernel.inner_kernel)

test/test_kernels.py

@@ -48,6 +48,7 @@
     HelmholtzKernel,
     BiharmonicKernel,
     StokesletKernel,
+    GaussTransformKernel,
     AxisTargetDerivative,
     DirectionalSourceDerivative)
 
@@ -488,6 +489,8 @@ def test_p2e2p(actx_factory, base_knl, expn_class, order, with_source_derivative
      VolumeTaylorMultipoleExpansion),
     (StokesletKernel(2, 0, 0), LinearPDEConformingVolumeTaylorLocalExpansion,
      LinearPDEConformingVolumeTaylorMultipoleExpansion),
+    (GaussTransformKernel(2), VolumeTaylorLocalExpansion,
+     VolumeTaylorMultipoleExpansion),
     ])
 def test_translations(actx_factory, knl, local_expn_class, mpole_expn_class,
         visualize=False):
@@ -506,6 +509,8 @@ def test_translations(actx_factory, knl, local_expn_class, mpole_expn_class,
         extra_kwargs["k"] = 0.05
     if isinstance(knl, StokesletKernel):
         extra_kwargs["mu"] = 0.05
+    if isinstance(knl, GaussTransformKernel):
+        extra_kwargs["delta"] = 2.0
 
     # Just to make sure things also work away from the origin
     rng = np.random.default_rng(18)