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