From 5cca13c5ff18788bd2380e75849263596030e203 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 30 May 2017 18:01:31 -0400 Subject: [PATCH 1/4] Add biharmonic kernel --- sumpy/kernel.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index bb46b964..735ac0f6 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 -- GitLab From 23ed797ad94060cc0451d091c3c6d4d466d32903 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 30 May 2017 18:02:11 -0400 Subject: [PATCH 2/4] Some docs for the FieldPlotter --- sumpy/visualization.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/sumpy/visualization.py b/sumpy/visualization.py index 1bc9cbe3..7c8737bf 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 -- GitLab From 2d3ed2564a29d57e44c319c27be2e724a46a6054 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 30 May 2017 18:02:22 -0400 Subject: [PATCH 3/4] Code formatting --- sumpy/tools.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/sumpy/tools.py b/sumpy/tools.py index f5014295..56c32bf3 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 -- GitLab From 227f8c8c4605be16dd446e621f4ca0b8ba65cd1f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 30 May 2017 18:03:47 -0400 Subject: [PATCH 4/4] Add, document, test potential-differentiating calculus patch --- doc/eval.rst | 12 +++ doc/index.rst | 1 + sumpy/point_calculus.py | 163 ++++++++++++++++++++++++++++++++++++++++ test/test_misc.py | 92 +++++++++++++++++++++++ 4 files changed, 268 insertions(+) create mode 100644 doc/eval.rst create mode 100644 sumpy/point_calculus.py create mode 100644 test/test_misc.py diff --git a/doc/eval.rst b/doc/eval.rst new file mode 100644 index 00000000..100255a3 --- /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 8a97711e..516a51a9 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -17,6 +17,7 @@ Contents expansion interactions codegen + eval misc diff --git a/sumpy/point_calculus.py b/sumpy/point_calculus.py new file mode 100644 index 00000000..6117002f --- /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/test/test_misc.py b/test/test_misc.py new file mode 100644 index 00000000..30170804 --- /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 -- GitLab