diff --git a/sumpy/e2e.py b/sumpy/e2e.py index 7243ea37be8e493e513396cf6f21d5a494357392..1d0a10bd1d5ee7d83a7a9554ef752c1a99d89d19 100644 --- a/sumpy/e2e.py +++ b/sumpy/e2e.py @@ -61,8 +61,8 @@ class E2EBase(KernelCacheWrapper): self.dim = src_expansion.dim def get_translation_loopy_insns(self): - from sumpy.symbolic import make_sym_vector - dvec = make_sym_vector("d", self.dim) + from sumpy.symbolic import make_sympy_vector + dvec = make_sympy_vector("d", self.dim) src_coeff_exprs = [sp.Symbol("src_coeff%d" % i) for i in xrange(len(self.src_expansion))] diff --git a/sumpy/e2p.py b/sumpy/e2p.py index 78c28c683d3f32a98be0077c198ebf3c17da7985..6c030f661f18a84c2423c5f5635178f3aa6f9810 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -59,8 +59,8 @@ class E2PBase(KernelCacheWrapper): assert tdr(knl) == expansion.kernel def get_loopy_insns_and_result_names(self): - from sumpy.symbolic import make_sym_vector - bvec = make_sym_vector("b", self.dim) + from sumpy.symbolic import make_sympy_vector + bvec = make_sympy_vector("b", self.dim) from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index f8e69506ddaee5efbfe7b2bb53ebc0330d6bf4b8..23ba8e5b22ab0bedce3e448ac772bd830771436d 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -31,7 +31,7 @@ class ExpansionBase(object): def __init__(self, kernel, order): # Don't be tempted to remove target derivatives here. # Line Taylor QBX can't do without them, because it can't - # take them afterwards. + # apply those derivatives to the expanded quantity. self.kernel = kernel self.order = order diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index bcb5d8a7ac67cdf31bb5bcd24b8c8ad85acf3f20..c89e4c8b9f499f822e339f868b484b4322a89434 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -49,8 +49,8 @@ class VolumeTaylorMultipoleExpansion( raise NotImplementedError("more than one source derivative " "not supported at present") - from sumpy.symbolic import make_sym_vector - dir_vec = make_sym_vector(kernel.dir_vec_name, kernel.dim) + from sumpy.symbolic import make_sympy_vector + dir_vec = make_sympy_vector(kernel.dir_vec_name, kernel.dim) coeff_identifiers = self.get_coefficient_identifiers() result = [0] * len(coeff_identifiers) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 19bc31f57fb307145ffca35b3ae78d4e2141f678..40cecbe72fd4c83a81ab1b7a478b242018476c9a 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -24,9 +24,11 @@ THE SOFTWARE. import loopy as lp -import sympy as sp import numpy as np from pymbolic.mapper import IdentityMapper +from sumpy.symbolic import pymbolic_real_norm_2 +from pymbolic.primitives import make_sym_vector +from pymbolic import var # {{{ basic kernel interface @@ -38,7 +40,7 @@ class Kernel(object): .. attribute:: dim *dim* is allowed to be *None* if the dimensionality is not yet - known. Use the :meth:`fix_dimensions` + known. """ def __init__(self, dim=None): @@ -91,13 +93,6 @@ class Kernel(object): return self._dim - def fix_dim(self, dim): - """Return a new :class:`Kernel` with :attr:`dim` set to - *dim*. - """ - - raise NotImplementedError - def get_base_kernel(self): return self @@ -154,52 +149,116 @@ class Kernel(object): """ return [] - def __sub__(self, other): - return DifferenceKernel(self, other) - # }}} -# {{{ "one" kernel (for testing) +class ExpressionKernel(Kernel): + is_complex_valued = False -class OneKernel(Kernel): - """A kernel whose value is 1 everywhere, for every source and target. - (for testing) - """ + init_arg_names = ("dim", "expression", "scaling", "is_complex_valued") - is_complex_valued = False + def __init__(self, dim, expression, scaling, is_complex_valued): + """ + :arg expression: A :mod:`pymbolic` expression depending on + variables *d_1* through *d_N* where *N* equals *dim*. + (These variables match what is returned from + :func:`pymbolic.primitives.make_sym_vector` with + argument `"d"`.) + :arg scaling: A :mod:`pymbolic` expression for the scaling + of the kernel. + """ - init_arg_names = ("dim",) + # expression and scaling are pymbolic objects because those pickle + # cleanly. D'oh, sympy! + + Kernel.__init__(self, dim) + + self.expression = expression + self.scaling = scaling + self.is_complex_valued = is_complex_valued def __getinitargs__(self): return (self._dim,) def __repr__(self): if self._dim is not None: - return "OneKnl%dD" % self.dim + return "ExprKnl%dD" % self.dim else: - return "OneKnl" + return "ExprKnl" def get_expression(self, dist_vec): - return sp.sympify(1) + if self.expression is None: + raise RuntimeError("expression in ExpressionKernel has not " + "been determined yet (this could be due to a PDE kernel " + "not having learned its dimensionality yet)") + + from sumpy.symbolic import PymbolicToSympyMapperWithSymbols + expr = PymbolicToSympyMapperWithSymbols()(self.expression) + + if self.dim != len(dist_vec): + raise ValueError("dist_vec length does not match expected dimension") + + expr = expr.subs([ + ("d_%d" % i, dist_vec_i) + for i, dist_vec_i in enumerate(dist_vec) + ]) + + return expr def get_scaling(self): """Return a global scaling of the kernel.""" - return 1 + if self.scaling is None: + raise RuntimeError("scaling in ExpressionKernel has not " + "been determined yet (this could be due to a PDE kernel " + "not having learned its dimensionality yet)") - mapper_method = "map_one_kernel" + from sumpy.symbolic import PymbolicToSympyMapperWithSymbols + return PymbolicToSympyMapperWithSymbols()(self.scaling) -# }}} + mapper_method = "map_expression_kernel" -# {{{ PDE kernels +one_kernel_2d = ExpressionKernel( + dim=2, + expression=1, + scaling=1, + is_complex_valued=False) +one_kernel_3d = ExpressionKernel( + dim=3, + expression=1, + scaling=1, + is_complex_valued=False) -class LaplaceKernel(Kernel): - is_complex_valued = False +# {{{ PDE kernels + +class LaplaceKernel(ExpressionKernel): init_arg_names = ("dim",) + def __init__(self, dim=None): + # See (Kress LIE, Thm 6.2) for scaling + if dim == 2: + r = pymbolic_real_norm_2(make_sym_vector("d", dim)) + expr = var("log")(r) + scaling = 1/(-2*var("pi")) + elif dim == 3: + r = pymbolic_real_norm_2(make_sym_vector("d", dim)) + expr = 1/r + scaling = 1/(4*var("pi")) + elif dim is None: + expr = None + scaling = None + else: + raise RuntimeError("unsupported dimensionality") + + ExpressionKernel.__init__( + self, + dim, + expression=expr, + scaling=scaling, + is_complex_valued=False) + def __getinitargs__(self): return (self._dim,) @@ -209,29 +268,6 @@ class LaplaceKernel(Kernel): else: return "LapKnl" - def get_expression(self, dist_vec): - assert self.dim == len(dist_vec) - from sumpy.symbolic import sympy_real_norm_2 - r = sympy_real_norm_2(dist_vec) - - if self.dim == 2: - return sp.log(r) - elif self.dim == 3: - return 1/r - else: - raise RuntimeError("unsupported dimensionality") - - def get_scaling(self): - """Return a global scaling of the kernel.""" - - # (Kress LIE, Thm 6.2) - if self.dim == 2: - return 1/(-2*sp.pi) - elif self.dim == 3: - return 1/(4*sp.pi) - else: - raise RuntimeError("unsupported dimensionality") - mapper_method = "map_laplace_kernel" @@ -240,7 +276,29 @@ class HelmholtzKernel(Kernel): def __init__(self, dim=None, helmholtz_k_name="k", allow_evanescent=False): - Kernel.__init__(self, dim) + k = var(helmholtz_k_name) + + if dim == 2: + r = pymbolic_real_norm_2(make_sym_vector("d", dim)) + expr = var("hankel_1")(0, k*r) + scaling = var("I")/4 + elif dim == 3: + r = pymbolic_real_norm_2(make_sym_vector("d", dim)) + expr = var("exp")(var("I")*k*r)/r + scaling = 1/(4*var("pi")) + elif dim is None: + expr = None + scaling = None + else: + raise RuntimeError("unsupported dimensionality") + + ExpressionKernel.__init__( + self, + dim, + expression=expr, + scaling=scaling, + is_complex_valued=True) + self.helmholtz_k_name = helmholtz_k_name self.allow_evanescent = allow_evanescent @@ -254,8 +312,6 @@ class HelmholtzKernel(Kernel): else: return "HelmKnl(%s)" % (self.helmholtz_k_name) - is_complex_valued = True - def prepare_loopy_kernel(self, loopy_knl): from sumpy.codegen import bessel_preamble_generator, bessel_mangler loopy_knl = lp.register_function_manglers(loopy_knl, @@ -265,31 +321,6 @@ class HelmholtzKernel(Kernel): return loopy_knl - def get_expression(self, dist_vec): - assert self.dim == len(dist_vec) - - from sumpy.symbolic import sympy_real_norm_2 - r = sympy_real_norm_2(dist_vec) - - k = sp.Symbol(self.helmholtz_k_name) - - if self.dim == 2: - return sp.Function("hankel_1")(0, k*r) - elif self.dim == 3: - return sp.exp(sp.I*k*r)/r - else: - raise RuntimeError("unsupported dimensionality") - - def get_scaling(self): - """Return a global scaling of the kernel.""" - - if self.dim == 2: - return sp.I/4 - elif self.dim == 3: - return 1/(4*sp.pi) - else: - raise RuntimeError("unsupported dimensionality") - def get_args(self): if self.allow_evanescent: k_dtype = np.complex128 @@ -300,41 +331,6 @@ class HelmholtzKernel(Kernel): mapper_method = "map_helmholtz_kernel" - -class DifferenceKernel(Kernel): - init_arg_names = ("kernel_plus", "kernel_minus") - - def __init__(self, kernel_plus, kernel_minus): - if (kernel_plus.get_base_kernel() is not kernel_plus - or kernel_minus.get_base_kernel() is not kernel_minus): - raise ValueError("kernels in difference kernel " - "must be basic, unwrapped PDE kernels") - - if kernel_plus.dim != kernel_minus.dim: - raise ValueError( - "kernels in difference kernel have different dim") - - self.kernel_plus = kernel_plus - self.kernel_minus = kernel_minus - - Kernel.__init__(self, self.kernel_plus.dim) - - # FIXME - raise NotImplementedError("difference kernel mostly unimplemented") - - def __getinitargs__(self): - return (self.kernel_plus, self.kernel_minus) - - def fix_dimensions(self, dim): - """Return a new :class:`Kernel` with :attr:`dim` set to - *dim*. - """ - return DifferenceKernel( - self.kernel_plus.fix_dimensions(dim), - self.kernel_minus.fix_dimensions(dim)) - - mapper_method = "map_difference_kernel" - # }}} @@ -481,8 +477,8 @@ class DirectionalTargetDerivative(DirectionalDerivative): dim = len(bvec) assert dim == self.dim - from sumpy.symbolic import make_sym_vector - dir_vec = make_sym_vector(self.dir_vec_name, dim) + from sumpy.symbolic import make_sympy_vector + dir_vec = make_sympy_vector(self.dir_vec_name, dim) # bvec = tgt-center return sum(dir_vec[axis]*expr.diff(bvec[axis]) @@ -507,8 +503,8 @@ class DirectionalSourceDerivative(DirectionalDerivative): dimensions = len(avec) assert dimensions == self.dim - from sumpy.symbolic import make_sym_vector - dir_vec = make_sym_vector(self.dir_vec_name, dimensions) + from sumpy.symbolic import make_sympy_vector + dir_vec = make_sympy_vector(self.dir_vec_name, dimensions) # avec = center-src -> minus sign from chain rule return sum(-dir_vec[axis]*expr.diff(avec[axis]) @@ -548,18 +544,11 @@ class KernelCombineMapper(KernelMapper): class KernelIdentityMapper(KernelMapper): - def map_laplace_kernel(self, kernel): + def map_expression_kernel(self, kernel): return kernel - def map_helmholtz_kernel(self, kernel): - return kernel - - def map_one_kernel(self, kernel): - return kernel - - def map_difference_kernel(self, kernel): - return DifferenceKernel( - self.rec(kernel.kernel_plus), self.rec(kernel.kernel_minus)) + map_laplace_kernel = map_expression_kernel + map_helmholtz_kernel = map_expression_kernel def map_axis_target_derivative(self, kernel): return AxisTargetDerivative(kernel.axis, self.rec(kernel.inner_kernel)) @@ -586,14 +575,11 @@ class DerivativeCounter(KernelCombineMapper): def combine(self, values): return max(values) - def map_laplace_kernel(self, kernel): - return 0 - - def map_helmholtz_kernel(self, kernel): + def map_expression_kernel(self, kernel): return 0 - def map_one_kernel(self, kernel): - return 0 + map_laplace_kernel = map_expression_kernel + map_helmholtz_kernel = map_expression_kernel def map_axis_target_derivative(self, kernel): return 1 + self.rec(kernel.inner_kernel) @@ -606,16 +592,31 @@ class KernelDimensionSetter(KernelIdentityMapper): def __init__(self, dim): self.dim = dim + def map_expression_kernel(self, kernel): + if kernel._dim is not None and self.dim != kernel.dim: + raise RuntimeError("cannot set kernel dimension to new value (%d)" + "different from existing one (%d)" + % (self.dim, kernel.dim)) + + return kernel + def map_laplace_kernel(self, kernel): + if kernel._dim is not None and self.dim != kernel.dim: + raise RuntimeError("cannot set kernel dimension to new value (%d)" + "different from existing one (%d)" + % (self.dim, kernel.dim)) + return LaplaceKernel(self.dim) def map_helmholtz_kernel(self, kernel): + if kernel._dim is not None and self.dim != kernel.dim: + raise RuntimeError("cannot set kernel dimension to new value (%d)" + "different from existing one (%d)" + % (self.dim, kernel.dim)) + return HelmholtzKernel(self.dim, kernel.helmholtz_k_name, kernel.allow_evanescent) - def map_one_kernel(self, kernel): - return OneKernel(self.dim) - # }}} diff --git a/sumpy/p2e.py b/sumpy/p2e.py index ce12d0f59e3910e8154de14dcb1d262927be26ca..ebb841f59b502009ca4a55a12a5a75b609ec769f 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -55,8 +55,8 @@ class P2EBase(KernelCacheWrapper): self.dim = expansion.dim def get_looy_instructions(self): - from sumpy.symbolic import make_sym_vector - avec = make_sym_vector("a", self.dim) + from sumpy.symbolic import make_sympy_vector + avec = make_sympy_vector("a", self.dim) from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() diff --git a/sumpy/p2p.py b/sumpy/p2p.py index fc7feeec2ff653c9458147001f99a5e574a2f5f1..b798039aa62b4b5383e50386389b0a3b77329c92 100644 --- a/sumpy/p2p.py +++ b/sumpy/p2p.py @@ -54,8 +54,8 @@ class P2PBase(KernelComputation, KernelCacheWrapper): self.dim = single_valued(knl.dim for knl in self.kernels) def get_loopy_insns_and_result_names(self): - from sumpy.symbolic import make_sym_vector - dvec = make_sym_vector("d", self.dim) + from sumpy.symbolic import make_sympy_vector + dvec = make_sympy_vector("d", self.dim) from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() diff --git a/sumpy/qbx.py b/sumpy/qbx.py index e5d0c24da3b4a623f1b9b31f77af293c882522cb..a7663f79519d4c0969911f864f01332ca912d914 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -81,10 +81,10 @@ class LayerPotentialBase(KernelComputation): @memoize_method def get_kernel(self): - from sumpy.symbolic import make_sym_vector + from sumpy.symbolic import make_sympy_vector - avec = make_sym_vector("a", self.dim) - bvec = make_sym_vector("b", self.dim) + avec = make_sympy_vector("a", self.dim) + bvec = make_sympy_vector("b", self.dim) from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() @@ -136,8 +136,8 @@ class LayerPotentialBase(KernelComputation): "{[isrc,itgt,idim]: 0<=itgt a[idim] = center[idim,itgt] - src[idim,isrc] {id=compute_a}", - "<> b[idim] = tgt[idim,itgt] - center[idim,itgt] {id=compute_b}", + "<> a[idim] = center[idim,itgt] - src[idim,isrc] {id=compute_a}", + "<> b[idim] = tgt[idim,itgt] - center[idim,itgt] {id=compute_b}", ]+self.get_kernel_scaling_assignments()+loopy_insns+[ lp.ExpressionInstruction(id=None, assignee="pair_result_%d" % i, expression=expr, diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py index 82c5890434e9ff0802826cd7f8f2cbc680638f23..b0166298729f48f0ecd4ba4f56777dae200b295e 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -24,7 +24,9 @@ THE SOFTWARE. import sympy as sp +import numpy as np from pymbolic.mapper import IdentityMapper as IdentityMapperBase +from pymbolic.sympy_interface import PymbolicToSympyMapper # {{{ trivial assignment elimination @@ -160,7 +162,12 @@ def sympy_real_norm_2(x): return sp.sqrt((x.T*x)[0, 0]) -def make_sym_vector(name, components): +def pymbolic_real_norm_2(x): + from pymbolic import var + return var("sqrt")(np.dot(x, x)) + + +def make_sympy_vector(name, components): return sp.Matrix( [sp.Symbol("%s%d" % (name, i)) for i in range(components)]) @@ -181,4 +188,14 @@ def find_power_of(base, prod): return 0 return result[power] + +class PymbolicToSympyMapperWithSymbols(PymbolicToSympyMapper): + def map_variable(self, expr): + if expr.name == "I": + return sp.I + elif expr.name == "pi": + return sp.pi + else: + return PymbolicToSympyMapper.map_variable(expr) + # vim: fdm=marker diff --git a/sumpy/version.py b/sumpy/version.py index 3247a986f662dabd59ec586cc988c8beb3bad360..d311f19b4bad064bda0fe5e8a23365145faa80d9 100644 --- a/sumpy/version.py +++ b/sumpy/version.py @@ -22,5 +22,5 @@ THE SOFTWARE. VERSION = (2014, 1) -VERSION_STATUS = "beta" +VERSION_STATUS = "beta1" VERSION_TEXT = ".".join(str(x) for x in VERSION) + VERSION_STATUS