From 51dd80425d0ffc3106fda51f760a527f85ff0d74 Mon Sep 17 00:00:00 2001
From: Alexandru Fikl <alexfikl@gmail.com>
Date: Wed, 1 Apr 2020 15:15:39 -0500
Subject: [PATCH] kernel: add 3d Yukawa kernel

---
 sumpy/kernel.py   | 24 ++++++++++++++++++++----
 test/test_misc.py |  1 +
 2 files changed, 21 insertions(+), 4 deletions(-)

diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index 6383d8d5..53b8a550 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -559,14 +559,30 @@ class YukawaKernel(ExpressionKernel):
         """
         lam = var(yukawa_lambda_name)
 
-        if dim == 2:
-            r = pymbolic_real_norm_2(make_sym_vector("d", dim))
+        # NOTE: The Yukawa kernel is given by [1]
+        #   -1/(2 pi)**(n/2) * (lam/r)**(n/2-1) * K(n/2-1, lam r)
+        # where K is a modified Bessel function of the second kind.
+        #
+        # [1] https://en.wikipedia.org/wiki/Green%27s_function
+        # [2] http://dlmf.nist.gov/10.27#E8
+        # [3] https://dlmf.nist.gov/10.47#E2
+        # [4] https://dlmf.nist.gov/10.49
 
-            # http://dlmf.nist.gov/10.27#E8
+        r = pymbolic_real_norm_2(make_sym_vector("d", dim))
+        if dim == 2:
+            # NOTE: transform K(0, lam r) into a Hankel function using [2]
             expr = var("hankel_1")(0, var("I")*lam*r)
-            scaling_for_K0 = 1/2*var("pi")*var("I")  # noqa: N806
+            scaling_for_K0 = var("pi")/2*var("I")       # noqa: N806
 
             scaling = -1/(2*var("pi")) * scaling_for_K0
+        elif dim == 3:
+            # NOTE: to get the expression, we do the following and simplify
+            # 1. express K(1/2, lam r) as a modified spherical Bessel function
+            #   k(0, lam r) using [3] and use expression for k(0, lam r) from [4]
+            # 2. or use (AS 10.2.17) directly
+            expr = var("exp")(-lam*r) / r
+
+            scaling = -1/(4 * var("pi")**2)
         else:
             raise RuntimeError("unsupported dimensionality")
 
diff --git a/test/test_misc.py b/test/test_misc.py
index a25b2e80..77f887c1 100644
--- a/test/test_misc.py
+++ b/test/test_misc.py
@@ -68,6 +68,7 @@ class YukawaKernelInfo:
     BiharmonicKernelInfo(2),
     BiharmonicKernelInfo(3),
     YukawaKernelInfo(2, 5),
+    YukawaKernelInfo(3, 5),
     ])
 def test_pde_check_kernels(ctx_factory, knl_info, order=5):
     dim = knl_info.kernel.dim
-- 
GitLab