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/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/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