From de7c1fd2f08541a5223b6205fad0bccab79eddeb Mon Sep 17 00:00:00 2001
From: Alex Fikl <alexfikl@gmail.com>
Date: Sun, 24 Jun 2018 20:24:39 -0500
Subject: [PATCH 1/2] qbx: fix block generator with additional source vectors

---
 sumpy/qbx.py           |  3 +++
 test/test_matrixgen.py | 49 ++++++++++++++++++++++++++++--------------
 2 files changed, 36 insertions(+), 16 deletions(-)

diff --git a/sumpy/qbx.py b/sumpy/qbx.py
index 985565f1..b67c1042 100644
--- a/sumpy/qbx.py
+++ b/sumpy/qbx.py
@@ -354,6 +354,9 @@ class LayerPotentialMatrixBlockGenerator(LayerPotentialBase):
             self.get_kernel_scaling_assignments()
             + ["""
                 for imat
+                    <> itgt = tgtindices[imat]
+                    <> isrc = srcindices[imat]
+
                     <> a[idim] = center[idim, tgtindices[imat]] - \
                                  src[idim, srcindices[imat]] {dup=idim}
                     <> b[idim] = tgt[idim, tgtindices[imat]] - \
diff --git a/test/test_matrixgen.py b/test/test_matrixgen.py
index 892fef63..41a1ac09 100644
--- a/test/test_matrixgen.py
+++ b/test/test_matrixgen.py
@@ -26,12 +26,17 @@ import sys
 import numpy as np
 import numpy.linalg as la
 
-import pytest
 import pyopencl as cl
+import pyopencl.array
+
+from pytools.obj_array import make_obj_array
+
+from sumpy.tools import MatrixBlockIndex
+
+import pytest
 from pyopencl.tools import (  # noqa
         pytest_generate_tests_for_pyopencl as pytest_generate_tests)
 
-from sumpy.tools import MatrixBlockIndex
 
 import logging
 logger = logging.getLogger(__name__)
@@ -44,7 +49,7 @@ else:
     faulthandler.enable()
 
 
-def create_arguments(n, mode, target_radius=1.0):
+def _build_geometry(n, mode, target_radius=1.0):
     # parametrize circle
     t = np.linspace(0.0, 2.0 * np.pi, n, endpoint=False)
     unit_circle = np.exp(1j * t)
@@ -65,7 +70,7 @@ def create_arguments(n, mode, target_radius=1.0):
     return targets, sources, centers, sigma, expansion_radii
 
 
-def create_index_subset(nnodes, nblks, factor):
+def _build_block_index(nnodes, nblks, factor):
     indices = np.arange(0, nnodes)
     ranges = np.arange(0, nnodes + 1, nnodes // nblks)
 
@@ -86,15 +91,22 @@ def create_index_subset(nnodes, nblks, factor):
 
 
 @pytest.mark.parametrize('factor', [1.0, 0.6])
-def test_qbx_direct(ctx_getter, factor):
+@pytest.mark.parametrize('lpot_id', [1, 2])
+def test_qbx_direct(ctx_getter, factor, lpot_id):
     # This evaluates a single layer potential on a circle.
     logging.basicConfig(level=logging.INFO)
 
     ctx = ctx_getter()
     queue = cl.CommandQueue(ctx)
 
-    from sumpy.kernel import LaplaceKernel
-    lknl = LaplaceKernel(2)
+    from sumpy.kernel import LaplaceKernel, DirectionalSourceDerivative
+    if lpot_id == 1:
+        lknl = LaplaceKernel(2)
+    elif lpot_id == 2:
+        lknl = LaplaceKernel(2)
+        lknl = DirectionalSourceDerivative(lknl, dir_vec_name="dsource_vec")
+    else:
+        raise ValueError("unknow lpot_id")
 
     nblks = 10
     order = 12
@@ -112,21 +124,26 @@ def test_qbx_direct(ctx_getter, factor):
 
     for n in [200, 300, 400]:
         targets, sources, centers, sigma, expansion_radii = \
-                create_arguments(n, mode_nr)
+                _build_geometry(n, mode_nr)
 
         h = 2 * np.pi / n
         strengths = (sigma * h,)
 
-        tgtindices, tgtranges = create_index_subset(n, nblks, factor)
-        srcindices, srcranges = create_index_subset(n, nblks, factor)
+        tgtindices, tgtranges = _build_block_index(n, nblks, factor)
+        srcindices, srcranges = _build_block_index(n, nblks, factor)
         assert tgtranges.shape == srcranges.shape
 
+        knl_kwargs = {}
+        if lpot_id == 2:
+            knl_kwargs["dsource_vec"] = make_obj_array([
+                np.ones(n) for _ in range(2)])
+
         _, (mat,) = mat_gen(queue, targets, sources, centers,
-                expansion_radii)
+                expansion_radii, **knl_kwargs)
         result_mat = mat.dot(strengths[0])
 
         _, (result_lpot,) = lpot(queue, targets, sources, centers, strengths,
-                expansion_radii)
+                expansion_radii, **knl_kwargs)
 
         eps = 1.0e-10 * la.norm(result_lpot)
         assert la.norm(result_mat - result_lpot) < eps
@@ -134,7 +151,7 @@ def test_qbx_direct(ctx_getter, factor):
         index_set = MatrixBlockIndex(queue,
             tgtindices, srcindices, tgtranges, srcranges)
         _, (blk,) = blk_gen(queue, targets, sources, centers, expansion_radii,
-                            index_set)
+                            index_set, **knl_kwargs)
 
         rowindices, colindices = index_set.linear_indices()
         eps = 1.0e-10 * la.norm(mat)
@@ -163,13 +180,13 @@ def test_p2p_direct(ctx_getter, exclude_self, factor):
 
     for n in [200, 300, 400]:
         targets, sources, _, sigma, _ = \
-            create_arguments(n, mode_nr, target_radius=1.2)
+            _build_geometry(n, mode_nr, target_radius=1.2)
 
         h = 2 * np.pi / n
         strengths = (sigma * h,)
 
-        tgtindices, tgtranges = create_index_subset(n, nblks, factor)
-        srcindices, srcranges = create_index_subset(n, nblks, factor)
+        tgtindices, tgtranges = _build_block_index(n, nblks, factor)
+        srcindices, srcranges = _build_block_index(n, nblks, factor)
         assert tgtranges.shape == srcranges.shape
 
         extra_kwargs = {}
-- 
GitLab


From 6da3039b443ecfa921f7fa97c2e63f2f373c1624 Mon Sep 17 00:00:00 2001
From: Alex Fikl <alexfikl@gmail.com>
Date: Sun, 24 Jun 2018 20:52:10 -0500
Subject: [PATCH 2/2] flake8: unused import

---
 test/test_matrixgen.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/test_matrixgen.py b/test/test_matrixgen.py
index 41a1ac09..5fde01ed 100644
--- a/test/test_matrixgen.py
+++ b/test/test_matrixgen.py
@@ -27,7 +27,6 @@ import numpy as np
 import numpy.linalg as la
 
 import pyopencl as cl
-import pyopencl.array
 
 from pytools.obj_array import make_obj_array
 
-- 
GitLab