diff --git a/doc/eval.rst b/doc/eval.rst
new file mode 100644
index 0000000000000000000000000000000000000000..100255a39c42ead04dcba8c4e26eaef804e53a94
--- /dev/null
+++ b/doc/eval.rst
@@ -0,0 +1,12 @@
+Differentiation and Evaluation
+==============================
+
+Visualization of Potentials
+---------------------------
+
+.. automodule:: sympy.visualization
+
+Differentiation of Potentials
+-----------------------------
+
+.. automodule:: sympy.pde_check
diff --git a/doc/index.rst b/doc/index.rst
index 8a97711e08df824b705b8a19202d61a693d0e804..516a51a9217fa1c095d000f23a37a9b35e8bd326 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -17,6 +17,7 @@ Contents
expansion
interactions
codegen
+ eval
misc
diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index bb46b9647029c5f89719af3c998f4dc80e5ecd95..735ac0f6b0d7fb285ce81b8467f64b170b0419db 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -348,6 +348,39 @@ class LaplaceKernel(ExpressionKernel):
mapper_method = "map_laplace_kernel"
+class BiharmonicKernel(ExpressionKernel):
+ init_arg_names = ("dim",)
+
+ def __init__(self, dim=None):
+ r = pymbolic_real_norm_2(make_sym_vector("d", dim))
+ if dim == 2:
+ expr = r**2 * var("log")(r)
+ scaling = 1/(8*var("pi"))
+ elif dim == 3:
+ expr = r
+ scaling = 1 # FIXME: Unknown
+ else:
+ raise RuntimeError("unsupported dimensionality")
+
+ ExpressionKernel.__init__(
+ self,
+ dim,
+ expression=expr,
+ scaling=scaling,
+ is_complex_valued=False)
+
+ def __getinitargs__(self):
+ return (self._dim,)
+
+ def __repr__(self):
+ if self._dim is not None:
+ return "BiharmKnl%dD" % self.dim
+ else:
+ return "BiharmKnl"
+
+ mapper_method = "map_biharmonic_kernel"
+
+
class HelmholtzKernel(ExpressionKernel):
init_arg_names = ("dim", "helmholtz_k_name", "allow_evanescent")
@@ -792,6 +825,7 @@ class KernelIdentityMapper(KernelMapper):
return kernel
map_laplace_kernel = map_expression_kernel
+ map_biharmonic_kernel = map_expression_kernel
map_helmholtz_kernel = map_expression_kernel
map_stokeslet_kernel = map_expression_kernel
map_stresslet_kernel = map_expression_kernel
@@ -825,6 +859,7 @@ class DerivativeCounter(KernelCombineMapper):
return 0
map_laplace_kernel = map_expression_kernel
+ map_biharmonic_kernel = map_expression_kernel
map_helmholtz_kernel = map_expression_kernel
map_stokeslet_kernel = map_expression_kernel
map_stresslet_kernel = map_expression_kernel
@@ -837,6 +872,9 @@ class DerivativeCounter(KernelCombineMapper):
class KernelDimensionSetter(KernelIdentityMapper):
+ """Deprecated: This is no longer used and will be removed in 2018.
+ """
+
def __init__(self, dim):
self.dim = dim
diff --git a/sumpy/point_calculus.py b/sumpy/point_calculus.py
new file mode 100644
index 0000000000000000000000000000000000000000..6117002fb93562dab11f1240e55f5d86006cd4be
--- /dev/null
+++ b/sumpy/point_calculus.py
@@ -0,0 +1,163 @@
+from __future__ import division, absolute_import
+
+__copyright__ = "Copyright (C) 2017 Andreas Kloeckner"
+
+__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 numpy.linalg as la
+from pytools import memoize_method
+
+___doc___ = """
+.. autoclass:: CalculusPatch
+"""
+
+
+class CalculusPatch(object):
+ """Sets up a grid of points on which derivatives can be calculated. Useful
+ to verify that an evaluated potential actually solves a PDE.
+
+ .. attribute: dim
+
+ .. attribute: points
+
+ shape: ``(dim, npoints_total)``
+
+ .. autofunction:: diff
+ .. autofunction:: dx
+ .. autofunction:: dy
+ .. autofunction:: dy
+ .. autofunction:: laplace
+ .. autofunction:: eval_at_0
+ .. autoattribute:: x
+ .. autoattribute:: y
+ .. autoattribute:: z
+ """
+ def __init__(self, center, h=1e-3, order=4):
+ self.center = center
+
+ npoints = order + 1
+ points_1d = np.linspace(-h/2, h/2, npoints)
+ self._points_1d = points_1d
+
+ self.dim = dim = len(self.center)
+ self.center = center
+
+ points_shaped = np.array(np.meshgrid(
+ *[center[i] + points_1d for i in range(dim)],
+ indexing="ij"))
+
+ self._points_shaped = points_shaped
+ self.points = points_shaped.reshape(dim, -1)
+
+ self._pshape = points_shaped.shape[1:]
+
+ @memoize_method
+ def _vandermonde_1d(self):
+ points_1d = self._points_1d
+
+ npoints = len(self._points_1d)
+ vandermonde = np.zeros((npoints, npoints))
+ for i in range(npoints):
+ vandermonde[:, i] = points_1d**i
+
+ return vandermonde
+
+ @memoize_method
+ def _zero_eval_vec_1d(self):
+ # The zeroth coefficient--all others involve x=0.
+ return self._vandermonde_1d()[0]
+
+ @memoize_method
+ def _diff_mat_1d(self, nderivs):
+ npoints = len(self._points_1d)
+
+ vandermonde = self._vandermonde_1d()
+ coeff_diff_mat = np.diag(np.arange(1, npoints), 1)
+
+ n_diff_mat = np.eye(npoints)
+ for i in range(nderivs):
+ n_diff_mat = n_diff_mat.dot(coeff_diff_mat)
+
+ deriv_coeffs_mat = la.solve(vandermonde.T, n_diff_mat.T).T
+ return vandermonde.dot(deriv_coeffs_mat)
+
+ def diff(self, axis, f_values, nderivs=1):
+ """Return the derivative along *axis* of *f_values*.
+
+ :arg f_values: an array of shape ``(dim, npoints_total)``
+ :returns: an array of shape ``(dim, npoints_total)``
+ """
+
+ dim = len(self.center)
+
+ assert axis < dim
+
+ axes = "klmno"
+ src_axes = (axes[:axis] + "j" + axes[axis:])[:dim]
+ tgt_axes = (axes[:axis] + "i" + axes[axis:])[:dim]
+
+ return np.einsum(
+ "ij,%s->%s" % (src_axes, tgt_axes),
+ self._diff_mat_1d(nderivs),
+ f_values.reshape(*self._pshape)).reshape(-1)
+
+ def dx(self, f_values):
+ return self.diff(0, f_values)
+
+ def dy(self, f_values):
+ return self.diff(1, f_values)
+
+ def dz(self, f_values):
+ return self.diff(2, f_values)
+
+ def laplace(self, f_values):
+ """Return the Laplacian of *f_values*.
+
+ :arg f_values: an array of shape ``(dim, npoints_total)``
+ :returns: an array of shape ``(dim, npoints_total)``
+ """
+
+ return sum(self.diff(iaxis, f_values, 2) for iaxis in range(self.dim))
+
+ def eval_at_center(self, f_values):
+ """Interpolate *f_values* to the center point.
+
+ :arg f_values: an array of shape ``(dim, npoints_total)``
+ :returns: a scalar.
+ """
+ f_values = f_values.reshape(*self._pshape)
+ for i in range(self.dim):
+ f_values = self._zero_eval_vec_1d.dot(f_values)
+
+ return f_values
+
+ @property
+ def x(self):
+ return self.points[0]
+
+ @property
+ def y(self):
+ return self.points[1]
+
+ @property
+ def z(self):
+ return self.points[2]
diff --git a/sumpy/tools.py b/sumpy/tools.py
index f5014295c9f9f560470a85ec2f8d6a3dcc3a42ae..56c32bf3dbd80ecabcb02a553097deb8e9d85161 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -1,8 +1,4 @@
-from __future__ import division
-from __future__ import absolute_import
-import six
-from six.moves import range
-from six.moves import zip
+from __future__ import division, absolute_import
__copyright__ = "Copyright (C) 2012 Andreas Kloeckner"
@@ -26,6 +22,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
+import six
+from six.moves import range, zip
from pytools import memoize_method
import numpy as np
import sumpy.symbolic as sym
diff --git a/sumpy/visualization.py b/sumpy/visualization.py
index 1bc9cbe3ec96e484b5cce4baf3da71af9aef3ac9..7c8737bfc9b8a94598eeb0df035b896b0a8993c3 100644
--- a/sumpy/visualization.py
+++ b/sumpy/visualization.py
@@ -23,6 +23,10 @@ THE SOFTWARE.
"""
+___doc___ = """
+.. autoclass:: FieldPlotter
+"""
+
import numpy as np
from six.moves import range
@@ -59,6 +63,12 @@ def separate_by_real_and_imag(data, real_only):
class FieldPlotter:
+ """
+ .. automethod:: set_matplotlib_limits
+ .. automethod:: show_scalar_in_matplotlib
+ .. automethod:: show_scalar_in_mayavi
+ .. automethod:: write_vtk_file
+ """
def __init__(self, center, extent=1, npoints=1000):
center = np.asarray(center)
self.dimensions, = dim, = center.shape
diff --git a/test/test_misc.py b/test/test_misc.py
new file mode 100644
index 0000000000000000000000000000000000000000..301708043418ebd00da604d845ae8692418fa00e
--- /dev/null
+++ b/test/test_misc.py
@@ -0,0 +1,92 @@
+from __future__ import division, absolute_import, print_function
+
+__copyright__ = "Copyright (C) 2017 Andreas Kloeckner"
+
+__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 numpy.linalg as la
+import sys
+import sumpy.toys as t
+
+import pytest
+import pyopencl as cl # noqa: F401
+from pyopencl.tools import ( # noqa
+ pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+
+
+@pytest.mark.parametrize("dim", [2, 3])
+def test_pde_check_biharmonic(ctx_factory, dim, order=5):
+ from sumpy.kernel import BiharmonicKernel
+ tctx = t.ToyContext(ctx_factory(), BiharmonicKernel(dim))
+
+ pt_src = t.PointSources(
+ tctx,
+ np.random.rand(dim, 50) - 0.5,
+ np.ones(50))
+
+ from pytools.convergence import EOCRecorder
+ from sumpy.point_calculus import CalculusPatch
+ eoc_rec = EOCRecorder()
+
+ for h in [0.1, 0.05, 0.025]:
+ cp = CalculusPatch(np.array([1, 0, 0])[:dim], h=h, order=order)
+ pot = pt_src.eval(cp.points)
+
+ deriv = cp.laplace(cp.laplace(pot))
+
+ err = la.norm(deriv)
+ eoc_rec.add_data_point(h, err)
+
+ print(eoc_rec)
+ assert eoc_rec.order_estimate() > order-3-0.1
+
+
+@pytest.mark.parametrize("dim", [1, 2, 3])
+def test_pde_check(dim, order=4):
+ from sumpy.point_calculus import CalculusPatch
+ from pytools.convergence import EOCRecorder
+
+ for iaxis in range(dim):
+ eoc_rec = EOCRecorder()
+ for h in [0.1, 0.01, 0.001]:
+ cp = CalculusPatch(np.array([3, 0, 0])[:dim], h=h, order=order)
+ df_num = cp.diff(iaxis, np.sin(10*cp.points[iaxis]))
+ df_true = 10*np.cos(10*cp.points[iaxis])
+
+ err = la.norm(df_num-df_true)
+ eoc_rec.add_data_point(h, err)
+
+ print(eoc_rec)
+ assert eoc_rec.order_estimate() > order-2-0.1
+
+
+# You can test individual routines by typing
+# $ python test_misc.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