From 11f1697650e2458cbd10362ffea2939cb5472a1e Mon Sep 17 00:00:00 2001
From: Matt Wala <wala1@illinois.edu>
Date: Fri, 19 May 2017 15:32:58 -0500
Subject: [PATCH] LineTaylorLocalExpansion: Reinstate the old derivative taking
behavior when using sympy, so code generation doesn't blow up.
Also add a test the uses LineTaylorLocalExpansion.
---
sumpy/expansion/local.py | 33 ++++++++++---
test/test_qbx.py | 104 +++++++++++++++++++++++++++++++++++++++
2 files changed, 129 insertions(+), 8 deletions(-)
create mode 100644 test/test_qbx.py
diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py
index f9ec106a..02a47fed 100644
--- a/sumpy/expansion/local.py
+++ b/sumpy/expansion/local.py
@@ -68,16 +68,33 @@ class LineTaylorLocalExpansion(LocalExpansionBase):
line_kernel = self.kernel.get_expression(avec_line)
- from sumpy.tools import MiDerivativeTaker, my_syntactic_subs
- deriv_taker = MiDerivativeTaker(line_kernel, (tau,))
+ from sumpy.symbolic import USE_SYMENGINE
+
+ if USE_SYMENGINE:
+ from sumpy.tools import MiDerivativeTaker, my_syntactic_subs
+ deriv_taker = MiDerivativeTaker(line_kernel, (tau,))
+
+ return [my_syntactic_subs(
+ self.kernel.postprocess_at_target(
+ self.kernel.postprocess_at_source(
+ deriv_taker.diff(i),
+ avec), bvec),
+ {tau: 0})
+ for i in self.get_coefficient_identifiers()]
+ else:
+ # Workaround for sympy. The automatic distribution after
+ # single-variable diff makes the expressions very large
+ # (https://github.com/sympy/sympy/issues/4596), so avoid doing
+ # single variable diff.
+ #
+ # See also https://gitlab.tiker.net/inducer/pytential/merge_requests/12
- return [my_syntactic_subs(
- self.kernel.postprocess_at_target(
+ return [self.kernel.postprocess_at_target(
self.kernel.postprocess_at_source(
- deriv_taker.diff(i),
- avec), bvec),
- {tau: 0})
- for i in self.get_coefficient_identifiers()]
+ line_kernel.diff("tau", i), avec),
+ bvec)
+ .subs("tau", 0)
+ for i in self.get_coefficient_identifiers()]
def evaluate(self, coeffs, bvec):
from pytools import factorial
diff --git a/test/test_qbx.py b/test/test_qbx.py
new file mode 100644
index 00000000..a6c2319b
--- /dev/null
+++ b/test/test_qbx.py
@@ -0,0 +1,104 @@
+from __future__ import division, absolute_import, print_function
+
+__copyright__ = "Copyright (C) 2017 Matt Wala"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import numpy as np
+import sys
+
+import pyopencl as cl
+
+import logging
+logger = logging.getLogger(__name__)
+
+
+try:
+ import faulthandler
+except ImportError:
+ pass
+else:
+ faulthandler.enable()
+
+
+def test_direct(ctx_getter):
+ # 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)
+
+ order = 12
+
+ from sumpy.qbx import LayerPotential
+ from sumpy.expansion.local import LineTaylorLocalExpansion
+
+ lpot = LayerPotential(ctx, [LineTaylorLocalExpansion(lknl, order)])
+
+ mode_nr = 25
+
+ from pytools.convergence import EOCRecorder
+
+ eocrec = EOCRecorder()
+
+ for n in [200, 300, 400]:
+ t = np.linspace(0, 2 * np.pi, n, endpoint=False)
+ unit_circle = np.exp(1j * t)
+ unit_circle = np.array([unit_circle.real, unit_circle.imag])
+
+ sigma = np.cos(mode_nr * t)
+ eigval = 1/(2*mode_nr)
+
+ result_ref = eigval * sigma
+
+ h = 2 * np.pi / n
+
+ targets = unit_circle
+ sources = unit_circle
+
+ radius = 7 * h
+ centers = unit_circle * (1 - radius)
+
+ strengths = (sigma * h,)
+ evt, (result_qbx,) = lpot(queue, targets, sources, centers, strengths)
+
+ eocrec.add_data_point(h, np.max(np.abs(result_ref - result_qbx)))
+
+ print(eocrec)
+
+ slack = 1.5
+ assert eocrec.order_estimate() > order - slack
+
+
+# You can test individual routines by typing
+# $ python test_kernels.py 'test_p2p(cl.create_some_context)'
+
+if __name__ == "__main__":
+ if len(sys.argv) > 1:
+ exec(sys.argv[1])
+ else:
+ from py.test.cmdline import main
+ main([__file__])
+
+# vim: fdm=marker
--
GitLab