From 691c0fbedb6c24fe3778cb89015c0c48b4ba95d7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 11:05:26 -0500 Subject: [PATCH 01/60] Minor code cleanups --- sumpy/__init__.py | 3 +-- sumpy/codegen.py | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/sumpy/__init__.py b/sumpy/__init__.py index b7d4293e..c770e802 100644 --- a/sumpy/__init__.py +++ b/sumpy/__init__.py @@ -1,5 +1,4 @@ -from __future__ import division -from __future__ import absolute_import +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2013 Andreas Kloeckner" diff --git a/sumpy/codegen.py b/sumpy/codegen.py index ca608690..5d27708f 100644 --- a/sumpy/codegen.py +++ b/sumpy/codegen.py @@ -477,7 +477,7 @@ class PowerRewriter(CSECachingMapperMixin, IdentityMapper): if exp > 1 and exp % 2 == 0: square = prim.wrap_in_cse(new_base*new_base) return self.rec(prim.wrap_in_cse(square**(exp//2))) - if exp > 1 and exp % 2 == 1: + elif exp > 1 and exp % 2 == 1: square = prim.wrap_in_cse(new_base*new_base) return self.rec(prim.wrap_in_cse(square**((exp-1)//2))*new_base) elif exp == 1: -- GitLab From 57caeabb394fbc79da409ce360f2a6cb71b8e478 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 11:07:18 -0500 Subject: [PATCH 02/60] Another minor cleanup --- sumpy/expansion/local.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 76f3d1eb..a39fe99b 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -1,5 +1,4 @@ from __future__ import division, absolute_import -from six.moves import range, zip __copyright__ = "Copyright (C) 2012 Andreas Kloeckner" @@ -23,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from six.moves import range, zip import sumpy.symbolic as sym from sumpy.expansion import ( -- GitLab From fb5eb7fcea47bea73c175dc0b3f069c07212934f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 13:05:28 -0500 Subject: [PATCH 03/60] Minor doc fix for sumpy.point_calculus --- sumpy/point_calculus.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/point_calculus.py b/sumpy/point_calculus.py index 2706aa98..e938a9d8 100644 --- a/sumpy/point_calculus.py +++ b/sumpy/point_calculus.py @@ -46,7 +46,7 @@ class CalculusPatch(object): .. automethod:: dy .. automethod:: dy .. automethod:: laplace - .. automethod:: eval_at_0 + .. automethod:: eval_at_center .. autoattribute:: x .. autoattribute:: y .. autoattribute:: z -- GitLab From 9727e70457fa4b360bfdc2f2df1bcfa3a1769afa Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 13:05:47 -0500 Subject: [PATCH 04/60] Fix doc build with python from virtualenv --- doc/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/Makefile b/doc/Makefile index cae96bdb..cd9f9109 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -3,7 +3,7 @@ # You can set these variables from the command line. SPHINXOPTS = -SPHINXBUILD = sphinx-build +SPHINXBUILD = python `which sphinx-build` PAPER = BUILDDIR = _build -- GitLab From 59b9b1cef64e1cb6df44853675337ef156fb02fb Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 13:06:17 -0500 Subject: [PATCH 05/60] Minor post-modernize cleanup --- sumpy/expansion/__init__.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 08754a10..4c76ea45 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -1,6 +1,4 @@ -from __future__ import division -from __future__ import absolute_import -from six.moves import range +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2012 Andreas Kloeckner" @@ -24,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from six.moves import range import numpy as np import logging from pytools import memoize_method -- GitLab From 73598e259dbe0739057739ded2398abe4101128e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 15:39:07 -0500 Subject: [PATCH 06/60] Make expansion toys work with target derivative kernels --- sumpy/toys.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/sumpy/toys.py b/sumpy/toys.py index 8caa657c..e90bf5dc 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -27,6 +27,7 @@ from six.moves import range # noqa: F401 from pytools import memoize_method from numbers import Number +from sumpy.kernel import TargetDerivativeRemover import numpy as np # noqa: F401 import loopy as lp # noqa: F401 @@ -74,7 +75,8 @@ class ToyContext(object): def get_p2m(self, order): from sumpy import P2EFromSingleBox return P2EFromSingleBox(self.cl_context, - self.mpole_expn_class(self.kernel, order), + self.mpole_expn_class( + TargetDerivativeRemover()(self.kernel), order), [self.kernel]) @memoize_method @@ -88,7 +90,8 @@ class ToyContext(object): def get_m2p(self, order): from sumpy import E2PFromSingleBox return E2PFromSingleBox(self.cl_context, - self.mpole_expn_class(self.kernel, order), + self.mpole_expn_class( + TargetDerivativeRemover()(self.kernel), order), [self.kernel]) @memoize_method -- GitLab From 28a552531dc06a4a48c6d6acc0d596c811736170 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 15:40:07 -0500 Subject: [PATCH 07/60] Get rscale to work with multipole expansions, target derivatives, and Laplace and Helmholtz --- examples/coeff-scaling-experiment.py | 55 +++++++ sumpy/e2p.py | 16 +- sumpy/expansion/__init__.py | 6 +- sumpy/expansion/multipole.py | 27 ++-- sumpy/kernel.py | 224 ++++++++++++++++++--------- sumpy/p2e.py | 16 +- sumpy/p2p.py | 8 +- sumpy/tools.py | 12 +- sumpy/toys.py | 23 +-- sumpy/version.py | 2 +- 10 files changed, 287 insertions(+), 102 deletions(-) create mode 100644 examples/coeff-scaling-experiment.py diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py new file mode 100644 index 00000000..3c14dd45 --- /dev/null +++ b/examples/coeff-scaling-experiment.py @@ -0,0 +1,55 @@ +import pyopencl as cl +import sumpy.toys as t +import numpy as np +from sumpy.visualization import FieldPlotter +import matplotlib.pyplot as plt + + +# FIXME: Get rid of this once everything is working + +def main(): + dim = 2 + + from sumpy.kernel import ( # noqa: F401 + YukawaKernel, HelmholtzKernel, LaplaceKernel, + AxisTargetDerivative) + tctx = t.ToyContext( + cl.create_some_context(), + #LaplaceKernel(dim), + #AxisTargetDerivative(0, LaplaceKernel(dim)), + #YukawaKernel(dim), extra_source_kwargs={"lam": 5}, + HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3}, + ) + + np.random.seed(12) + scale = 2**(-14) + #scale = 1 + pts = np.random.rand(dim, 50) - 0.5 + pt_src = t.PointSources( + tctx, + scale * pts, + np.ones(50)) + + mctr = scale*np.array([0., 0, 0])[:dim] + mexp = t.multipole_expand(pt_src, mctr, order=4, rscale=scale) + + lctr = scale*np.array([2.5, 0, 0])[:dim] + #lexp = t.local_expand(mexp, lctr) + print(mexp.coeffs) + #print(lexp.coeffs) + + diff = mexp - pt_src + + diag = np.sqrt(dim) + print(t.l_inf(diff, scale*0.5*diag, center=lctr) + / t.l_inf(pt_src, scale*0.5*diag, center=lctr)) + if 0: + fp = FieldPlotter(lexp.center, extent=scale*0.5) + + t.logplot(fp, diff, cmap="jet", vmin=-16, vmax=-7) + plt.colorbar() + plt.show() + + +if __name__ == "__main__": + main() diff --git a/sumpy/e2p.py b/sumpy/e2p.py index 6a77a5cc..ea5c382f 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -77,12 +77,15 @@ class E2PBase(KernelCacheWrapper): from sumpy.symbolic import make_sym_vector bvec = make_sym_vector("b", self.dim) + import sympy as sp + rscale = sp.Symbol("rscale") + from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() coeff_exprs = [sym.Symbol("coeff%d" % i) for i in range(len(self.expansion.get_coefficient_identifiers()))] - value = self.expansion.evaluate(coeff_exprs, bvec) + value = self.expansion.evaluate(coeff_exprs, bvec, rscale) result_names = [ sac.assign_unique("result_%d_p" % i, @@ -108,7 +111,8 @@ class E2PBase(KernelCacheWrapper): sympy_conv = SympyToPymbolicMapper() return [lp.Assignment(id=None, assignee="kernel_scaling", - expression=sympy_conv(self.expansion.kernel.get_scaling()), + expression=sympy_conv( + self.expansion.kernel.get_global_scaling_const()), temp_var_type=lp.auto)] def get_cache_key(self): @@ -165,6 +169,7 @@ class E2PFromSingleBox(E2PBase): lp.GlobalArg("box_target_starts,box_target_counts_nonchild", None, shape=None), lp.GlobalArg("centers", None, shape="dim, naligned_boxes"), + lp.ValueArg("rscale", None), lp.GlobalArg("result", None, shape="nresults, ntargets", dim_tags="sep,C"), lp.GlobalArg("src_expansions", None, @@ -205,7 +210,12 @@ class E2PFromSingleBox(E2PBase): """ knl = self.get_cached_optimized_kernel() - return knl(queue, **kwargs) + centers = kwargs.pop("centers") + # "1" may be passed for rscale, which won't have its type + # meaningfully inferred. Make the type of rscale explicit. + rscale = centers.dtype.type(kwargs.pop("rscale")) + + return knl(queue, centers=centers, rscale=rscale, **kwargs) # }}} diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 4c76ea45..38107482 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -494,8 +494,7 @@ class DefaultExpansionFactory(ExpansionFactoryBase): and base_kernel.dim == 2): from sumpy.expansion.local import Y2DLocalExpansion return Y2DLocalExpansion - elif (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel) - and base_kernel.dim == 3): + elif isinstance(base_kernel.get_base_kernel(), HelmholtzKernel): from sumpy.expansion.local import \ HelmholtzConformingVolumeTaylorLocalExpansion return HelmholtzConformingVolumeTaylorLocalExpansion @@ -523,8 +522,7 @@ class DefaultExpansionFactory(ExpansionFactoryBase): from sumpy.expansion.multipole import ( LaplaceConformingVolumeTaylorMultipoleExpansion) return LaplaceConformingVolumeTaylorMultipoleExpansion - elif (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel) - and base_kernel.dim == 3): + elif isinstance(base_kernel.get_base_kernel(), HelmholtzKernel): from sumpy.expansion.multipole import ( HelmholtzConformingVolumeTaylorMultipoleExpansion) return HelmholtzConformingVolumeTaylorMultipoleExpansion diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 958bca6d..50d18142 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -25,6 +25,7 @@ THE SOFTWARE. from six.moves import range, zip import sumpy.symbolic as sym # noqa +from sumpy.tools import sympy_vec_subs from sumpy.expansion import ( ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, HelmholtzConformingVolumeTaylorExpansion) @@ -53,12 +54,14 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): Coefficients represent the terms in front of the kernel derivatives. """ - def coefficients_from_source(self, avec, bvec): + def coefficients_from_source(self, avec, bvec, rscale): from sumpy.kernel import DirectionalSourceDerivative kernel = self.kernel from sumpy.tools import mi_power, mi_factorial + avec = avec/rscale + if isinstance(kernel, DirectionalSourceDerivative): if kernel.get_base_kernel() is not kernel.inner_kernel: raise NotImplementedError("more than one source derivative " @@ -90,18 +93,24 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): return ( self.derivative_wrangler.get_stored_mpole_coefficients_from_full(result)) - def evaluate(self, coeffs, bvec): - taker = self.get_kernel_derivative_taker(bvec) + def evaluate(self, coeffs, bvec, rscale): + taker = self.get_kernel_derivative_taker(bvec, rscale) result = sym.Add(*tuple( - coeff * taker.diff(mi) + coeff + * self.kernel.adjust_proxy_expression( + sympy_vec_subs( + bvec, bvec/rscale, + taker.diff(mi)), + rscale, + sum(mi)) + * rscale**sum(mi) for coeff, mi in zip(coeffs, self.get_coefficient_identifiers()))) - return result - def get_kernel_derivative_taker(self, bvec): - ppkernel = self.kernel.postprocess_at_target( - self.kernel.get_expression(bvec), bvec) + return result - return self.derivative_wrangler.get_derivative_taker(ppkernel, bvec) + def get_kernel_derivative_taker(self, bvec, rscale): + return (self.derivative_wrangler.get_derivative_taker( + self.kernel.get_proxy_expression(bvec), bvec)) def translate_from(self, src_expansion, src_coeff_exprs, dvec): if not isinstance(src_expansion, type(self)): diff --git a/sumpy/kernel.py b/sumpy/kernel.py index aed3d931..f93a7575 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -31,8 +31,10 @@ from sumpy.symbolic import pymbolic_real_norm_2 from pymbolic.primitives import make_sym_vector from pymbolic import var -__doc__ = """ +# FIXME: Rename ProxyExpressionKernel back to ExpressionKernel once we're done. + +__doc__ = """ Kernel interface ---------------- @@ -48,6 +50,7 @@ PDE kernels ----------- .. autoclass:: LaplaceKernel +.. autoclass:: BiharmonicKernel .. autoclass:: HelmholtzKernel .. autoclass:: YukawaKernel .. autoclass:: StokesletKernel @@ -101,12 +104,20 @@ class Kernel(object): .. attribute:: is_complex_valued .. attribute:: dim - *dim* is allowed to be *None* if the dimensionality is not yet - known. + .. automethod:: get_base_kernel + .. automethod:: prepare_loopy_kernel + .. automethod:: get_code_transformer + .. automethod:: get_proxy_expression + .. automethod:: adjust_proxy_expression + .. automethod:: postprocess_at_source + .. automethod:: postprocess_at_target + .. automethod:: get_global_scaling_const + .. automethod:: get_args + .. automethod:: get_source_args """ def __init__(self, dim=None): - self._dim = dim + self.dim = dim # {{{ hashing/pickling/equality @@ -140,21 +151,10 @@ class Kernel(object): # Can't use trivial pickling: hash_value cache must stay unset assert len(self.init_arg_names) == len(state) for name, value in zip(self.init_arg_names, state): - if name == "dim": - name = "_dim" - setattr(self, name, value) # }}} - @property - def dim(self): - if self._dim is None: - raise RuntimeError("the number of dimensions for this kernel " - "has not yet been set") - - return self._dim - def get_base_kernel(self): """Return the kernel being wrapped by this one, or else *self*. @@ -170,16 +170,62 @@ class Kernel(object): def get_code_transformer(self): """Return a function to postprocess the :mod:`pymbolic` expression generated from the result of - :meth:`get_expression` on the way to code generation. + :meth:`get_proxy_expression` on the way to code generation. """ return lambda expr: expr - def get_expression(self, dist_vec): - """Return a :mod:`sympy` expression for the kernel. - - :arg dist_vec: target - source + def get_proxy_expression(self, scaled_dist_vec): + r"""Return a proxy expression :math:`\tilde G(\vec \rho)` where + :math:`\vec \rho` is *scaled_dist_vec*. Here, + :math:`\vec \rho =\vec x/\alpha`, + where :math:`\alpha` is a scaling factor (often written as *rscale*) + used to ensure that :math:`\rho` has roughly unit magnitude and + :math:`\vec x=\text{target} - \text{source}`. (In + an FMM, :math:`\alpha` would be chosen proportional to the + box size.) :math:`\tilde G` should be chosen so that (a) + no powers of :math:`\alpha` appear when symbolic derivatives + of it are taken with respect to :math:`\vec \rho` and + (b) the correct partial derivatives (with input scaling by :math:`\alpha` + can nonetheless be recovered) by :meth:`adjust_proxy_expression`. + + Let :math:`\hat G` be the result of calling + :meth:`adjust_proxy_expression` with + :math:`\partial_{\vec \rho}^\mu \tilde G` and with + *nderivatives* set equal to :math:`|\mu|=\sum_i \mu_i`, + where :math:`\mu` is a multi-index indicating + how many derivatives are being taken along each input axis. + + Then the following equality is required to hold: + + .. math:: + + \hat G(\vec \rho) \alpha^{|\mu|} + = \partial_{\vec \rho}^\mu \tilde G(\vec \rho) \alpha^{|\mu|} + = \partial_{\vec x}^\mu G(\vec x). + + In other words, :meth:`adjust_proxy_expression` is *not* + responsible for applying the scaling :math:`\alpha^{|\mu|}` + resulting from the chain rule (this is managed separately within + :mod:`sumpy` in an effort to control expansion coefficient scaling), + but it should compensate for any other effects of the + scaling of the input variable. + + As a worked example, consider a kernel of :math:`\log(x)`. The + proxy expression would be :math:`\log(\rho)`. + + For ``nderivatives == 0``, :meth:`adjust_proxy_expression` + would return :math:`\tilde G + \log(\alpha)`. + + For ``nderivatives > 0``, :meth:`adjust_proxy_expression` + would return :math:`\tilde G/\alpha^{\text{nderivatives}}`. + + In each case, the added terms merely compensate for the missing + input scaling. + """ + raise NotImplementedError - (Assumes translation invariance of the kernel.) + def adjust_proxy_expression(self, expr, rscale, nderivatives): + r"""See :meth:`get_proxy_expression`. """ raise NotImplementedError @@ -187,6 +233,9 @@ class Kernel(object): """Transform a kernel evaluation or expansion expression in a place where the vector a (something - source) is known. ("something" may be an expansion center or a target.) + + The typical use of this function is to apply source-variable + derivatives to the kernel. """ return expr @@ -194,11 +243,20 @@ class Kernel(object): """Transform a kernel evaluation or expansion expression in a place where the vector b (target - something) is known. ("something" may be an expansion center or a target.) + + The typical use of this function is to apply target-variable + derivatives to the kernel. """ return expr - def get_scaling(self): - """Return a global scaling of the kernel.""" + def get_global_scaling_const(self): + """Return a global scaling constant of the kernel. + Typically, this ensures that the kernel is scaled so that + :math:`\mathcal L(G)(x)=C\delta(x)` with a constant of 1, where + :math:`\mathcal L` is the PDE operator associated with the kernel. + Not to be confused with *rscale*, which keeps expansion + coefficients benignly scaled. + """ raise NotImplementedError def get_args(self): @@ -218,32 +276,42 @@ class Kernel(object): class ExpressionKernel(Kernel): + def get_expression(self, dist_vec): + raise NotImplementedError + + +class ProxyExpressionKernel(Kernel): is_complex_valued = False - init_arg_names = ("dim", "expression", "scaling", "is_complex_valued") + init_arg_names = ("dim", "expression", "kernel_scaling", "is_complex_valued") - def __init__(self, dim, expression, scaling, is_complex_valued): + def __init__(self, dim, proxy_expression, global_scaling_const, + 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. + :arg global_scaling_const: A constant :mod:`pymbolic` expression for the + global scaling of the kernel. Typically, this ensures that + the kernel is scaled so that :math:`\mathcal L(G)(x)=C\delta(x)` + with a constant of 1, where :math:`\mathcal L` is the PDE + operator associated with the kernel. Not to be confused with + *rscale*, which keeps expansion coefficients benignly scaled. """ - # expression and scaling are pymbolic objects because those pickle - # cleanly. D'oh, sympy! + # expression and global_scaling_const are pymbolic objects because + # those pickle cleanly. D'oh, sympy! Kernel.__init__(self, dim) - self.expression = expression - self.scaling = scaling + self.proxy_expression = proxy_expression + self.global_scaling_const = global_scaling_const self.is_complex_valued = is_complex_valued def __getinitargs__(self): - return (self._dim, self.expression, self.scaling, + return (self._dim, self.expression, self.global_scaling_const, self.is_complex_valued) def __repr__(self): @@ -252,36 +320,27 @@ class ExpressionKernel(Kernel): else: return "ExprKnl" - def get_expression(self, dist_vec): - 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)") - + def get_proxy_expression(self, scaled_dist_vec): from sumpy.symbolic import PymbolicToSympyMapperWithSymbols - expr = PymbolicToSympyMapperWithSymbols()(self.expression) + expr = PymbolicToSympyMapperWithSymbols()(self.proxy_expression) - if self.dim != len(dist_vec): + if self.dim != len(scaled_dist_vec): raise ValueError("dist_vec length does not match expected dimension") from sumpy.symbolic import Symbol expr = expr.subs(dict( (Symbol("d%d" % i), dist_vec_i) - for i, dist_vec_i in enumerate(dist_vec) + for i, dist_vec_i in enumerate(scaled_dist_vec) )) return expr - def get_scaling(self): + def get_global_scaling_const(self): """Return a global scaling of the kernel.""" - 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)") - from sumpy.symbolic import PymbolicToSympyMapperWithSymbols - return PymbolicToSympyMapperWithSymbols()(self.scaling) + return PymbolicToSympyMapperWithSymbols()( + self.global_scaling_const) def update_persistent_hash(self, key_hash, key_builder): key_hash.update(type(self).__name__.encode("utf8")) @@ -296,21 +355,21 @@ class ExpressionKernel(Kernel): mapper_method = "map_expression_kernel" -one_kernel_2d = ExpressionKernel( +one_kernel_2d = ProxyExpressionKernel( dim=2, - expression=1, - scaling=1, + proxy_expression=1, + global_scaling_const=1, is_complex_valued=False) -one_kernel_3d = ExpressionKernel( +one_kernel_3d = ProxyExpressionKernel( dim=3, - expression=1, - scaling=1, + proxy_expression=1, + global_scaling_const=1, is_complex_valued=False) # {{{ PDE kernels -class LaplaceKernel(ExpressionKernel): +class LaplaceKernel(ProxyExpressionKernel): init_arg_names = ("dim",) def __init__(self, dim=None): @@ -324,15 +383,28 @@ class LaplaceKernel(ExpressionKernel): expr = 1/r scaling = 1/(4*var("pi")) else: - raise RuntimeError("unsupported dimensionality") + raise NotImplementedError("unsupported dimensionality") - ExpressionKernel.__init__( - self, + super(LaplaceKernel, self).__init__( dim, - expression=expr, - scaling=scaling, + proxy_expression=expr, + global_scaling_const=scaling, is_complex_valued=False) + def adjust_proxy_expression(self, expr, rscale, nderivatives): + if self.dim == 2: + if nderivatives == 0: + import sympy as sp + return expr + sp.log(rscale) + else: + return expr / rscale**nderivatives + + elif self.dim == 3: + return expr / rscale**(nderivatives+1) + + else: + raise NotImplementedError("unsupported dimensionality") + def __getinitargs__(self): return (self._dim,) @@ -354,6 +426,7 @@ class BiharmonicKernel(ExpressionKernel): expr = r**2 * var("log")(r) scaling = 1/(8*var("pi")) elif dim == 3: + raise RuntimeError("unsupported dimensionality") expr = r scaling = 1 # FIXME: Unknown else: @@ -378,7 +451,7 @@ class BiharmonicKernel(ExpressionKernel): mapper_method = "map_biharmonic_kernel" -class HelmholtzKernel(ExpressionKernel): +class HelmholtzKernel(ProxyExpressionKernel): init_arg_names = ("dim", "helmholtz_k_name", "allow_evanescent") def __init__(self, dim=None, helmholtz_k_name="k", @@ -403,16 +476,26 @@ class HelmholtzKernel(ExpressionKernel): else: raise RuntimeError("unsupported dimensionality") - ExpressionKernel.__init__( - self, + super(HelmholtzKernel, self).__init__( dim, - expression=expr, - scaling=scaling, + proxy_expression=expr, + global_scaling_const=scaling, is_complex_valued=True) self.helmholtz_k_name = helmholtz_k_name self.allow_evanescent = allow_evanescent + def adjust_proxy_expression(self, expr, rscale, nderivatives): + import sympy as sp + result = expr.subs( + sp.Symbol(self.helmholtz_k_name), + sp.Symbol(self.helmholtz_k_name) * rscale) + + if self.dim == 2: + return result + else: + return result / rscale**(nderivatives+1) + def __getinitargs__(self): return (self._dim, self.helmholtz_k_name, self.allow_evanescent) @@ -665,8 +748,11 @@ class KernelWrapper(Kernel): def is_complex_valued(self): return self.inner_kernel.is_complex_valued - def get_expression(self, dist_vec): - return self.inner_kernel.get_expression(dist_vec) + def get_proxy_expression(self, scaled_dist_vec): + return self.inner_kernel.get_proxy_expression(scaled_dist_vec) + + def adjust_proxy_expression(self, expr, rscale, nderivatives): + return self.inner_kernel.adjust_proxy_expression(expr, rscale, nderivatives) def postprocess_at_source(self, expr, avec): return self.inner_kernel.postprocess_at_source(expr, avec) @@ -674,8 +760,8 @@ class KernelWrapper(Kernel): def postprocess_at_target(self, expr, avec): return self.inner_kernel.postprocess_at_target(expr, avec) - def get_scaling(self): - return self.inner_kernel.get_scaling() + def get_global_scaling_const(self): + return self.inner_kernel.get_global_scaling_const() def get_code_transformer(self): return self.inner_kernel.get_code_transformer() diff --git a/sumpy/p2e.py b/sumpy/p2e.py index 6a8fcfa9..ce60eea9 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -73,13 +73,16 @@ class P2EBase(KernelCacheWrapper): from sumpy.symbolic import make_sym_vector avec = make_sym_vector("a", self.dim) + import sympy as sp + rscale = sp.Symbol("rscale") + from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() coeff_names = [ sac.assign_unique("coeff%d" % i, coeff_i) for i, coeff_i in enumerate( - self.expansion.coefficients_from_source(avec, None))] + self.expansion.coefficients_from_source(avec, None, rscale))] sac.run_global_cse() @@ -93,7 +96,7 @@ class P2EBase(KernelCacheWrapper): ) def get_cache_key(self): - return (type(self).__name__, self.expansion) + return (type(self).__name__, self.name, self.expansion) # }}} @@ -140,6 +143,7 @@ class P2EFromSingleBox(P2EBase): lp.GlobalArg("box_source_starts,box_source_counts_nonchild", None, shape=None), lp.GlobalArg("centers", None, shape="dim, aligned_nboxes"), + lp.ValueArg("rscale", None), lp.GlobalArg("tgt_expansions", None, shape=("nboxes", ncoeffs), offset=lp.auto), lp.ValueArg("nboxes,aligned_nboxes,tgt_base_ibox", np.int32), @@ -174,10 +178,16 @@ class P2EFromSingleBox(P2EBase): :arg centers: :arg sources: :arg strengths: + :arg rscale: """ + centers = kwargs.pop("centers") + # "1" may be passed for rscale, which won't have its type + # meaningfully inferred. Make the type of rscale explicit. + rscale = centers.dtype.type(kwargs.pop("rscale")) + knl = self.get_cached_optimized_kernel() - return knl(queue, **kwargs) + return knl(queue, centers=centers, rscale=rscale, **kwargs) # }}} diff --git a/sumpy/p2p.py b/sumpy/p2p.py index 9e07fbe0..c7af8280 100644 --- a/sumpy/p2p.py +++ b/sumpy/p2p.py @@ -79,8 +79,12 @@ class P2PBase(KernelComputation, KernelCacheWrapper): sac.assign_unique("knl%d" % i, knl.postprocess_at_target( knl.postprocess_at_source( - knl.get_expression(dvec), dvec), - dvec)) + knl.adjust_proxy_expression( + knl.get_proxy_expression(dvec), + rscale=1, nderivatives=0), + dvec), + dvec) + ) for i, knl in enumerate(self.kernels)] sac.run_global_cse() diff --git a/sumpy/tools.py b/sumpy/tools.py index 2aa4732c..6bda64f6 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -271,7 +271,7 @@ class KernelComputation(object): return [ lp.Assignment(id=None, assignee="knl_%d_scaling" % i, - expression=sympy_conv(kernel.get_scaling()), + expression=sympy_conv(kernel.get_global_scaling_const()), temp_var_type=dtype) for i, (kernel, dtype) in enumerate( zip(self.kernels, self.value_dtypes))] @@ -431,4 +431,14 @@ def my_syntactic_subs(expr, subst_dict): return expr +def sympy_vec_subs(from_vec, to_vec, expr): + subs_pairs = [] + assert (from_vec.rows, from_vec.cols) == (to_vec.rows, to_vec.cols) + for irow in range(from_vec.rows): + for icol in range(from_vec.cols): + subs_pairs.append((from_vec[irow, icol], to_vec[irow, icol])) + + return expr.subs(subs_pairs) + + # vim: fdm=marker diff --git a/sumpy/toys.py b/sumpy/toys.py index e90bf5dc..fae5eaf2 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -127,7 +127,7 @@ class ToyContext(object): # {{{ helpers -def _p2e(psource, center, order, p2e, expn_class, expn_kwargs): +def _p2e(psource, center, rscale, order, p2e, expn_class, expn_kwargs): source_boxes = np.array([0], dtype=np.int32) box_source_starts = np.array([0], dtype=np.int32) box_source_counts_nonchild = np.array( @@ -145,14 +145,15 @@ def _p2e(psource, center, order, p2e, expn_class, expn_kwargs): centers=centers, sources=psource.points, strengths=psource.weights, + rscale=rscale, nboxes=1, tgt_base_ibox=0, #flags="print_hl_cl", out_host=True, **toy_ctx.extra_source_kwargs) - return expn_class(toy_ctx, center, order, coeffs[0], derived_from=psource, - **expn_kwargs) + return expn_class(toy_ctx, center, rscale, order, coeffs[0], + derived_from=psource, **expn_kwargs) def _e2p(psource, targets, e2p): @@ -176,6 +177,7 @@ def _e2p(psource, targets, e2p): box_target_starts=box_target_starts, box_target_counts_nonchild=box_target_counts_nonchild, centers=centers, + rscale=psource.rscale, targets=targets, #flags="print_hl_cl", out_host=True, **toy_ctx.extra_source_kwargs) @@ -183,7 +185,7 @@ def _e2p(psource, targets, e2p): return pot -def _e2e(psource, to_center, to_order, e2e, expn_class, expn_kwargs): +def _e2e(psource, to_center, to_rscale, to_order, e2e, expn_class, expn_kwargs): toy_ctx = psource.toy_ctx target_boxes = np.array([1], dtype=np.int32) @@ -217,7 +219,7 @@ def _e2e(psource, to_center, to_order, e2e, expn_class, expn_kwargs): #flags="print_hl_cl", out_host=True, **toy_ctx.extra_source_kwargs) - return expn_class(toy_ctx, to_center, to_order, to_coeffs[1], + return expn_class(toy_ctx, to_center, to_rscale, to_order, to_coeffs[1], derived_from=psource, **expn_kwargs) # }}} @@ -319,10 +321,11 @@ class ExpansionPotentialSource(PotentialSource): Not used mathematically. Just for visualization, purely advisory. """ - def __init__(self, toy_ctx, center, order, coeffs, derived_from, + def __init__(self, toy_ctx, center, rscale, order, coeffs, derived_from, radius=None, expn_style=None): super(ExpansionPotentialSource, self).__init__(toy_ctx) self.center = np.asarray(center) + self.rscale = rscale self.order = order self.coeffs = coeffs @@ -380,12 +383,12 @@ class Product(PotentialExpressionNode): # }}} -def multipole_expand(psource, center, order=None, **expn_kwargs): +def multipole_expand(psource, center, order=None, rscale=1, **expn_kwargs): if isinstance(psource, PointSources): if order is None: raise ValueError("order may not be None") - return _p2e(psource, center, order, psource.toy_ctx.get_p2m(order), + return _p2e(psource, center, rscale, order, psource.toy_ctx.get_p2m(order), MultipoleExpansion, expn_kwargs) elif isinstance(psource, MultipoleExpansion): @@ -401,12 +404,12 @@ def multipole_expand(psource, center, order=None, **expn_kwargs): % type(psource).__name__) -def local_expand(psource, center, order=None, **expn_kwargs): +def local_expand(psource, center, order=None, rscale=1, **expn_kwargs): if isinstance(psource, PointSources): if order is None: raise ValueError("order may not be None") - return _p2e(psource, center, order, psource.toy_ctx.get_p2l(order), + return _p2e(psource, center, rscale, order, psource.toy_ctx.get_p2l(order), LocalExpansion, expn_kwargs) elif isinstance(psource, MultipoleExpansion): diff --git a/sumpy/version.py b/sumpy/version.py index 6333fc0c..72c4bd46 100644 --- a/sumpy/version.py +++ b/sumpy/version.py @@ -25,4 +25,4 @@ VERSION = (2016, 1) VERSION_STATUS = "beta1" VERSION_TEXT = ".".join(str(x) for x in VERSION) + VERSION_STATUS -KERNEL_VERSION = 15 +KERNEL_VERSION = 16 -- GitLab From 08de50062f9185b3fca7a640c9c21c8cd29f5ab1 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 16:26:05 -0500 Subject: [PATCH 08/60] Remove remaining uses of _dim in sumpy.kernel --- sumpy/kernel.py | 29 ++++++++++------------------- 1 file changed, 10 insertions(+), 19 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index f93a7575..32123dae 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -311,14 +311,11 @@ class ProxyExpressionKernel(Kernel): self.is_complex_valued = is_complex_valued def __getinitargs__(self): - return (self._dim, self.expression, self.global_scaling_const, + return (self.dim, self.expression, self.global_scaling_const, self.is_complex_valued) def __repr__(self): - if self._dim is not None: - return "ExprKnl%dD" % self.dim - else: - return "ExprKnl" + return "ExprKnl%dD" % self.dim def get_proxy_expression(self, scaled_dist_vec): from sumpy.symbolic import PymbolicToSympyMapperWithSymbols @@ -406,13 +403,10 @@ class LaplaceKernel(ProxyExpressionKernel): raise NotImplementedError("unsupported dimensionality") def __getinitargs__(self): - return (self._dim,) + return (self.dim,) def __repr__(self): - if self._dim is not None: - return "LapKnl%dD" % self.dim - else: - return "LapKnl" + return "LapKnl%dD" % self.dim mapper_method = "map_laplace_kernel" @@ -440,13 +434,10 @@ class BiharmonicKernel(ExpressionKernel): is_complex_valued=False) def __getinitargs__(self): - return (self._dim,) + return (self.dim,) def __repr__(self): - if self._dim is not None: - return "BiharmKnl%dD" % self.dim - else: - return "BiharmKnl" + return "BiharmKnl%dD" % self.dim mapper_method = "map_biharmonic_kernel" @@ -497,7 +488,7 @@ class HelmholtzKernel(ProxyExpressionKernel): return result / rscale**(nderivatives+1) def __getinitargs__(self): - return (self._dim, self.helmholtz_k_name, + return (self.dim, self.helmholtz_k_name, self.allow_evanescent) def update_persistent_hash(self, key_hash, key_builder): @@ -563,7 +554,7 @@ class YukawaKernel(ExpressionKernel): self.yukawa_lambda_name = yukawa_lambda_name def __getinitargs__(self): - return (self._dim, self.yukawa_lambda_name) + return (self.dim, self.yukawa_lambda_name) def update_persistent_hash(self, key_hash, key_builder): key_hash.update(type(self).__name__.encode("utf8")) @@ -640,7 +631,7 @@ class StokesletKernel(ExpressionKernel): is_complex_valued=False) def __getinitargs__(self): - return (self._dim, self.icomp, self.jcomp, self.viscosity_mu_name) + return (self.dim, self.icomp, self.jcomp, self.viscosity_mu_name) def update_persistent_hash(self, key_hash, key_builder): key_hash.update(type(self).__name__.encode()) @@ -706,7 +697,7 @@ class StressletKernel(ExpressionKernel): is_complex_valued=False) def __getinitargs__(self): - return (self._dim, self.icomp, self.jcomp, self.kcomp, + return (self.dim, self.icomp, self.jcomp, self.kcomp, self.viscosity_mu_name) def update_persistent_hash(self, key_hash, key_builder): -- GitLab From 7bc070b286f83f952a6be14dab0f398d693dd866 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 16:51:00 -0500 Subject: [PATCH 09/60] Push target-deriv removal through more of the expansion toys --- sumpy/toys.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/sumpy/toys.py b/sumpy/toys.py index fae5eaf2..537d9288 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -49,6 +49,8 @@ class ToyContext(object): self.queue = cl.CommandQueue(self.cl_context) self.kernel = kernel + self.no_target_deriv_kernel = TargetDerivativeRemover()(kernel) + if expansion_factory is None: from sumpy.expansion import DefaultExpansionFactory expansion_factory = DefaultExpansionFactory() @@ -75,52 +77,50 @@ class ToyContext(object): def get_p2m(self, order): from sumpy import P2EFromSingleBox return P2EFromSingleBox(self.cl_context, - self.mpole_expn_class( - TargetDerivativeRemover()(self.kernel), order), + self.mpole_expn_class(self.no_target_deriv_kernel, order), [self.kernel]) @memoize_method def get_p2l(self, order): from sumpy import P2EFromSingleBox return P2EFromSingleBox(self.cl_context, - self.local_expn_class(self.kernel, order), + self.local_expn_class(self.no_target_deriv_kernel, order), [self.kernel]) @memoize_method def get_m2p(self, order): from sumpy import E2PFromSingleBox return E2PFromSingleBox(self.cl_context, - self.mpole_expn_class( - TargetDerivativeRemover()(self.kernel), order), + self.mpole_expn_class(self.no_target_deriv_kernel, order), [self.kernel]) @memoize_method def get_l2p(self, order): from sumpy import E2PFromSingleBox return E2PFromSingleBox(self.cl_context, - self.local_expn_class(self.kernel, order), + self.local_expn_class(self.no_target_deriv_kernel, order), [self.kernel]) @memoize_method def get_m2m(self, from_order, to_order): from sumpy import E2EFromCSR return E2EFromCSR(self.cl_context, - self.mpole_expn_class(self.kernel, from_order), - self.mpole_expn_class(self.kernel, to_order)) + self.mpole_expn_class(self.no_target_deriv_kernel, from_order), + self.mpole_expn_class(self.no_target_deriv_kernel, to_order)) @memoize_method def get_m2l(self, from_order, to_order): from sumpy import E2EFromCSR return E2EFromCSR(self.cl_context, - self.mpole_expn_class(self.kernel, from_order), - self.local_expn_class(self.kernel, to_order)) + self.mpole_expn_class(self.no_target_deriv_kernel, from_order), + self.local_expn_class(self.no_target_deriv_kernel, to_order)) @memoize_method def get_l2l(self, from_order, to_order): from sumpy import E2EFromCSR return E2EFromCSR(self.cl_context, - self.local_expn_class(self.kernel, from_order), - self.local_expn_class(self.kernel, to_order)) + self.local_expn_class(self.no_target_deriv_kernel, from_order), + self.local_expn_class(self.no_target_deriv_kernel, to_order)) # }}} -- GitLab From 250e06a33acfbc77d7a0066ebf6d1633b4236c75 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 16:53:06 -0500 Subject: [PATCH 10/60] Implement rscaling for Taylor local --- examples/coeff-scaling-experiment.py | 13 +++++++------ sumpy/expansion/local.py | 14 +++++++++----- 2 files changed, 16 insertions(+), 11 deletions(-) diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py index 3c14dd45..9cb622ec 100644 --- a/examples/coeff-scaling-experiment.py +++ b/examples/coeff-scaling-experiment.py @@ -9,16 +9,17 @@ import matplotlib.pyplot as plt def main(): dim = 2 + order = 7 from sumpy.kernel import ( # noqa: F401 YukawaKernel, HelmholtzKernel, LaplaceKernel, AxisTargetDerivative) tctx = t.ToyContext( cl.create_some_context(), - #LaplaceKernel(dim), + LaplaceKernel(dim), #AxisTargetDerivative(0, LaplaceKernel(dim)), #YukawaKernel(dim), extra_source_kwargs={"lam": 5}, - HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3}, + #HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3}, ) np.random.seed(12) @@ -31,14 +32,14 @@ def main(): np.ones(50)) mctr = scale*np.array([0., 0, 0])[:dim] - mexp = t.multipole_expand(pt_src, mctr, order=4, rscale=scale) + #mexp = t.multipole_expand(pt_src, mctr, order=order, rscale=scale) lctr = scale*np.array([2.5, 0, 0])[:dim] - #lexp = t.local_expand(mexp, lctr) - print(mexp.coeffs) + lexp = t.local_expand(pt_src, lctr, order=order, rscale=scale) + print(lexp.coeffs) #print(lexp.coeffs) - diff = mexp - pt_src + diff = lexp - pt_src diag = np.sqrt(dim) print(t.l_inf(diff, scale*0.5*diag, center=lctr) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index a39fe99b..1f2d3859 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -113,21 +113,25 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): Coefficients represent derivative values of the kernel. """ - def coefficients_from_source(self, avec, bvec): + def coefficients_from_source(self, avec, bvec, rscale): from sumpy.tools import MiDerivativeTaker ppkernel = self.kernel.postprocess_at_source( - self.kernel.get_expression(avec), avec) + self.kernel.get_proxy_expression(avec), avec) taker = MiDerivativeTaker(ppkernel, avec) - return [taker.diff(mi) for mi in self.get_coefficient_identifiers()] + return [ + self.kernel.adjust_proxy_expression(taker.diff(mi), 1, 0) + * rscale**sum(mi) + for mi in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec): + def evaluate(self, coeffs, bvec, rscale): from sumpy.tools import mi_power, mi_factorial evaluated_coeffs = ( self.derivative_wrangler.get_full_kernel_derivatives_from_stored(coeffs)) + bvec = bvec / rscale result = sum( coeff - * self.kernel.postprocess_at_target(mi_power(bvec, mi), bvec) + * mi_power(bvec, mi) / mi_factorial(mi) for coeff, mi in zip( evaluated_coeffs, self.get_full_coefficient_identifiers())) -- GitLab From 7614dc4163f0a014be90e9ee809871b865038330 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 17:03:36 -0500 Subject: [PATCH 11/60] Fix Fourier-Bessel form/eval local/mpole for rscale --- examples/coeff-scaling-experiment.py | 14 ++++++++------ sumpy/expansion/local.py | 6 ++++-- sumpy/expansion/multipole.py | 14 +++++++++----- 3 files changed, 21 insertions(+), 13 deletions(-) diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py index 9cb622ec..55c37056 100644 --- a/examples/coeff-scaling-experiment.py +++ b/examples/coeff-scaling-experiment.py @@ -16,10 +16,10 @@ def main(): AxisTargetDerivative) tctx = t.ToyContext( cl.create_some_context(), - LaplaceKernel(dim), + #LaplaceKernel(dim), #AxisTargetDerivative(0, LaplaceKernel(dim)), #YukawaKernel(dim), extra_source_kwargs={"lam": 5}, - #HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3}, + HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3}, ) np.random.seed(12) @@ -32,14 +32,16 @@ def main(): np.ones(50)) mctr = scale*np.array([0., 0, 0])[:dim] - #mexp = t.multipole_expand(pt_src, mctr, order=order, rscale=scale) + mexp = t.multipole_expand(pt_src, mctr, order=order, rscale=scale) lctr = scale*np.array([2.5, 0, 0])[:dim] - lexp = t.local_expand(pt_src, lctr, order=order, rscale=scale) - print(lexp.coeffs) + #lexp = t.local_expand(pt_src, lctr, order=order, rscale=scale) + + exp = mexp + print(exp.coeffs) #print(lexp.coeffs) - diff = lexp - pt_src + diff = exp - pt_src diag = np.sqrt(dim) print(t.l_inf(diff, scale*0.5*diag, center=lctr) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 1f2d3859..3507509f 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -216,7 +216,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): def get_coefficient_identifiers(self): return list(range(-self.order, self.order+1)) - def coefficients_from_source(self, avec, bvec): + def coefficients_from_source(self, avec, bvec, rscale): from sumpy.symbolic import sym_real_norm_2 hankel_1 = sym.Function("hankel_1") @@ -227,10 +227,11 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): avec_len = sym_real_norm_2(avec) return [self.kernel.postprocess_at_source( hankel_1(l, arg_scale * avec_len) + * rscale ** abs(l) * sym.exp(sym.I * l * source_angle_rel_center), avec) for l in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec): + def evaluate(self, coeffs, bvec, rscale): from sumpy.symbolic import sym_real_norm_2 bessel_j = sym.Function("bessel_j") bvec_len = sym_real_norm_2(bvec) @@ -241,6 +242,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): return sum(coeffs[self.get_storage_index(l)] * self.kernel.postprocess_at_target( bessel_j(l, arg_scale * bvec_len) + / rscale ** abs(l) * sym.exp(sym.I * l * -target_angle_rel_center), bvec) for l in self.get_coefficient_identifiers()) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 50d18142..20e1cb7d 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -204,7 +204,7 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): def get_coefficient_identifiers(self): return list(range(-self.order, self.order+1)) - def coefficients_from_source(self, avec, bvec): + def coefficients_from_source(self, avec, bvec, rscale): from sumpy.symbolic import sym_real_norm_2 bessel_j = sym.Function("bessel_j") avec_len = sym_real_norm_2(avec) @@ -213,12 +213,15 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): # The coordinates are negated since avec points from source to center. source_angle_rel_center = sym.atan2(-avec[1], -avec[0]) - return [self.kernel.postprocess_at_source( - bessel_j(l, arg_scale * avec_len) * - sym.exp(sym.I * l * -source_angle_rel_center), avec) + return [ + self.kernel.postprocess_at_source( + bessel_j(l, arg_scale * avec_len) + / rscale ** abs(l) + * sym.exp(sym.I * l * -source_angle_rel_center), + avec) for l in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec): + def evaluate(self, coeffs, bvec, rscale): from sumpy.symbolic import sym_real_norm_2 hankel_1 = sym.Function("hankel_1") bvec_len = sym_real_norm_2(bvec) @@ -229,6 +232,7 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): return sum(coeffs[self.get_storage_index(l)] * self.kernel.postprocess_at_target( hankel_1(l, arg_scale * bvec_len) + * rscale ** abs(l) * sym.exp(sym.I * l * target_angle_rel_center), bvec) for l in self.get_coefficient_identifiers()) -- GitLab From dda85e19ce0af59d305e50eb72e3ae035ae9e955 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 18 Aug 2017 17:39:55 -0500 Subject: [PATCH 12/60] Bring remaining kernels back into the fold by simply turning off rscaling --- examples/coeff-scaling-experiment.py | 8 ++- sumpy/expansion/local.py | 6 +++ sumpy/expansion/multipole.py | 6 +++ sumpy/kernel.py | 81 ++++++++++++++++++---------- 4 files changed, 71 insertions(+), 30 deletions(-) diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py index 55c37056..8cb01622 100644 --- a/examples/coeff-scaling-experiment.py +++ b/examples/coeff-scaling-experiment.py @@ -13,13 +13,17 @@ def main(): from sumpy.kernel import ( # noqa: F401 YukawaKernel, HelmholtzKernel, LaplaceKernel, + BiharmonicKernel, StokesletKernel, StressletKernel, AxisTargetDerivative) tctx = t.ToyContext( cl.create_some_context(), - #LaplaceKernel(dim), + LaplaceKernel(dim), #AxisTargetDerivative(0, LaplaceKernel(dim)), #YukawaKernel(dim), extra_source_kwargs={"lam": 5}, - HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3}, + #HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3}, + #BiharmonicKernel(dim), + #StokesletKernel(dim, 1, 1), extra_source_kwargs={"mu": 0.3}, + #StressletKernel(dim, 1, 1, 0), extra_source_kwargs={"mu": 0.3}, ) np.random.seed(12) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 3507509f..f2aec706 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -217,6 +217,9 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): return list(range(-self.order, self.order+1)) def coefficients_from_source(self, avec, bvec, rscale): + if not self.kernel.supports_rscale: + rscale = 1 + from sumpy.symbolic import sym_real_norm_2 hankel_1 = sym.Function("hankel_1") @@ -232,6 +235,9 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): for l in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec, rscale): + if not self.kernel.supports_rscale: + rscale = 1 + from sumpy.symbolic import sym_real_norm_2 bessel_j = sym.Function("bessel_j") bvec_len = sym_real_norm_2(bvec) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 20e1cb7d..ca091397 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -60,6 +60,9 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): from sumpy.tools import mi_power, mi_factorial + if not self.kernel.supports_rscale: + rscale = 1 + avec = avec/rscale if isinstance(kernel, DirectionalSourceDerivative): @@ -94,6 +97,9 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): self.derivative_wrangler.get_stored_mpole_coefficients_from_full(result)) def evaluate(self, coeffs, bvec, rscale): + if not self.kernel.supports_rscale: + rscale = 1 + taker = self.get_kernel_derivative_taker(bvec, rscale) result = sym.Add(*tuple( coeff diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 32123dae..4ae60dbb 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -32,8 +32,6 @@ from pymbolic.primitives import make_sym_vector from pymbolic import var -# FIXME: Rename ProxyExpressionKernel back to ExpressionKernel once we're done. - __doc__ = """ Kernel interface ---------------- @@ -102,6 +100,7 @@ class Kernel(object): """Basic kernel interface. .. attribute:: is_complex_valued + .. attribute:: supports_rscale .. attribute:: dim .. automethod:: get_base_kernel @@ -116,6 +115,8 @@ class Kernel(object): .. automethod:: get_source_args """ + supports_rscale = False + def __init__(self, dim=None): self.dim = dim @@ -276,11 +277,6 @@ class Kernel(object): class ExpressionKernel(Kernel): - def get_expression(self, dist_vec): - raise NotImplementedError - - -class ProxyExpressionKernel(Kernel): is_complex_valued = False init_arg_names = ("dim", "expression", "kernel_scaling", "is_complex_valued") @@ -352,12 +348,12 @@ class ProxyExpressionKernel(Kernel): mapper_method = "map_expression_kernel" -one_kernel_2d = ProxyExpressionKernel( +one_kernel_2d = ExpressionKernel( dim=2, proxy_expression=1, global_scaling_const=1, is_complex_valued=False) -one_kernel_3d = ProxyExpressionKernel( +one_kernel_3d = ExpressionKernel( dim=3, proxy_expression=1, global_scaling_const=1, @@ -366,8 +362,9 @@ one_kernel_3d = ProxyExpressionKernel( # {{{ PDE kernels -class LaplaceKernel(ProxyExpressionKernel): +class LaplaceKernel(ExpressionKernel): init_arg_names = ("dim",) + supports_rscale = True def __init__(self, dim=None): # See (Kress LIE, Thm 6.2) for scaling @@ -413,6 +410,7 @@ class LaplaceKernel(ProxyExpressionKernel): class BiharmonicKernel(ExpressionKernel): init_arg_names = ("dim",) + supports_rscale = False def __init__(self, dim=None): r = pymbolic_real_norm_2(make_sym_vector("d", dim)) @@ -426,13 +424,16 @@ class BiharmonicKernel(ExpressionKernel): else: raise RuntimeError("unsupported dimensionality") - ExpressionKernel.__init__( - self, + super(BiharmonicKernel, self).__init__( dim, - expression=expr, - scaling=scaling, + proxy_expression=expr, + global_scaling_const=scaling, is_complex_valued=False) + def adjust_proxy_expression(self, expr, rscale, nderivatives): + assert rscale == 1 + return expr + def __getinitargs__(self): return (self.dim,) @@ -442,8 +443,9 @@ class BiharmonicKernel(ExpressionKernel): mapper_method = "map_biharmonic_kernel" -class HelmholtzKernel(ProxyExpressionKernel): +class HelmholtzKernel(ExpressionKernel): init_arg_names = ("dim", "helmholtz_k_name", "allow_evanescent") + supports_rscale = True def __init__(self, dim=None, helmholtz_k_name="k", allow_evanescent=False): @@ -484,8 +486,10 @@ class HelmholtzKernel(ProxyExpressionKernel): if self.dim == 2: return result - else: + elif self.dim == 3: return result / rscale**(nderivatives+1) + else: + raise RuntimeError("unsupported dimensionality") def __getinitargs__(self): return (self.dim, self.helmholtz_k_name, @@ -525,6 +529,7 @@ class HelmholtzKernel(ProxyExpressionKernel): class YukawaKernel(ExpressionKernel): init_arg_names = ("dim", "yukawa_lambda_name") + supports_rscale = True def __init__(self, dim=None, yukawa_lambda_name="lam"): """ @@ -544,15 +549,27 @@ class YukawaKernel(ExpressionKernel): else: raise RuntimeError("unsupported dimensionality") - ExpressionKernel.__init__( - self, + super(YukawaKernel, self).__init__( dim, - expression=expr, - scaling=scaling, + proxy_expression=expr, + global_scaling_const=scaling, is_complex_valued=True) self.yukawa_lambda_name = yukawa_lambda_name + def adjust_proxy_expression(self, expr, rscale, nderivatives): + import sympy as sp + result = expr.subs( + sp.Symbol(self.yukawa_lambda_name), + sp.Symbol(self.yukawa_lambda_name) * rscale) + + if self.dim == 2: + return result + elif self.dim == 3: + return result / rscale**(nderivatives+1) + else: + raise RuntimeError("unsupported dimensionality") + def __getinitargs__(self): return (self.dim, self.yukawa_lambda_name) @@ -584,6 +601,7 @@ class YukawaKernel(ExpressionKernel): class StokesletKernel(ExpressionKernel): init_arg_names = ("dim", "icomp", "jcomp", "viscosity_mu_name") + supports_rscale = False def __init__(self, dim, icomp, jcomp, viscosity_mu_name="mu"): r""" @@ -623,13 +641,16 @@ class StokesletKernel(ExpressionKernel): self.icomp = icomp self.jcomp = jcomp - ExpressionKernel.__init__( - self, + super(StokesletKernel, self).__init__( dim, - expression=expr, - scaling=scaling, + proxy_expression=expr, + global_scaling_const=scaling, is_complex_valued=False) + def adjust_proxy_expression(self, expr, rscale, nderivatives): + assert rscale == 1 + return expr + def __getinitargs__(self): return (self.dim, self.icomp, self.jcomp, self.viscosity_mu_name) @@ -652,6 +673,7 @@ class StokesletKernel(ExpressionKernel): class StressletKernel(ExpressionKernel): init_arg_names = ("dim", "icomp", "jcomp", "kcomp", "viscosity_mu_name") + supports_rscale = False def __init__(self, dim=None, icomp=None, jcomp=None, kcomp=None, viscosity_mu_name="mu"): @@ -689,13 +711,16 @@ class StressletKernel(ExpressionKernel): self.jcomp = jcomp self.kcomp = kcomp - ExpressionKernel.__init__( - self, + super(StressletKernel, self).__init__( dim, - expression=expr, - scaling=scaling, + proxy_expression=expr, + global_scaling_const=scaling, is_complex_valued=False) + def adjust_proxy_expression(self, expr, rscale, nderivatives): + assert rscale == 1 + return expr + def __getinitargs__(self): return (self.dim, self.icomp, self.jcomp, self.kcomp, self.viscosity_mu_name) -- GitLab From 5fe535668f157c5079e5caf263421bbe30033c40 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 19 Aug 2017 15:00:29 -0500 Subject: [PATCH 13/60] Track move of pyopencl back to inducer namespace --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 6e23ae95..6690e3d9 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,7 +2,7 @@ numpy sympy==1.0 git+https://github.com/inducer/pymbolic git+https://github.com/inducer/islpy -git+https://github.com/pyopencl/pyopencl +git+https://github.com/inducer/pyopencl git+https://github.com/inducer/boxtree git+https://github.com/inducer/loopy git+https://github.com/inducer/pyfmmlib -- GitLab From ae6fc02fbcaa992f2ecb36e49cd3db7af88a875a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 14:26:31 -0500 Subject: [PATCH 14/60] Remove extraneous 'rscale' argument form get_kernel_derivative_taker --- sumpy/expansion/multipole.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index ca091397..a252a308 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -100,21 +100,20 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): if not self.kernel.supports_rscale: rscale = 1 - taker = self.get_kernel_derivative_taker(bvec, rscale) + taker = self.get_kernel_derivative_taker(bvec) result = sym.Add(*tuple( coeff * self.kernel.adjust_proxy_expression( sympy_vec_subs( bvec, bvec/rscale, taker.diff(mi)), - rscale, - sum(mi)) + rscale, sum(mi)) * rscale**sum(mi) for coeff, mi in zip(coeffs, self.get_coefficient_identifiers()))) return result - def get_kernel_derivative_taker(self, bvec, rscale): + def get_kernel_derivative_taker(self, bvec): return (self.derivative_wrangler.get_derivative_taker( self.kernel.get_proxy_expression(bvec), bvec)) -- GitLab From a6a0c3de067db577bf0abc45ddf2b5b631eff07b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 14:27:06 -0500 Subject: [PATCH 15/60] Add rscale plumbing to e2e --- sumpy/e2e.py | 48 ++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 44 insertions(+), 4 deletions(-) diff --git a/sumpy/e2e.py b/sumpy/e2e.py index 66e5e975..572caa00 100644 --- a/sumpy/e2e.py +++ b/sumpy/e2e.py @@ -82,6 +82,9 @@ class E2EBase(KernelCacheWrapper): src_coeff_exprs = [sym.Symbol("src_coeff%d" % i) for i in range(len(self.src_expansion))] + src_rscale = sym.Symbol("src_rscale") + + tgt_rscale = sym.Symbol("tgt_rscale") from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() @@ -89,7 +92,8 @@ class E2EBase(KernelCacheWrapper): sac.assign_unique("coeff%d" % i, coeff_i) for i, coeff_i in enumerate( self.tgt_expansion.translate_from( - self.src_expansion, src_coeff_exprs, dvec))] + self.src_expansion, src_coeff_exprs, src_rscale, + dvec, tgt_rscale))] sac.run_global_cse() @@ -180,6 +184,7 @@ class E2EFromCSR(E2EBase): """], [ lp.GlobalArg("centers", None, shape="dim, aligned_nboxes"), + lp.ValueArg("src_rscale,tgt_rscale", None), lp.GlobalArg("src_box_starts, src_box_lists", None, shape=None, strides=(1,), offset=lp.auto), lp.ValueArg("aligned_nboxes,tgt_base_ibox,src_base_ibox", @@ -212,11 +217,22 @@ class E2EFromCSR(E2EBase): :arg src_expansions: :arg src_box_starts: :arg src_box_lists: + :arg src_rscale: + :arg tgt_rscale: :arg centers: """ knl = self.get_cached_optimized_kernel() - return knl(queue, **kwargs) + centers = kwargs.pop("centers") + # "1" may be passed for rscale, which won't have its type + # meaningfully inferred. Make the type of rscale explicit. + src_rscale = centers.dtype.type(kwargs.pop("src_rscale")) + tgt_rscale = centers.dtype.type(kwargs.pop("tgt_rscale")) + + return knl(queue, + centers=centers, + src_rscale=src_rscale, tgt_rscale=tgt_rscale, + **kwargs) # }}} @@ -290,6 +306,7 @@ class E2EFromChildren(E2EBase): lp.GlobalArg("target_boxes", None, shape=lp.auto, offset=lp.auto), lp.GlobalArg("centers", None, shape="dim, aligned_nboxes"), + lp.ValueArg("src_rscale,tgt_rscale", None), lp.GlobalArg("box_child_ids", None, shape="nchildren, aligned_nboxes"), lp.GlobalArg("tgt_expansions", None, @@ -321,11 +338,22 @@ class E2EFromChildren(E2EBase): :arg src_expansions: :arg src_box_starts: :arg src_box_lists: + :arg src_rscale: + :arg tgt_rscale: :arg centers: """ knl = self.get_cached_optimized_kernel() - return knl(queue, **kwargs) + centers = kwargs.pop("centers") + # "1" may be passed for rscale, which won't have its type + # meaningfully inferred. Make the type of rscale explicit. + src_rscale = centers.dtype.type(kwargs.pop("src_rscale")) + tgt_rscale = centers.dtype.type(kwargs.pop("tgt_rscale")) + + return knl(queue, + centers=centers, + src_rscale=src_rscale, tgt_rscale=tgt_rscale, + **kwargs) # }}} @@ -386,6 +414,7 @@ class E2EFromParent(E2EBase): lp.GlobalArg("target_boxes", None, shape=lp.auto, offset=lp.auto), lp.GlobalArg("centers", None, shape="dim, naligned_boxes"), + lp.ValueArg("src_rscale,tgt_rscale", None), lp.ValueArg("naligned_boxes,nboxes", np.int32), lp.ValueArg("tgt_base_ibox,src_base_ibox", np.int32), lp.ValueArg("ntgt_level_boxes,nsrc_level_boxes", np.int32), @@ -414,11 +443,22 @@ class E2EFromParent(E2EBase): :arg src_expansions: :arg src_box_starts: :arg src_box_lists: + :arg src_rscale: + :arg tgt_rscale: :arg centers: """ knl = self.get_cached_optimized_kernel() - return knl(queue, **kwargs) + centers = kwargs.pop("centers") + # "1" may be passed for rscale, which won't have its type + # meaningfully inferred. Make the type of rscale explicit. + src_rscale = centers.dtype.type(kwargs.pop("src_rscale")) + tgt_rscale = centers.dtype.type(kwargs.pop("tgt_rscale")) + + return knl(queue, + centers=centers, + src_rscale=src_rscale, tgt_rscale=tgt_rscale, + **kwargs) # }}} -- GitLab From eac1fd522a0d697e835f67a59081757d6fad3a5a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 14:28:12 -0500 Subject: [PATCH 16/60] Add rscale plumbing for toy translations --- sumpy/toys.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/sumpy/toys.py b/sumpy/toys.py index 537d9288..9000da9e 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -216,6 +216,10 @@ def _e2e(psource, to_center, to_rscale, to_order, e2e, expn_class, expn_kwargs): src_box_starts=src_box_starts, src_box_lists=src_box_lists, centers=centers, + + src_rscale=psource.rscale, + tgt_rscale=to_rscale, + #flags="print_hl_cl", out_host=True, **toy_ctx.extra_source_kwargs) @@ -395,7 +399,7 @@ def multipole_expand(psource, center, order=None, rscale=1, **expn_kwargs): if order is None: order = psource.order - return _e2e(psource, center, order, + return _e2e(psource, center, rscale, order, psource.toy_ctx.get_m2m(psource.order, order), MultipoleExpansion, expn_kwargs) @@ -416,7 +420,7 @@ def local_expand(psource, center, order=None, rscale=1, **expn_kwargs): if order is None: order = psource.order - return _e2e(psource, center, order, + return _e2e(psource, center, rscale, order, psource.toy_ctx.get_m2l(psource.order, order), LocalExpansion, expn_kwargs) @@ -424,7 +428,7 @@ def local_expand(psource, center, order=None, rscale=1, **expn_kwargs): if order is None: order = psource.order - return _e2e(psource, center, order, + return _e2e(psource, center, rscale, order, psource.toy_ctx.get_l2l(psource.order, order), LocalExpansion, expn_kwargs) -- GitLab From 3fafe1cf1afd99b4ca6b630b24a0ac2d9003de56 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 14:35:42 -0500 Subject: [PATCH 17/60] Get Taylor M2L rscaling to work --- examples/coeff-scaling-experiment.py | 13 ++++++------- sumpy/expansion/local.py | 23 ++++++++++++++++++++--- 2 files changed, 26 insertions(+), 10 deletions(-) diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py index 8cb01622..c942f08f 100644 --- a/examples/coeff-scaling-experiment.py +++ b/examples/coeff-scaling-experiment.py @@ -9,7 +9,7 @@ import matplotlib.pyplot as plt def main(): dim = 2 - order = 7 + order = 20 from sumpy.kernel import ( # noqa: F401 YukawaKernel, HelmholtzKernel, LaplaceKernel, @@ -36,16 +36,15 @@ def main(): np.ones(50)) mctr = scale*np.array([0., 0, 0])[:dim] - mexp = t.multipole_expand(pt_src, mctr, order=order, rscale=scale) + mexp = t.multipole_expand(pt_src, mctr, order=order, rscale=1) lctr = scale*np.array([2.5, 0, 0])[:dim] - #lexp = t.local_expand(pt_src, lctr, order=order, rscale=scale) + lexp = t.local_expand(mexp, lctr, order=order, rscale=1) - exp = mexp - print(exp.coeffs) - #print(lexp.coeffs) + print(mexp.coeffs) + print(lexp.coeffs) - diff = exp - pt_src + diff = lexp - pt_src diag = np.sqrt(dim) print(t.l_inf(diff, scale*0.5*diag, center=lctr) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index f2aec706..2975203a 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -24,6 +24,7 @@ THE SOFTWARE. from six.moves import range, zip import sumpy.symbolic as sym +from sumpy.tools import sympy_vec_subs from sumpy.expansion import ( ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, @@ -137,13 +138,18 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): evaluated_coeffs, self.get_full_coefficient_identifiers())) return result - def translate_from(self, src_expansion, src_coeff_exprs, dvec): + def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, + dvec, tgt_rscale): logger.info("building translation operator: %s(%d) -> %s(%d): start" % (type(src_expansion).__name__, src_expansion.order, type(self).__name__, self.order)) + if not self.kernel.supports_rscale: + src_rscale = 1 + tgt_rscale = 1 + from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansionBase if isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase): # We know the general form of the multipole expansion is: @@ -155,6 +161,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # # This code speeds up derivative taking by caching all kernel # derivatives. + taker = src_expansion.get_kernel_derivative_taker(dvec) from sumpy.tools import add_mi @@ -165,8 +172,18 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): for coeff, term in zip( src_coeff_exprs, src_expansion.get_coefficient_identifiers()): - kernel_deriv = taker.diff(add_mi(deriv, term)) - local_result.append(coeff * kernel_deriv) + + kernel_deriv = self.kernel.adjust_proxy_expression( + sympy_vec_subs( + dvec, dvec/src_rscale, + taker.diff(add_mi(deriv, term)), + ), + src_rscale, sum(deriv) + sum(term)) + + local_result.append( + coeff * kernel_deriv + * src_rscale**sum(term) + * tgt_rscale**sum(deriv)) result.append(sym.Add(*local_result)) else: from sumpy.tools import MiDerivativeTaker -- GitLab From bae6dcd6a21847046f845e112f3cb0d53e2e88d1 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 14:44:36 -0500 Subject: [PATCH 18/60] Change x -> d, rho -> delta notation in proyx expressions [ci skip] --- sumpy/kernel.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 4ae60dbb..574b3e71 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -176,22 +176,22 @@ class Kernel(object): return lambda expr: expr def get_proxy_expression(self, scaled_dist_vec): - r"""Return a proxy expression :math:`\tilde G(\vec \rho)` where - :math:`\vec \rho` is *scaled_dist_vec*. Here, - :math:`\vec \rho =\vec x/\alpha`, + r"""Return a proxy expression :math:`\tilde G(\vec \delta)` where + :math:`\vec \delta` is *scaled_dist_vec*. Here, + :math:`\vec \delta =\vec d/\alpha`, where :math:`\alpha` is a scaling factor (often written as *rscale*) - used to ensure that :math:`\rho` has roughly unit magnitude and - :math:`\vec x=\text{target} - \text{source}`. (In + used to ensure that :math:`\delta` has roughly unit magnitude and + :math:`\vec d=\text{target} - \text{source}`. (In an FMM, :math:`\alpha` would be chosen proportional to the box size.) :math:`\tilde G` should be chosen so that (a) no powers of :math:`\alpha` appear when symbolic derivatives - of it are taken with respect to :math:`\vec \rho` and + of it are taken with respect to :math:`\vec \delta` and (b) the correct partial derivatives (with input scaling by :math:`\alpha` can nonetheless be recovered) by :meth:`adjust_proxy_expression`. Let :math:`\hat G` be the result of calling :meth:`adjust_proxy_expression` with - :math:`\partial_{\vec \rho}^\mu \tilde G` and with + :math:`\partial_{\vec \delta}^\mu \tilde G` and with *nderivatives* set equal to :math:`|\mu|=\sum_i \mu_i`, where :math:`\mu` is a multi-index indicating how many derivatives are being taken along each input axis. @@ -200,9 +200,9 @@ class Kernel(object): .. math:: - \hat G(\vec \rho) \alpha^{|\mu|} - = \partial_{\vec \rho}^\mu \tilde G(\vec \rho) \alpha^{|\mu|} - = \partial_{\vec x}^\mu G(\vec x). + \hat G(\vec \delta) \alpha^{|\mu|} + = \partial_{\vec \delta}^\mu \tilde G(\vec \delta) \alpha^{|\mu|} + = \partial_{\vec d}^\mu G(\vec d). In other words, :meth:`adjust_proxy_expression` is *not* responsible for applying the scaling :math:`\alpha^{|\mu|}` @@ -212,7 +212,7 @@ class Kernel(object): scaling of the input variable. As a worked example, consider a kernel of :math:`\log(x)`. The - proxy expression would be :math:`\log(\rho)`. + proxy expression would be :math:`\log(\delta)`. For ``nderivatives == 0``, :meth:`adjust_proxy_expression` would return :math:`\tilde G + \log(\alpha)`. -- GitLab From ca09dd95f5837751a69de035e5b74983834df674 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 16:15:43 -0500 Subject: [PATCH 19/60] Add 'factor' to adjust_proxy_expression to tweak symbolic order of operations for 'maximum simplification opportunity' --- examples/coeff-scaling-experiment.py | 4 +-- sumpy/expansion/local.py | 9 +++--- sumpy/expansion/multipole.py | 3 +- sumpy/kernel.py | 44 ++++++++++++++++------------ sumpy/p2p.py | 2 +- 5 files changed, 34 insertions(+), 28 deletions(-) diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py index c942f08f..5d72d911 100644 --- a/examples/coeff-scaling-experiment.py +++ b/examples/coeff-scaling-experiment.py @@ -36,10 +36,10 @@ def main(): np.ones(50)) mctr = scale*np.array([0., 0, 0])[:dim] - mexp = t.multipole_expand(pt_src, mctr, order=order, rscale=1) + mexp = t.multipole_expand(pt_src, mctr, order=order, rscale=scale) lctr = scale*np.array([2.5, 0, 0])[:dim] - lexp = t.local_expand(mexp, lctr, order=order, rscale=1) + lexp = t.local_expand(mexp, lctr, order=order, rscale=scale) print(mexp.coeffs) print(lexp.coeffs) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 2975203a..8b5d84bd 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -121,7 +121,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): taker = MiDerivativeTaker(ppkernel, avec) return [ - self.kernel.adjust_proxy_expression(taker.diff(mi), 1, 0) + self.kernel.adjust_proxy_expression(taker.diff(mi), 1, 0, factor=1) * rscale**sum(mi) for mi in self.get_coefficient_identifiers()] @@ -178,12 +178,11 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): dvec, dvec/src_rscale, taker.diff(add_mi(deriv, term)), ), - src_rscale, sum(deriv) + sum(term)) + src_rscale, sum(deriv) + sum(term), + factor=src_rscale**sum(term)) local_result.append( - coeff * kernel_deriv - * src_rscale**sum(term) - * tgt_rscale**sum(deriv)) + coeff * kernel_deriv * tgt_rscale**sum(deriv)) result.append(sym.Add(*local_result)) else: from sumpy.tools import MiDerivativeTaker diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index a252a308..1047a8a0 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -107,8 +107,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): sympy_vec_subs( bvec, bvec/rscale, taker.diff(mi)), - rscale, sum(mi)) - * rscale**sum(mi) + rscale, sum(mi), factor=rscale**sum(mi)) for coeff, mi in zip(coeffs, self.get_coefficient_identifiers()))) return result diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 574b3e71..474efc91 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -225,8 +225,15 @@ class Kernel(object): """ raise NotImplementedError - def adjust_proxy_expression(self, expr, rscale, nderivatives): + def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): r"""See :meth:`get_proxy_expression`. + + :arg factor: a factor by which the result is to be mutliplied. + The idea here is that :meth:`adjust_proxy_expression` will + often result in the application of a scale factor. Supplying + an additional scale factors gives :mod:`sympy` the opportunity + to first simplify the scale factor before applying it to the + result. """ raise NotImplementedError @@ -385,13 +392,13 @@ class LaplaceKernel(ExpressionKernel): global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives): + def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): if self.dim == 2: if nderivatives == 0: import sympy as sp - return expr + sp.log(rscale) + return (expr + sp.log(rscale)) * factor else: - return expr / rscale**nderivatives + return expr * (factor / rscale**nderivatives) elif self.dim == 3: return expr / rscale**(nderivatives+1) @@ -430,9 +437,9 @@ class BiharmonicKernel(ExpressionKernel): global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives): + def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): assert rscale == 1 - return expr + return expr * factor def __getinitargs__(self): return (self.dim,) @@ -478,16 +485,16 @@ class HelmholtzKernel(ExpressionKernel): self.helmholtz_k_name = helmholtz_k_name self.allow_evanescent = allow_evanescent - def adjust_proxy_expression(self, expr, rscale, nderivatives): + def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): import sympy as sp result = expr.subs( sp.Symbol(self.helmholtz_k_name), sp.Symbol(self.helmholtz_k_name) * rscale) if self.dim == 2: - return result + return result * factor elif self.dim == 3: - return result / rscale**(nderivatives+1) + return result * (factor / rscale**(nderivatives+1)) else: raise RuntimeError("unsupported dimensionality") @@ -557,16 +564,16 @@ class YukawaKernel(ExpressionKernel): self.yukawa_lambda_name = yukawa_lambda_name - def adjust_proxy_expression(self, expr, rscale, nderivatives): + def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): import sympy as sp result = expr.subs( sp.Symbol(self.yukawa_lambda_name), sp.Symbol(self.yukawa_lambda_name) * rscale) if self.dim == 2: - return result + return result * factor elif self.dim == 3: - return result / rscale**(nderivatives+1) + return result * (factor / rscale**(nderivatives+1)) else: raise RuntimeError("unsupported dimensionality") @@ -647,9 +654,9 @@ class StokesletKernel(ExpressionKernel): global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives): + def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): assert rscale == 1 - return expr + return expr * factor def __getinitargs__(self): return (self.dim, self.icomp, self.jcomp, self.viscosity_mu_name) @@ -717,9 +724,9 @@ class StressletKernel(ExpressionKernel): global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives): + def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): assert rscale == 1 - return expr + return expr * factor def __getinitargs__(self): return (self.dim, self.icomp, self.jcomp, self.kcomp, @@ -767,8 +774,9 @@ class KernelWrapper(Kernel): def get_proxy_expression(self, scaled_dist_vec): return self.inner_kernel.get_proxy_expression(scaled_dist_vec) - def adjust_proxy_expression(self, expr, rscale, nderivatives): - return self.inner_kernel.adjust_proxy_expression(expr, rscale, nderivatives) + def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): + return self.inner_kernel.adjust_proxy_expression( + expr, rscale, nderivatives, factor) def postprocess_at_source(self, expr, avec): return self.inner_kernel.postprocess_at_source(expr, avec) diff --git a/sumpy/p2p.py b/sumpy/p2p.py index c7af8280..cb94934d 100644 --- a/sumpy/p2p.py +++ b/sumpy/p2p.py @@ -81,7 +81,7 @@ class P2PBase(KernelComputation, KernelCacheWrapper): knl.postprocess_at_source( knl.adjust_proxy_expression( knl.get_proxy_expression(dvec), - rscale=1, nderivatives=0), + rscale=1, nderivatives=0, factor=1), dvec), dvec) ) -- GitLab From 85a11b88333ab78fcfe432a79d9edbaad70762b4 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 16:33:59 -0500 Subject: [PATCH 20/60] Rscaling: Get Taylor L2L to work --- examples/coeff-scaling-experiment.py | 14 ++++++++------ sumpy/expansion/local.py | 7 +++++-- 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py index 5d72d911..9e30f7f4 100644 --- a/examples/coeff-scaling-experiment.py +++ b/examples/coeff-scaling-experiment.py @@ -9,7 +9,7 @@ import matplotlib.pyplot as plt def main(): dim = 2 - order = 20 + order = 10 from sumpy.kernel import ( # noqa: F401 YukawaKernel, HelmholtzKernel, LaplaceKernel, @@ -35,16 +35,18 @@ def main(): scale * pts, np.ones(50)) - mctr = scale*np.array([0., 0, 0])[:dim] - mexp = t.multipole_expand(pt_src, mctr, order=order, rscale=scale) + #mctr = scale*np.array([0., 0, 0])[:dim] + #mexp = t.multipole_expand(pt_src, mctr, order=order, rscale=scale) + lctr1 = scale*np.array([2.8, 0, 0])[:dim] lctr = scale*np.array([2.5, 0, 0])[:dim] - lexp = t.local_expand(mexp, lctr, order=order, rscale=scale) + lexp1 = t.local_expand(pt_src, lctr1, order=order, rscale=scale) + lexp = t.local_expand(lexp1, lctr, order=order, rscale=scale) - print(mexp.coeffs) + #print(mexp.coeffs) print(lexp.coeffs) - diff = lexp - pt_src + diff = lexp1 - pt_src diag = np.sqrt(dim) print(t.l_inf(diff, scale*0.5*diag, center=lctr) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 8b5d84bd..a22b0605 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -186,9 +186,12 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): result.append(sym.Add(*local_result)) else: from sumpy.tools import MiDerivativeTaker - expr = src_expansion.evaluate(src_coeff_exprs, dvec) + expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=1) taker = MiDerivativeTaker(expr, dvec) - result = [taker.diff(mi) for mi in self.get_coefficient_identifiers()] + rscale = tgt_rscale/src_rscale + result = [ + taker.diff(mi) * rscale**sum(mi) + for mi in self.get_coefficient_identifiers()] logger.info("building translation operator: done") return result -- GitLab From 940496106621721fd108bbe6764894742051f625 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 16:51:33 -0500 Subject: [PATCH 21/60] Fix Taylor L2L [ci skip] --- examples/coeff-scaling-experiment.py | 6 ++++-- sumpy/expansion/local.py | 5 ++--- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py index 9e30f7f4..3cf1cf65 100644 --- a/examples/coeff-scaling-experiment.py +++ b/examples/coeff-scaling-experiment.py @@ -36,7 +36,9 @@ def main(): np.ones(50)) #mctr = scale*np.array([0., 0, 0])[:dim] - #mexp = t.multipole_expand(pt_src, mctr, order=order, rscale=scale) + #mctr1 = scale*np.array([0., 0.2, 0])[:dim] + #mexp1 = t.multipole_expand(pt_src, mctr1, order=order, rscale=scale) + #mexp = t.multipole_expand(mexp1, mctr, order=order, rscale=scale) lctr1 = scale*np.array([2.8, 0, 0])[:dim] lctr = scale*np.array([2.5, 0, 0])[:dim] @@ -46,7 +48,7 @@ def main(): #print(mexp.coeffs) print(lexp.coeffs) - diff = lexp1 - pt_src + diff = lexp - pt_src diag = np.sqrt(dim) print(t.l_inf(diff, scale*0.5*diag, center=lctr) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index a22b0605..0d37ef27 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -186,11 +186,10 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): result.append(sym.Add(*local_result)) else: from sumpy.tools import MiDerivativeTaker - expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=1) + expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=src_rscale) taker = MiDerivativeTaker(expr, dvec) - rscale = tgt_rscale/src_rscale result = [ - taker.diff(mi) * rscale**sum(mi) + (taker.diff(mi) * tgt_rscale**sum(mi)).expand() for mi in self.get_coefficient_identifiers()] logger.info("building translation operator: done") -- GitLab From 88a6765f46e35d4ecb0b9857d9351896262741d3 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 17:06:22 -0500 Subject: [PATCH 22/60] Rscaling: get Taylor M2M to work [ci skip] --- examples/coeff-scaling-experiment.py | 20 ++++++++++---------- sumpy/expansion/multipole.py | 12 ++++++++++-- 2 files changed, 20 insertions(+), 12 deletions(-) diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py index 3cf1cf65..9dc4a6d5 100644 --- a/examples/coeff-scaling-experiment.py +++ b/examples/coeff-scaling-experiment.py @@ -35,20 +35,20 @@ def main(): scale * pts, np.ones(50)) - #mctr = scale*np.array([0., 0, 0])[:dim] - #mctr1 = scale*np.array([0., 0.2, 0])[:dim] - #mexp1 = t.multipole_expand(pt_src, mctr1, order=order, rscale=scale) - #mexp = t.multipole_expand(mexp1, mctr, order=order, rscale=scale) + mctr = scale*np.array([0., 0, 0])[:dim] + mctr1 = scale*np.array([0.2, 0, 0])[:dim] + mexp1 = t.multipole_expand(pt_src, mctr1, order=order, rscale=scale) + mexp = t.multipole_expand(mexp1, mctr, order=order, rscale=scale) - lctr1 = scale*np.array([2.8, 0, 0])[:dim] + #lctr1 = scale*np.array([2.8, 0, 0])[:dim] lctr = scale*np.array([2.5, 0, 0])[:dim] - lexp1 = t.local_expand(pt_src, lctr1, order=order, rscale=scale) - lexp = t.local_expand(lexp1, lctr, order=order, rscale=scale) + #lexp1 = t.local_expand(pt_src, lctr1, order=order, rscale=scale) + #lexp = t.local_expand(lexp1, lctr, order=order, rscale=scale) - #print(mexp.coeffs) - print(lexp.coeffs) + print(mexp.coeffs) + #print(lexp.coeffs) - diff = lexp - pt_src + diff = mexp - pt_src diag = np.sqrt(dim) print(t.l_inf(diff, scale*0.5*diag, center=lctr) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 1047a8a0..8fcfc645 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -116,12 +116,17 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): return (self.derivative_wrangler.get_derivative_taker( self.kernel.get_proxy_expression(bvec), bvec)) - def translate_from(self, src_expansion, src_coeff_exprs, dvec): + def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, + dvec, tgt_rscale): if not isinstance(src_expansion, type(self)): raise RuntimeError("do not know how to translate %s to " "Taylor multipole expansion" % type(src_expansion).__name__) + if not self.kernel.supports_rscale: + src_rscale = 1 + tgt_rscale = 1 + logger.info("building translation operator: %s(%d) -> %s(%d): start" % (type(src_expansion).__name__, src_expansion.order, @@ -161,7 +166,10 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): contrib *= (binomial(n, k) * dvec[idim]**(n-k)) - result[i] += contrib + result[i] += ( + contrib + * src_rscale**sum(src_mi) + / tgt_rscale**sum(tgt_mi)) result[i] /= mi_factorial(tgt_mi) -- GitLab From dcb8aeea368e7353895d244f3370ad04e8458197 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 21:34:32 -0500 Subject: [PATCH 23/60] Rscale L2L: Add comment regarding expand/order of ops [ci skip] --- sumpy/expansion/local.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 0d37ef27..9bbc76d1 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -188,6 +188,11 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): from sumpy.tools import MiDerivativeTaker expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=src_rscale) taker = MiDerivativeTaker(expr, dvec) + + # Rscale/operand magnitude is fairly sensitive to the order of + # operations--which is something we don't have fantastic control + # over at the symbolic level. This expand moves the two canceling + # "rscales" closer to each other in the hope of helping with that. result = [ (taker.diff(mi) * tgt_rscale**sum(mi)).expand() for mi in self.get_coefficient_identifiers()] -- GitLab From 5f70ddd84633676f450438b0c7320cac04220315 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 22:13:58 -0500 Subject: [PATCH 24/60] Rscaling: Add support for Fourier-Bessel translations [ci skip] --- examples/coeff-scaling-experiment.py | 14 +++++++------- sumpy/expansion/local.py | 11 +++++++++-- sumpy/expansion/multipole.py | 5 ++++- 3 files changed, 20 insertions(+), 10 deletions(-) diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py index 9dc4a6d5..a9bbf32f 100644 --- a/examples/coeff-scaling-experiment.py +++ b/examples/coeff-scaling-experiment.py @@ -38,17 +38,17 @@ def main(): mctr = scale*np.array([0., 0, 0])[:dim] mctr1 = scale*np.array([0.2, 0, 0])[:dim] mexp1 = t.multipole_expand(pt_src, mctr1, order=order, rscale=scale) - mexp = t.multipole_expand(mexp1, mctr, order=order, rscale=scale) + mexp = t.multipole_expand(mexp1, mctr, order=order, rscale=2*scale) - #lctr1 = scale*np.array([2.8, 0, 0])[:dim] + lctr1 = scale*np.array([2.8, 0, 0])[:dim] lctr = scale*np.array([2.5, 0, 0])[:dim] - #lexp1 = t.local_expand(pt_src, lctr1, order=order, rscale=scale) - #lexp = t.local_expand(lexp1, lctr, order=order, rscale=scale) + lexp1 = t.local_expand(mexp, lctr1, order=order, rscale=scale) + lexp = t.local_expand(lexp1, lctr, order=order, rscale=2*scale) - print(mexp.coeffs) - #print(lexp.coeffs) + #print(mexp.coeffs) + print(lexp.coeffs) - diff = mexp - pt_src + diff = lexp - pt_src diag = np.sqrt(dim) print(t.l_inf(diff, scale*0.5*diag, center=lctr) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 9bbc76d1..18593de8 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -275,7 +275,8 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): * sym.exp(sym.I * l * -target_angle_rel_center), bvec) for l in self.get_coefficient_identifiers()) - def translate_from(self, src_expansion, src_coeff_exprs, dvec): + def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, + dvec, tgt_rscale): from sumpy.symbolic import sym_real_norm_2 arg_scale = self.get_bessel_arg_scaling() @@ -289,6 +290,8 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): translated_coeffs.append( sum(src_coeff_exprs[src_expansion.get_storage_index(m)] * bessel_j(m - l, arg_scale * dvec_len) + / src_rscale ** abs(m) + * tgt_rscale ** abs(l) * sym.exp(sym.I * (m - l) * -new_center_angle_rel_old_center) for m in src_expansion.get_coefficient_identifiers())) return translated_coeffs @@ -300,7 +303,11 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): translated_coeffs = [] for l in self.get_coefficient_identifiers(): translated_coeffs.append( - sum((-1) ** l * hankel_1(m + l, arg_scale * dvec_len) + sum( + (-1) ** l + * hankel_1(m + l, arg_scale * dvec_len) + * src_rscale ** abs(m) + * tgt_rscale ** abs(l) * sym.exp(sym.I * (m + l) * new_center_angle_rel_old_center) * src_coeff_exprs[src_expansion.get_storage_index(m)] for m in src_expansion.get_coefficient_identifiers())) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 8fcfc645..286a99a0 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -248,7 +248,8 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): * sym.exp(sym.I * l * target_angle_rel_center), bvec) for l in self.get_coefficient_identifiers()) - def translate_from(self, src_expansion, src_coeff_exprs, dvec): + def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, + dvec, tgt_rscale): if not isinstance(src_expansion, type(self)): raise RuntimeError("do not know how to translate %s to %s" % (type(src_expansion).__name__, @@ -265,6 +266,8 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): translated_coeffs.append( sum(src_coeff_exprs[src_expansion.get_storage_index(m)] * bessel_j(m - l, arg_scale * dvec_len) + * src_rscale ** abs(m) + / tgt_rscale ** abs(l) * sym.exp(sym.I * (m - l) * new_center_angle_rel_old_center) for m in src_expansion.get_coefficient_identifiers())) return translated_coeffs -- GitLab From aba94058d1550756a24358df216d4d34e3f91a0d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 20 Aug 2017 23:23:46 -0500 Subject: [PATCH 25/60] Fix indentation [ci skip] --- sumpy/expansion/local.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 18593de8..9a4de3d8 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -310,7 +310,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): * tgt_rscale ** abs(l) * sym.exp(sym.I * (m + l) * new_center_angle_rel_old_center) * src_coeff_exprs[src_expansion.get_storage_index(m)] - for m in src_expansion.get_coefficient_identifiers())) + for m in src_expansion.get_coefficient_identifiers())) return translated_coeffs raise RuntimeError("do not know how to translate %s to %s" -- GitLab From 612e9f5297c597734e6094ad5c1229417c36250a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 23 Aug 2017 17:07:25 -0500 Subject: [PATCH 26/60] Adapt expansion toys to be able to deal with source kernel kwargs for source derivatives [ci skip] --- sumpy/toys.py | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/sumpy/toys.py b/sumpy/toys.py index 9000da9e..251c0803 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -44,7 +44,8 @@ class ToyContext(object): mpole_expn_class=None, local_expn_class=None, expansion_factory=None, - extra_source_kwargs=None): + extra_source_kwargs=None, + extra_kernel_kwargs=None): self.cl_context = cl_context self.queue = cl.CommandQueue(self.cl_context) self.kernel = kernel @@ -61,12 +62,20 @@ class ToyContext(object): local_expn_class = \ expansion_factory.get_local_expansion_class(kernel) + self.mpole_expn_class = mpole_expn_class + self.local_expn_class = local_expn_class + if extra_source_kwargs is None: extra_source_kwargs = {} + if extra_kernel_kwargs is None: + extra_kernel_kwargs = {} - self.mpole_expn_class = mpole_expn_class - self.local_expn_class = local_expn_class self.extra_source_kwargs = extra_source_kwargs + self.extra_kernel_kwargs = extra_kernel_kwargs + + extra_source_and_kernel_kwargs = extra_source_kwargs.copy() + extra_source_and_kernel_kwargs.update(extra_kernel_kwargs) + self.extra_source_and_kernel_kwargs = extra_source_and_kernel_kwargs @memoize_method def get_p2p(self): @@ -150,7 +159,8 @@ def _p2e(psource, center, rscale, order, p2e, expn_class, expn_kwargs): tgt_base_ibox=0, #flags="print_hl_cl", - out_host=True, **toy_ctx.extra_source_kwargs) + out_host=True, + **toy_ctx.extra_source_and_kernel_kwargs) return expn_class(toy_ctx, center, rscale, order, coeffs[0], derived_from=psource, **expn_kwargs) @@ -180,7 +190,7 @@ def _e2p(psource, targets, e2p): rscale=psource.rscale, targets=targets, #flags="print_hl_cl", - out_host=True, **toy_ctx.extra_source_kwargs) + out_host=True, **toy_ctx.extra_kernel_kwargs) return pot @@ -221,7 +231,7 @@ def _e2e(psource, to_center, to_rscale, to_order, e2e, expn_class, expn_kwargs): tgt_rscale=to_rscale, #flags="print_hl_cl", - out_host=True, **toy_ctx.extra_source_kwargs) + out_host=True, **toy_ctx.extra_kernel_kwargs) return expn_class(toy_ctx, to_center, to_rscale, to_order, to_coeffs[1], derived_from=psource, **expn_kwargs) @@ -307,7 +317,9 @@ class PointSources(PotentialSource): def eval(self, targets): evt, (potential,) = self.toy_ctx.get_p2p()( self.toy_ctx.queue, targets, self.points, [self.weights], - out_host=True, **self.toy_ctx.extra_source_kwargs) + out_host=True, + **self.toy_ctx.extra_source_and_kernel_kwargs, + ) return potential -- GitLab From e1681877ba660438d5c181bcc86351671ffcbe6b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 23 Aug 2017 19:02:35 -0500 Subject: [PATCH 27/60] Add some measure of explanation of rscaling --- sumpy/kernel.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 474efc91..ba04f84c 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -175,6 +175,44 @@ class Kernel(object): """ return lambda expr: expr + # How does rscaling work for symbolic expressions? + # ------------------------------------------------ + # + # The thing is that all this is just a big no-op. We are trying to + # differentiate + # + # $G(\frac{d}\alpha \alpha)$ + # + # while shuffling around powers of $\alpha$ to keep all terms as close as + # possible to unit-scale. + # + # We call $\delta = d/\alpha$ the unit-scaled geometry variable. + # + # We could just differentiate + # + # $ \delta \mapsto G( \delta \alpha),$ + # + # which would in principle tell us everything we need to know about the + # scaling behavior of the kernel. In this iteration I chose not to do that + # because I felt the extra $\alpha$ would put an unnecessary "clutter" + # burden on the symbolic machinery, making it scale even less well than it + # already does. I don't know whether that's true and/or conceptually + # mistaken. + # + # Instead, I decided that we'll differentiate "just $G$" and have a + # fix-up routine that provides the result of substituting $\alpha \delta$ + # into $\partial^\mu G$ after the fact by modifying the expression "from + # the outside". We then supply $\delta$ as input. The thought behind this + # is that any arithmetic in $\partial^\mu G$ is carried out on (close-to) + # unit-scale quantities. + # + # But since we've differentiated a $G$ "with the wrong scaling", there is + # also a chain-rule term $\alpha^{|\mu|}$ that needs to be multiplied in. + # This is deferred until the evaluation of the expansion. + # + # See also the MR that introduced this: + # https://gitlab.tiker.net/inducer/sumpy/merge_requests/46 + def get_proxy_expression(self, scaled_dist_vec): r"""Return a proxy expression :math:`\tilde G(\vec \delta)` where :math:`\vec \delta` is *scaled_dist_vec*. Here, -- GitLab From 9c83549af40d3c232f2f6e6170cd62bad7a15e77 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 23 Aug 2017 19:03:05 -0500 Subject: [PATCH 28/60] Adapt test_kernels to rscale --- test/test_kernels.py | 36 +++++++++++++++++++++++++++++++----- 1 file changed, 31 insertions(+), 5 deletions(-) diff --git a/test/test_kernels.py b/test/test_kernels.py index c32435d3..65295e29 100644 --- a/test/test_kernels.py +++ b/test/test_kernels.py @@ -214,6 +214,8 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative): targets = fp.points + rscale = 0.5 # pick something non-1 + # {{{ apply p2e evt, (mpoles,) = p2e(queue, @@ -225,6 +227,7 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative): strengths=strengths, nboxes=1, tgt_base_ibox=0, + rscale=rscale, #flags="print_hl_cl", out_host=True, **extra_source_kwargs) @@ -247,6 +250,8 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative): box_target_counts_nonchild=box_target_counts_nonchild, centers=centers, targets=targets, + rscale=rscale, + #flags="print_hl_cl", out_host=True, **extra_kwargs) @@ -402,7 +407,7 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class): orders = [2, 3, 4] nboxes = centers.shape[-1] - def eval_at(e2p, source_box_nr): + def eval_at(e2p, source_box_nr, rscale): e2p_target_boxes = np.array([source_box_nr], dtype=np.int32) # These are indexed by global box numbers. @@ -422,6 +427,9 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class): box_target_counts_nonchild=e2p_box_target_counts_nonchild, centers=centers, targets=targets, + + rscale=rscale, + out_host=True, **extra_kwargs ) @@ -452,6 +460,11 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class): # }}} + m1_rscale = 0.5 + m2_rscale = 0.25 + l1_rscale = 0.5 + l2_rscale = 0.25 + # {{{ apply P2M p2m_source_boxes = np.array([0], dtype=np.int32) @@ -469,6 +482,7 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class): sources=sources, strengths=strengths, nboxes=nboxes, + rscale=m1_rscale, tgt_base_ibox=0, @@ -479,7 +493,7 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class): ntargets = targets.shape[-1] - pot = eval_at(m2p, 0) + pot = eval_at(m2p, 0, m1_rscale) err = la.norm((pot - pot_direct)/res**2) err = err / (la.norm(pot_direct) / res**2) @@ -503,12 +517,16 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class): src_box_starts=m2m_src_box_starts, src_box_lists=m2m_src_box_lists, centers=centers, + + src_rscale=m1_rscale, + tgt_rscale=m2_rscale, + #flags="print_hl_cl", out_host=True, **extra_kwargs) # }}} - pot = eval_at(m2p, 1) + pot = eval_at(m2p, 1, m2_rscale) err = la.norm((pot - pot_direct)/res**2) err = err / (la.norm(pot_direct) / res**2) @@ -531,12 +549,16 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class): src_box_starts=m2l_src_box_starts, src_box_lists=m2l_src_box_lists, centers=centers, + + src_rscale=m2_rscale, + tgt_rscale=l1_rscale, + #flags="print_hl_cl", out_host=True, **extra_kwargs) # }}} - pot = eval_at(l2p, 2) + pot = eval_at(l2p, 2, l1_rscale) err = la.norm((pot - pot_direct)/res**2) err = err / (la.norm(pot_direct) / res**2) @@ -559,12 +581,16 @@ def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class): src_box_starts=l2l_src_box_starts, src_box_lists=l2l_src_box_lists, centers=centers, + + src_rscale=l1_rscale, + tgt_rscale=l2_rscale, + #flags="print_hl_wrapper", out_host=True, **extra_kwargs) # }}} - pot = eval_at(l2p, 3) + pot = eval_at(l2p, 3, l2_rscale) err = la.norm((pot - pot_direct)/res**2) err = err / (la.norm(pot_direct) / res**2) -- GitLab From 2fa84960dd6a4e2a0322e974c02f1aceb97cf7b0 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 10:59:32 -0500 Subject: [PATCH 29/60] Fix 2D Helmholtz rscaling --- sumpy/kernel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index ba04f84c..23200773 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -530,7 +530,7 @@ class HelmholtzKernel(ExpressionKernel): sp.Symbol(self.helmholtz_k_name) * rscale) if self.dim == 2: - return result * factor + return result * (factor / rscale**nderivatives) elif self.dim == 3: return result * (factor / rscale**(nderivatives+1)) else: -- GitLab From a94b04e8e8848c2b91292c5677414119ff4bea93 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 12:28:12 -0500 Subject: [PATCH 30/60] Make formation/evaluation of rscale conforming expansions work. M2L doesn't yet work for those --- sumpy/expansion/__init__.py | 77 +++++++++++++++++++++++++----------- sumpy/expansion/local.py | 16 +++++--- sumpy/expansion/multipole.py | 8 ++-- sumpy/tools.py | 26 ++++++++++-- 4 files changed, 91 insertions(+), 36 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 38107482..ef174099 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -23,11 +23,15 @@ THE SOFTWARE. """ from six.moves import range +import six import numpy as np import logging from pytools import memoize_method import sumpy.symbolic as sym from sumpy.tools import MiDerivativeTaker +import sympy as sp +from collections import defaultdict + __doc__ = """ .. autoclass:: ExpansionBase @@ -164,7 +168,8 @@ class FullDerivativeWrangler(DerivativeWrangler): get_coefficient_identifiers = ( DerivativeWrangler.get_full_coefficient_identifiers) - def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives): + def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives, + rscale): return stored_kernel_derivatives get_stored_mpole_coefficients_from_full = ( @@ -207,28 +212,37 @@ def _spmv(spmat, x, sparse_vectors): class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): - - def __init__(self, order, dim): - DerivativeWrangler.__init__(self, order, dim) + _rscale_symbol = sp.Symbol("_sumpy_rscale_placeholder") def get_coefficient_identifiers(self): return self.stored_identifiers - def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives): - return _spmv(self.coeff_matrix, stored_kernel_derivatives, + def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives, + rscale): + coeff_matrix = self.get_coefficient_matrix(rscale) + return _spmv(coeff_matrix, stored_kernel_derivatives, sparse_vectors=False) - def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients): + def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients, + rscale): # = M^T x, where M = coeff matrix + + coeff_matrix = self.get_coefficient_matrix(rscale) result = [0] * len(self.stored_identifiers) for row, coeff in enumerate(full_mpole_coefficients): - for col, val in self.coeff_matrix[row]: + for col, val in coeff_matrix[row]: result[col] += coeff * val return result - def precompute_recurrences(self): + @property + def stored_identifiers(self): + stored_identifiers, coeff_matrix = self._get_stored_ids_and_coeff_mat() + return stored_identifiers + + @memoize_method + def get_coefficient_matrix(self, rscale): """ - Build up a matrix that expresses every derivative in terms of a + Return a matrix that expresses every derivative in terms of a set of "stored" derivatives. For example, for the recurrence @@ -245,6 +259,20 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): ^ rows = one for every derivative """ + stored_identifiers, coeff_matrix = self._get_stored_ids_and_coeff_mat() + + # substitute actual rscale for internal placeholder + return defaultdict(lambda: [], + ((irow, [ + (icol, + coeff.subs(self._rscale_symbol, rscale) + if isinstance(coeff, sp.Basic) + else coeff) + for icol, coeff in row]) + for irow, row in six.iteritems(coeff_matrix))) + + @memoize_method + def _get_stored_ids_and_coeff_mat(self): stored_identifiers = [] identifiers_so_far = {} @@ -253,7 +281,6 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): logger.debug("computing recurrence for Taylor coefficients: start") # Sparse matrix, indexed by row - from collections import defaultdict coeff_matrix_transpose = defaultdict(lambda: []) # Build up the matrix transpose by row. @@ -275,7 +302,7 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): identifiers_so_far[identifier] = i - self.stored_identifiers = stored_identifiers + stored_identifiers = stored_identifiers coeff_matrix = defaultdict(lambda: []) for i, row in iteritems(coeff_matrix_transpose): @@ -288,11 +315,18 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): logger.debug("number of Taylor coefficients was reduced from {orig} to {red}" .format(orig=len(self.get_full_coefficient_identifiers()), - red=len(self.stored_identifiers))) + red=len(stored_identifiers))) - self.coeff_matrix = coeff_matrix + return stored_identifiers, coeff_matrix def try_get_recurrence_for_derivative(self, deriv, in_terms_of): + """ + :arg deriv: a tuple of integers identifying a derivative for which + a recurrence is sought + :arg in_terms_of: a container supporting efficient containment testing + indicating availability of derivatives for use in recurrences + """ + raise NotImplementedError def get_derivative_taker(self, expr, var_list): @@ -302,14 +336,10 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): - def __init__(self, order, dim): - LinearRecurrenceBasedDerivativeWrangler.__init__(self, order, dim) - self.precompute_recurrences() - def try_get_recurrence_for_derivative(self, deriv, in_terms_of): deriv = np.array(deriv, dtype=int) - for dim in np.nonzero(2 <= deriv)[0]: + for dim in np.where(2 <= deriv)[0]: # Check if we can reduce this dimension in terms of the other # dimensions. @@ -335,14 +365,15 @@ class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): def __init__(self, order, dim, helmholtz_k_name): - LinearRecurrenceBasedDerivativeWrangler.__init__(self, order, dim) + super(HelmholtzDerivativeWrangler, self).__init__(order, dim) self.helmholtz_k_name = helmholtz_k_name - self.precompute_recurrences() def try_get_recurrence_for_derivative(self, deriv, in_terms_of): deriv = np.array(deriv, dtype=int) - for dim in np.nonzero(2 <= deriv)[0]: + rscale = self._rscale_symbol + + for dim in np.where(2 <= deriv)[0]: # Check if we can reduce this dimension in terms of the other # dimensions. @@ -363,7 +394,7 @@ class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): coeffs[needed_deriv] = -1 else: k = sym.Symbol(self.helmholtz_k_name) - coeffs[tuple(reduced_deriv)] = -k*k + coeffs[tuple(reduced_deriv)] = -k*k*rscale*rscale return coeffs # }}} diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 9a4de3d8..c72c6053 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -121,14 +121,17 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): taker = MiDerivativeTaker(ppkernel, avec) return [ - self.kernel.adjust_proxy_expression(taker.diff(mi), 1, 0, factor=1) + self.kernel.adjust_proxy_expression( + taker.diff(mi, rscale), + 1, 0, factor=1) * rscale**sum(mi) for mi in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec, rscale): from sumpy.tools import mi_power, mi_factorial evaluated_coeffs = ( - self.derivative_wrangler.get_full_kernel_derivatives_from_stored(coeffs)) + self.derivative_wrangler.get_full_kernel_derivatives_from_stored( + coeffs, rscale)) bvec = bvec / rscale result = sum( coeff @@ -176,7 +179,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): kernel_deriv = self.kernel.adjust_proxy_expression( sympy_vec_subs( dvec, dvec/src_rscale, - taker.diff(add_mi(deriv, term)), + taker.diff(add_mi(deriv, term), src_rscale), ), src_rscale, sum(deriv) + sum(term), factor=src_rscale**sum(term)) @@ -191,10 +194,11 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # Rscale/operand magnitude is fairly sensitive to the order of # operations--which is something we don't have fantastic control - # over at the symbolic level. This expand moves the two canceling - # "rscales" closer to each other in the hope of helping with that. + # over at the symbolic level. The '.expand()' below moves the two + # canceling "rscales" closer to each other in the hope of helping + # with that. result = [ - (taker.diff(mi) * tgt_rscale**sum(mi)).expand() + (taker.diff(mi, src_rscale) * tgt_rscale**sum(mi)).expand() for mi in self.get_coefficient_identifiers()] logger.info("building translation operator: done") diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 286a99a0..68fd416f 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -94,7 +94,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): mi_power(avec, mi) / mi_factorial(mi) for mi in self.get_full_coefficient_identifiers()] return ( - self.derivative_wrangler.get_stored_mpole_coefficients_from_full(result)) + self.derivative_wrangler.get_stored_mpole_coefficients_from_full( + result, rscale)) def evaluate(self, coeffs, bvec, rscale): if not self.kernel.supports_rscale: @@ -106,7 +107,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): * self.kernel.adjust_proxy_expression( sympy_vec_subs( bvec, bvec/rscale, - taker.diff(mi)), + taker.diff(mi, rscale)), rscale, sum(mi), factor=rscale**sum(mi)) for coeff, mi in zip(coeffs, self.get_coefficient_identifiers()))) @@ -175,7 +176,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): logger.info("building translation operator: done") return ( - self.derivative_wrangler.get_stored_mpole_coefficients_from_full(result)) + self.derivative_wrangler.get_stored_mpole_coefficients_from_full( + result, tgt_rscale)) class VolumeTaylorMultipoleExpansion( diff --git a/sumpy/tools.py b/sumpy/tools.py index 6bda64f6..5b1293ae 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -64,7 +64,11 @@ class MiDerivativeTaker(object): def mi_dist(self, a, b): return np.array(a, dtype=int) - np.array(b, dtype=int) - def diff(self, mi): + def diff(self, mi, rscale): + """ + :arg rscale: Only used if needed for recurrences. All actual coefficient + radius-scaling is handled outside of this routine. + """ try: expr = self.cache_by_mi[mi] except KeyError: @@ -95,14 +99,25 @@ class MiDerivativeTaker(object): class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker): """ - See :class:`sumpy.expansion.LinearRecurrenceBasedDerivativeWrangler` + The derivative taker for expansions that use + :class:`sumpy.expansion.LinearRecurrenceBasedDerivativeWrangler` """ def __init__(self, expr, var_list, wrangler): - MiDerivativeTaker.__init__(self, expr, var_list) + super(LinearRecurrenceBasedMiDerivativeTaker, self).__init__( + expr, var_list) self.wrangler = wrangler - def diff(self, mi): + @memoize_method + def diff(self, mi, rscale): + """ + :arg mi: a multi-index (tuple) indicating how many x/y derivatives are + to be taken. + :arg rscales: *rscale* is used if the derivative taking uses + recurrences across different derivative orders (which may have + different scaling). All actual coefficient radius-scaling is handled + outside of this routine. + """ try: expr = self.cache_by_mi[mi] except KeyError: @@ -129,6 +144,9 @@ class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker): expr = expr.diff(next_deriv) self.cache_by_mi[next_mi] = expr + + expr = expr.subs(self.wrangler._rscale_symbol, rscale) + return expr # }}} -- GitLab From 8aada7a34e89b3a3a5543495180380943c11b2cd Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 13:51:33 -0500 Subject: [PATCH 31/60] Get Helmholtz-conforming rscaled M2L to work --- sumpy/tools.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/tools.py b/sumpy/tools.py index 5b1293ae..3f06f3f1 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -145,7 +145,7 @@ class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker): self.cache_by_mi[next_mi] = expr - expr = expr.subs(self.wrangler._rscale_symbol, rscale) + expr = expr.subs(self.wrangler._rscale_symbol, 1) return expr -- GitLab From 5c66fcd8b79f6c396d52cebd941b5f7fcb9b3853 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 15:00:01 -0500 Subject: [PATCH 32/60] Remove nonsense rscale argument from MiDerivativeTaker.diff() --- sumpy/expansion/local.py | 6 +++--- sumpy/expansion/multipole.py | 2 +- sumpy/tools.py | 12 ++---------- 3 files changed, 6 insertions(+), 14 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index c72c6053..871717a7 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -122,7 +122,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): taker = MiDerivativeTaker(ppkernel, avec) return [ self.kernel.adjust_proxy_expression( - taker.diff(mi, rscale), + taker.diff(mi), 1, 0, factor=1) * rscale**sum(mi) for mi in self.get_coefficient_identifiers()] @@ -179,7 +179,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): kernel_deriv = self.kernel.adjust_proxy_expression( sympy_vec_subs( dvec, dvec/src_rscale, - taker.diff(add_mi(deriv, term), src_rscale), + taker.diff(add_mi(deriv, term)), ), src_rscale, sum(deriv) + sum(term), factor=src_rscale**sum(term)) @@ -198,7 +198,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # canceling "rscales" closer to each other in the hope of helping # with that. result = [ - (taker.diff(mi, src_rscale) * tgt_rscale**sum(mi)).expand() + (taker.diff(mi) * tgt_rscale**sum(mi)).expand() for mi in self.get_coefficient_identifiers()] logger.info("building translation operator: done") diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 68fd416f..1c0bb414 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -107,7 +107,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): * self.kernel.adjust_proxy_expression( sympy_vec_subs( bvec, bvec/rscale, - taker.diff(mi, rscale)), + taker.diff(mi)), rscale, sum(mi), factor=rscale**sum(mi)) for coeff, mi in zip(coeffs, self.get_coefficient_identifiers()))) diff --git a/sumpy/tools.py b/sumpy/tools.py index 3f06f3f1..64a509f3 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -64,11 +64,7 @@ class MiDerivativeTaker(object): def mi_dist(self, a, b): return np.array(a, dtype=int) - np.array(b, dtype=int) - def diff(self, mi, rscale): - """ - :arg rscale: Only used if needed for recurrences. All actual coefficient - radius-scaling is handled outside of this routine. - """ + def diff(self, mi): try: expr = self.cache_by_mi[mi] except KeyError: @@ -109,14 +105,10 @@ class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker): self.wrangler = wrangler @memoize_method - def diff(self, mi, rscale): + def diff(self, mi): """ :arg mi: a multi-index (tuple) indicating how many x/y derivatives are to be taken. - :arg rscales: *rscale* is used if the derivative taking uses - recurrences across different derivative orders (which may have - different scaling). All actual coefficient radius-scaling is handled - outside of this routine. """ try: expr = self.cache_by_mi[mi] -- GitLab From 9a0d1c27b9e36b6e064ed0bc9d19050201824e80 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 16:30:03 -0500 Subject: [PATCH 33/60] Move 1/rscale**nderivatives out of adjust_proxy_expression, remove kernel_support_rscale, add use_rscale flag, simplify. --- sumpy/expansion/__init__.py | 17 ++++--- sumpy/expansion/local.py | 21 ++++---- sumpy/expansion/multipole.py | 40 +++++++++------ sumpy/kernel.py | 98 ++++++++++++++++-------------------- sumpy/p2p.py | 4 +- sumpy/tools.py | 4 +- 6 files changed, 91 insertions(+), 93 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index ef174099..f8549244 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -51,7 +51,7 @@ logger = logging.getLogger(__name__) class ExpansionBase(object): - def __init__(self, kernel, order): + def __init__(self, kernel, order, use_rscale=None): # Don't be tempted to remove target derivatives here. # Line Taylor QBX can't do without them, because it can't # apply those derivatives to the expanded quantity. @@ -59,6 +59,11 @@ class ExpansionBase(object): self.kernel = kernel self.order = order + if use_rscale is None: + use_rscale = True + + self.use_rscale = use_rscale + # {{{ propagate kernel interface @property @@ -287,7 +292,7 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): from six import iteritems for i, identifier in enumerate(self.get_full_coefficient_identifiers()): expr = self.try_get_recurrence_for_derivative( - identifier, identifiers_so_far) + identifier, identifiers_so_far, rscale=self._rscale_symbol) if expr is None: # Identifier should be stored @@ -319,7 +324,7 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): return stored_identifiers, coeff_matrix - def try_get_recurrence_for_derivative(self, deriv, in_terms_of): + def try_get_recurrence_for_derivative(self, deriv, in_terms_of, rscale): """ :arg deriv: a tuple of integers identifying a derivative for which a recurrence is sought @@ -336,7 +341,7 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): - def try_get_recurrence_for_derivative(self, deriv, in_terms_of): + def try_get_recurrence_for_derivative(self, deriv, in_terms_of, rscale): deriv = np.array(deriv, dtype=int) for dim in np.where(2 <= deriv)[0]: @@ -368,11 +373,9 @@ class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): super(HelmholtzDerivativeWrangler, self).__init__(order, dim) self.helmholtz_k_name = helmholtz_k_name - def try_get_recurrence_for_derivative(self, deriv, in_terms_of): + def try_get_recurrence_for_derivative(self, deriv, in_terms_of, rscale): deriv = np.array(deriv, dtype=int) - rscale = self._rscale_symbol - for dim in np.where(2 <= deriv)[0]: # Check if we can reduce this dimension in terms of the other # dimensions. diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 871717a7..d911ab4f 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -117,14 +117,11 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): def coefficients_from_source(self, avec, bvec, rscale): from sumpy.tools import MiDerivativeTaker ppkernel = self.kernel.postprocess_at_source( - self.kernel.get_proxy_expression(avec), avec) + self.kernel.get_expression(avec), avec) taker = MiDerivativeTaker(ppkernel, avec) return [ - self.kernel.adjust_proxy_expression( - taker.diff(mi), - 1, 0, factor=1) - * rscale**sum(mi) + taker.diff(mi) * rscale ** sum(mi) for mi in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec, rscale): @@ -149,7 +146,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): type(self).__name__, self.order)) - if not self.kernel.supports_rscale: + if not self.use_rscale: src_rscale = 1 tgt_rscale = 1 @@ -181,8 +178,8 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): dvec, dvec/src_rscale, taker.diff(add_mi(deriv, term)), ), - src_rscale, sum(deriv) + sum(term), - factor=src_rscale**sum(term)) + src_rscale, sum(deriv) + sum(term) + ) / src_rscale**sum(deriv) local_result.append( coeff * kernel_deriv * tgt_rscale**sum(deriv)) @@ -244,7 +241,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): return list(range(-self.order, self.order+1)) def coefficients_from_source(self, avec, bvec, rscale): - if not self.kernel.supports_rscale: + if not self.use_rscale: rscale = 1 from sumpy.symbolic import sym_real_norm_2 @@ -262,7 +259,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): for l in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec, rscale): - if not self.kernel.supports_rscale: + if not self.use_rscale: rscale = 1 from sumpy.symbolic import sym_real_norm_2 @@ -283,6 +280,10 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): dvec, tgt_rscale): from sumpy.symbolic import sym_real_norm_2 + if not self.use_rscale: + src_rscale = 1 + tgt_rscale = 1 + arg_scale = self.get_bessel_arg_scaling() if isinstance(src_expansion, type(self)): diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 1c0bb414..20176f8b 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -60,11 +60,9 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): from sumpy.tools import mi_power, mi_factorial - if not self.kernel.supports_rscale: + if not self.use_rscale: rscale = 1 - avec = avec/rscale - if isinstance(kernel, DirectionalSourceDerivative): if kernel.get_base_kernel() is not kernel.inner_kernel: raise NotImplementedError("more than one source derivative " @@ -88,8 +86,10 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): - mi_power(avec, derivative_mi) * mi[idim] * dir_vec[idim]) for i, mi in enumerate(coeff_identifiers): - result[i] /= mi_factorial(mi) + result[i] /= (mi_factorial(mi) * rscale ** sum(mi)) else: + avec = avec/rscale + result = [ mi_power(avec, mi) / mi_factorial(mi) for mi in self.get_full_coefficient_identifiers()] @@ -98,7 +98,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): result, rscale)) def evaluate(self, coeffs, bvec, rscale): - if not self.kernel.supports_rscale: + if not self.use_rscale: rscale = 1 taker = self.get_kernel_derivative_taker(bvec) @@ -108,14 +108,14 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): sympy_vec_subs( bvec, bvec/rscale, taker.diff(mi)), - rscale, sum(mi), factor=rscale**sum(mi)) + rscale, sum(mi)) for coeff, mi in zip(coeffs, self.get_coefficient_identifiers()))) return result def get_kernel_derivative_taker(self, bvec): return (self.derivative_wrangler.get_derivative_taker( - self.kernel.get_proxy_expression(bvec), bvec)) + self.kernel.get_expression(bvec), bvec)) def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, dvec, tgt_rscale): @@ -124,7 +124,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): "Taylor multipole expansion" % type(src_expansion).__name__) - if not self.kernel.supports_rscale: + if not self.use_rscale: src_rscale = 1 tgt_rscale = 1 @@ -169,8 +169,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): result[i] += ( contrib - * src_rscale**sum(src_mi) - / tgt_rscale**sum(tgt_mi)) + * (src_rscale**sum(src_mi) / tgt_rscale**sum(tgt_mi))) result[i] /= mi_factorial(tgt_mi) @@ -219,6 +218,9 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): return list(range(-self.order, self.order+1)) def coefficients_from_source(self, avec, bvec, rscale): + if not self.use_rscale: + rscale = 1 + from sumpy.symbolic import sym_real_norm_2 bessel_j = sym.Function("bessel_j") avec_len = sym_real_norm_2(avec) @@ -236,6 +238,9 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): for l in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec, rscale): + if not self.use_rscale: + rscale = 1 + from sumpy.symbolic import sym_real_norm_2 hankel_1 = sym.Function("hankel_1") bvec_len = sym_real_norm_2(bvec) @@ -256,6 +261,11 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): raise RuntimeError("do not know how to translate %s to %s" % (type(src_expansion).__name__, type(self).__name__)) + + if not self.use_rscale: + src_rscale = 1 + tgt_rscale = 1 + from sumpy.symbolic import sym_real_norm_2 dvec_len = sym_real_norm_2(dvec) bessel_j = sym.Function("bessel_j") @@ -276,24 +286,26 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): class H2DMultipoleExpansion(_HankelBased2DMultipoleExpansion): - def __init__(self, kernel, order): + def __init__(self, kernel, order, use_rscale=None): from sumpy.kernel import HelmholtzKernel assert (isinstance(kernel.get_base_kernel(), HelmholtzKernel) and kernel.dim == 2) - super(H2DMultipoleExpansion, self).__init__(kernel, order) + super(H2DMultipoleExpansion, self).__init__( + kernel, order, use_rscale=use_rscale) def get_bessel_arg_scaling(self): return sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) class Y2DMultipoleExpansion(_HankelBased2DMultipoleExpansion): - def __init__(self, kernel, order): + def __init__(self, kernel, order, use_rscale=None): from sumpy.kernel import YukawaKernel assert (isinstance(kernel.get_base_kernel(), YukawaKernel) and kernel.dim == 2) - super(Y2DMultipoleExpansion, self).__init__(kernel, order) + super(Y2DMultipoleExpansion, self).__init__( + kernel, order, use_rscale=use_rscale) def get_bessel_arg_scaling(self): return sym.I * sym.Symbol(self.kernel.get_base_kernel().yukawa_lambda_name) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 23200773..f77e43f7 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -100,13 +100,12 @@ class Kernel(object): """Basic kernel interface. .. attribute:: is_complex_valued - .. attribute:: supports_rscale .. attribute:: dim .. automethod:: get_base_kernel .. automethod:: prepare_loopy_kernel .. automethod:: get_code_transformer - .. automethod:: get_proxy_expression + .. automethod:: get_expression .. automethod:: adjust_proxy_expression .. automethod:: postprocess_at_source .. automethod:: postprocess_at_target @@ -115,8 +114,6 @@ class Kernel(object): .. automethod:: get_source_args """ - supports_rscale = False - def __init__(self, dim=None): self.dim = dim @@ -171,7 +168,7 @@ class Kernel(object): def get_code_transformer(self): """Return a function to postprocess the :mod:`pymbolic` expression generated from the result of - :meth:`get_proxy_expression` on the way to code generation. + :meth:`get_expression` on the way to code generation. """ return lambda expr: expr @@ -213,7 +210,7 @@ class Kernel(object): # See also the MR that introduced this: # https://gitlab.tiker.net/inducer/sumpy/merge_requests/46 - def get_proxy_expression(self, scaled_dist_vec): + def get_expression(self, scaled_dist_vec): r"""Return a proxy expression :math:`\tilde G(\vec \delta)` where :math:`\vec \delta` is *scaled_dist_vec*. Here, :math:`\vec \delta =\vec d/\alpha`, @@ -263,16 +260,10 @@ class Kernel(object): """ raise NotImplementedError - def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): - r"""See :meth:`get_proxy_expression`. - - :arg factor: a factor by which the result is to be mutliplied. - The idea here is that :meth:`adjust_proxy_expression` will - often result in the application of a scale factor. Supplying - an additional scale factors gives :mod:`sympy` the opportunity - to first simplify the scale factor before applying it to the - result. + def adjust_proxy_expression(self, expr, rscale, nderivatives): + r"""See :meth:`get_expression`. """ + raise NotImplementedError def postprocess_at_source(self, expr, avec): @@ -324,9 +315,10 @@ class Kernel(object): class ExpressionKernel(Kernel): is_complex_valued = False - init_arg_names = ("dim", "expression", "kernel_scaling", "is_complex_valued") + init_arg_names = ("dim", "expression", "global_scaling_const", + "is_complex_valued") - def __init__(self, dim, proxy_expression, global_scaling_const, + def __init__(self, dim, expression, global_scaling_const, is_complex_valued): """ :arg expression: A :mod:`pymbolic` expression depending on @@ -347,7 +339,7 @@ class ExpressionKernel(Kernel): Kernel.__init__(self, dim) - self.proxy_expression = proxy_expression + self.expression = expression self.global_scaling_const = global_scaling_const self.is_complex_valued = is_complex_valued @@ -358,9 +350,9 @@ class ExpressionKernel(Kernel): def __repr__(self): return "ExprKnl%dD" % self.dim - def get_proxy_expression(self, scaled_dist_vec): + def get_expression(self, scaled_dist_vec): from sumpy.symbolic import PymbolicToSympyMapperWithSymbols - expr = PymbolicToSympyMapperWithSymbols()(self.proxy_expression) + expr = PymbolicToSympyMapperWithSymbols()(self.expression) if self.dim != len(scaled_dist_vec): raise ValueError("dist_vec length does not match expected dimension") @@ -395,12 +387,12 @@ class ExpressionKernel(Kernel): one_kernel_2d = ExpressionKernel( dim=2, - proxy_expression=1, + expression=1, global_scaling_const=1, is_complex_valued=False) one_kernel_3d = ExpressionKernel( dim=3, - proxy_expression=1, + expression=1, global_scaling_const=1, is_complex_valued=False) @@ -409,7 +401,6 @@ one_kernel_3d = ExpressionKernel( class LaplaceKernel(ExpressionKernel): init_arg_names = ("dim",) - supports_rscale = True def __init__(self, dim=None): # See (Kress LIE, Thm 6.2) for scaling @@ -426,20 +417,20 @@ class LaplaceKernel(ExpressionKernel): super(LaplaceKernel, self).__init__( dim, - proxy_expression=expr, + expression=expr, global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): + def adjust_proxy_expression(self, expr, rscale, nderivatives): if self.dim == 2: if nderivatives == 0: import sympy as sp - return (expr + sp.log(rscale)) * factor + return (expr + sp.log(rscale)) else: - return expr * (factor / rscale**nderivatives) + return expr elif self.dim == 3: - return expr / rscale**(nderivatives+1) + return expr / rscale else: raise NotImplementedError("unsupported dimensionality") @@ -455,7 +446,6 @@ class LaplaceKernel(ExpressionKernel): class BiharmonicKernel(ExpressionKernel): init_arg_names = ("dim",) - supports_rscale = False def __init__(self, dim=None): r = pymbolic_real_norm_2(make_sym_vector("d", dim)) @@ -471,13 +461,13 @@ class BiharmonicKernel(ExpressionKernel): super(BiharmonicKernel, self).__init__( dim, - proxy_expression=expr, + expression=expr, global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): + def adjust_proxy_expression(self, expr, rscale, nderivatives): assert rscale == 1 - return expr * factor + return expr def __getinitargs__(self): return (self.dim,) @@ -490,7 +480,6 @@ class BiharmonicKernel(ExpressionKernel): class HelmholtzKernel(ExpressionKernel): init_arg_names = ("dim", "helmholtz_k_name", "allow_evanescent") - supports_rscale = True def __init__(self, dim=None, helmholtz_k_name="k", allow_evanescent=False): @@ -516,23 +505,23 @@ class HelmholtzKernel(ExpressionKernel): super(HelmholtzKernel, self).__init__( dim, - proxy_expression=expr, + expression=expr, global_scaling_const=scaling, is_complex_valued=True) self.helmholtz_k_name = helmholtz_k_name self.allow_evanescent = allow_evanescent - def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): + def adjust_proxy_expression(self, expr, rscale, nderivatives): import sympy as sp result = expr.subs( sp.Symbol(self.helmholtz_k_name), sp.Symbol(self.helmholtz_k_name) * rscale) if self.dim == 2: - return result * (factor / rscale**nderivatives) + return result elif self.dim == 3: - return result * (factor / rscale**(nderivatives+1)) + return result / rscale else: raise RuntimeError("unsupported dimensionality") @@ -574,7 +563,6 @@ class HelmholtzKernel(ExpressionKernel): class YukawaKernel(ExpressionKernel): init_arg_names = ("dim", "yukawa_lambda_name") - supports_rscale = True def __init__(self, dim=None, yukawa_lambda_name="lam"): """ @@ -596,22 +584,22 @@ class YukawaKernel(ExpressionKernel): super(YukawaKernel, self).__init__( dim, - proxy_expression=expr, + expression=expr, global_scaling_const=scaling, is_complex_valued=True) self.yukawa_lambda_name = yukawa_lambda_name - def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): + def adjust_proxy_expression(self, expr, rscale, nderivatives): import sympy as sp result = expr.subs( - sp.Symbol(self.yukawa_lambda_name), - sp.Symbol(self.yukawa_lambda_name) * rscale) + sp.Symbol(self.helmholtz_k_name), + sp.Symbol(self.helmholtz_k_name) * rscale) if self.dim == 2: - return result * factor + return result elif self.dim == 3: - return result * (factor / rscale**(nderivatives+1)) + return result / rscale else: raise RuntimeError("unsupported dimensionality") @@ -646,7 +634,6 @@ class YukawaKernel(ExpressionKernel): class StokesletKernel(ExpressionKernel): init_arg_names = ("dim", "icomp", "jcomp", "viscosity_mu_name") - supports_rscale = False def __init__(self, dim, icomp, jcomp, viscosity_mu_name="mu"): r""" @@ -688,13 +675,13 @@ class StokesletKernel(ExpressionKernel): super(StokesletKernel, self).__init__( dim, - proxy_expression=expr, + expression=expr, global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): + def adjust_proxy_expression(self, expr, rscale, nderivatives): assert rscale == 1 - return expr * factor + return expr def __getinitargs__(self): return (self.dim, self.icomp, self.jcomp, self.viscosity_mu_name) @@ -718,7 +705,6 @@ class StokesletKernel(ExpressionKernel): class StressletKernel(ExpressionKernel): init_arg_names = ("dim", "icomp", "jcomp", "kcomp", "viscosity_mu_name") - supports_rscale = False def __init__(self, dim=None, icomp=None, jcomp=None, kcomp=None, viscosity_mu_name="mu"): @@ -758,13 +744,13 @@ class StressletKernel(ExpressionKernel): super(StressletKernel, self).__init__( dim, - proxy_expression=expr, + expression=expr, global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): + def adjust_proxy_expression(self, expr, rscale, nderivatives): assert rscale == 1 - return expr * factor + return expr def __getinitargs__(self): return (self.dim, self.icomp, self.jcomp, self.kcomp, @@ -809,12 +795,12 @@ class KernelWrapper(Kernel): def is_complex_valued(self): return self.inner_kernel.is_complex_valued - def get_proxy_expression(self, scaled_dist_vec): - return self.inner_kernel.get_proxy_expression(scaled_dist_vec) + def get_expression(self, scaled_dist_vec): + return self.inner_kernel.get_expression(scaled_dist_vec) - def adjust_proxy_expression(self, expr, rscale, nderivatives, factor): + def adjust_proxy_expression(self, expr, rscale, nderivatives): return self.inner_kernel.adjust_proxy_expression( - expr, rscale, nderivatives, factor) + expr, rscale, nderivatives) def postprocess_at_source(self, expr, avec): return self.inner_kernel.postprocess_at_source(expr, avec) diff --git a/sumpy/p2p.py b/sumpy/p2p.py index cb94934d..fb0a1cf5 100644 --- a/sumpy/p2p.py +++ b/sumpy/p2p.py @@ -79,9 +79,7 @@ class P2PBase(KernelComputation, KernelCacheWrapper): sac.assign_unique("knl%d" % i, knl.postprocess_at_target( knl.postprocess_at_source( - knl.adjust_proxy_expression( - knl.get_proxy_expression(dvec), - rscale=1, nderivatives=0, factor=1), + knl.get_expression(dvec), dvec), dvec) ) diff --git a/sumpy/tools.py b/sumpy/tools.py index 64a509f3..380a8d98 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -126,7 +126,7 @@ class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker): recurrence = ( self.wrangler.try_get_recurrence_for_derivative( - next_mi, self.cache_by_mi)) + next_mi, self.cache_by_mi, rscale=1)) if recurrence is not None: expr = Add(*tuple( @@ -137,8 +137,6 @@ class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker): self.cache_by_mi[next_mi] = expr - expr = expr.subs(self.wrangler._rscale_symbol, 1) - return expr # }}} -- GitLab From 3d77ef2651404ae9cba0676c57d11a31ef07a4e8 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 17:37:12 -0500 Subject: [PATCH 34/60] Push use_rscale argument through more expansion cosntructors --- sumpy/expansion/__init__.py | 6 +++--- sumpy/expansion/local.py | 20 +++++++++++--------- sumpy/expansion/multipole.py | 20 +++++++++++--------- 3 files changed, 25 insertions(+), 21 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index f8549244..a74be289 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -449,7 +449,7 @@ class VolumeTaylorExpansion(VolumeTaylorExpansionBase): derivative_wrangler_class = FullDerivativeWrangler derivative_wrangler_cache = {} - def __init__(self, kernel, order): + def __init__(self, kernel, order, use_rscale=None): self.derivative_wrangler_key = (order, kernel.dim) @@ -458,7 +458,7 @@ class LaplaceConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): derivative_wrangler_class = LaplaceDerivativeWrangler derivative_wrangler_cache = {} - def __init__(self, kernel, order): + def __init__(self, kernel, order, use_rscale): self.derivative_wrangler_key = (order, kernel.dim) @@ -467,7 +467,7 @@ class HelmholtzConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): derivative_wrangler_class = HelmholtzDerivativeWrangler derivative_wrangler_cache = {} - def __init__(self, kernel, order): + def __init__(self, kernel, order, use_rscale=None): helmholtz_k_name = kernel.get_base_kernel().helmholtz_k_name self.derivative_wrangler_key = (order, kernel.dim, helmholtz_k_name) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index d911ab4f..296e6c43 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -206,27 +206,29 @@ class VolumeTaylorLocalExpansion( VolumeTaylorLocalExpansionBase, VolumeTaylorExpansion): - def __init__(self, kernel, order): - VolumeTaylorLocalExpansionBase.__init__(self, kernel, order) - VolumeTaylorExpansion.__init__(self, kernel, order) + def __init__(self, kernel, order, use_rscale=None): + VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale) + VolumeTaylorExpansion.__init__(self, kernel, order, use_rscale) class LaplaceConformingVolumeTaylorLocalExpansion( VolumeTaylorLocalExpansionBase, LaplaceConformingVolumeTaylorExpansion): - def __init__(self, kernel, order): - VolumeTaylorLocalExpansionBase.__init__(self, kernel, order) - LaplaceConformingVolumeTaylorExpansion.__init__(self, kernel, order) + def __init__(self, kernel, order, use_rscale=None): + VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale) + LaplaceConformingVolumeTaylorExpansion.__init__( + self, kernel, order, use_rscale) class HelmholtzConformingVolumeTaylorLocalExpansion( VolumeTaylorLocalExpansionBase, HelmholtzConformingVolumeTaylorExpansion): - def __init__(self, kernel, order): - VolumeTaylorLocalExpansionBase.__init__(self, kernel, order) - HelmholtzConformingVolumeTaylorExpansion.__init__(self, kernel, order) + def __init__(self, kernel, order, use_rscale=None): + VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale) + HelmholtzConformingVolumeTaylorExpansion.__init__( + self, kernel, order, use_rscale) # }}} diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 20176f8b..40eab2dd 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -183,27 +183,29 @@ class VolumeTaylorMultipoleExpansion( VolumeTaylorMultipoleExpansionBase, VolumeTaylorExpansion): - def __init__(self, kernel, order): - VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order) - VolumeTaylorExpansion.__init__(self, kernel, order) + def __init__(self, kernel, order, use_rscale=None): + VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale) + VolumeTaylorExpansion.__init__(self, kernel, order, use_rscale) class LaplaceConformingVolumeTaylorMultipoleExpansion( VolumeTaylorMultipoleExpansionBase, LaplaceConformingVolumeTaylorExpansion): - def __init__(self, kernel, order): - VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order) - LaplaceConformingVolumeTaylorExpansion.__init__(self, kernel, order) + def __init__(self, kernel, order, use_rscale=None): + VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale) + LaplaceConformingVolumeTaylorExpansion.__init__( + self, kernel, order, use_rscale) class HelmholtzConformingVolumeTaylorMultipoleExpansion( VolumeTaylorMultipoleExpansionBase, HelmholtzConformingVolumeTaylorExpansion): - def __init__(self, kernel, order): - VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order) - HelmholtzConformingVolumeTaylorExpansion.__init__(self, kernel, order) + def __init__(self, kernel, order, use_rscale=None): + VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale) + HelmholtzConformingVolumeTaylorExpansion.__init__( + self, kernel, order, use_rscale) # }}} -- GitLab From fe3751e00c95381070ac99d1d945c2b2a16b0a1c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 17:37:56 -0500 Subject: [PATCH 35/60] Add some missing rscale plumbing to p2e and e2p --- sumpy/e2p.py | 8 +++++++- sumpy/p2e.py | 8 +++++++- 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index ea5c382f..d3345a6c 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -310,7 +310,13 @@ class E2PFromCSR(E2PBase): def __call__(self, queue, **kwargs): knl = self.get_cached_optimized_kernel() - return knl(queue, **kwargs) + + centers = kwargs.pop("centers") + # "1" may be passed for rscale, which won't have its type + # meaningfully inferred. Make the type of rscale explicit. + rscale = centers.dtype.type(kwargs.pop("rscale")) + + return knl(queue, centers=centers, rscale=rscale, **kwargs) # }}} diff --git a/sumpy/p2e.py b/sumpy/p2e.py index ce60eea9..187f23be 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -283,10 +283,16 @@ class P2EFromCSR(P2EBase): :arg centers: :arg sources: :arg strengths: + :arg rscale: """ knl = self.get_cached_optimized_kernel() - return knl(queue, **kwargs) + centers = kwargs.pop("centers") + # "1" may be passed for rscale, which won't have its type + # meaningfully inferred. Make the type of rscale explicit. + rscale = centers.dtype.type(kwargs.pop("rscale")) + + return knl(queue, centers=centers, rscale=rscale, **kwargs) # }}} -- GitLab From d38bf02cb6f5ca50a6465c42f542cc31c3372399 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 17:38:27 -0500 Subject: [PATCH 36/60] Declare 3D biharmonic kernel 'supported' again --- sumpy/kernel.py | 1 - 1 file changed, 1 deletion(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index f77e43f7..d443c791 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -453,7 +453,6 @@ class BiharmonicKernel(ExpressionKernel): expr = r**2 * var("log")(r) scaling = 1/(8*var("pi")) elif dim == 3: - raise RuntimeError("unsupported dimensionality") expr = r scaling = 1 # FIXME: Unknown else: -- GitLab From 67a7533090ed1a2ad1322d6add2cf9bb60bf9a28 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 17:40:56 -0500 Subject: [PATCH 37/60] Adapt FMM to pass rscale (works for locals but not multipoles) --- sumpy/fmm.py | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/sumpy/fmm.py b/sumpy/fmm.py index b20c9827..ee4ba146 100644 --- a/sumpy/fmm.py +++ b/sumpy/fmm.py @@ -43,6 +43,14 @@ from sumpy import ( E2EFromCSR, E2EFromChildren, E2EFromParent) +def mpole_level_to_rscale(tree, level): + return 1 # tree.root_extent * (2**-level) + + +def local_level_to_rscale(tree, level): + return tree.root_extent * (2**-level) + + # {{{ expansion wrangler code container class SumpyExpansionWranglerCodeContainer(object): @@ -56,7 +64,7 @@ class SumpyExpansionWranglerCodeContainer(object): def __init__(self, cl_context, multipole_expansion_factory, local_expansion_factory, - out_kernels, exclude_self=False): + out_kernels, exclude_self=False, use_rscale=None): """ :arg multipole_expansion_factory: a callable of a single argument (order) that returns a multipole expansion. @@ -69,16 +77,17 @@ class SumpyExpansionWranglerCodeContainer(object): self.local_expansion_factory = local_expansion_factory self.out_kernels = out_kernels self.exclude_self = exclude_self + self.use_rscale = use_rscale self.cl_context = cl_context @memoize_method def multipole_expansion(self, order): - return self.multipole_expansion_factory(order) + return self.multipole_expansion_factory(order, self.use_rscale) @memoize_method def local_expansion(self, order): - return self.local_expansion_factory(order) + return self.local_expansion_factory(order, self.use_rscale) @memoize_method def p2m(self, tgt_order): @@ -324,6 +333,7 @@ class SumpyExpansionWrangler(object): tgt_base_ibox=level_start_ibox, wait_for=src_weights.events, + rscale=mpole_level_to_rscale(self.tree, lev), **kwargs) @@ -375,6 +385,9 @@ class SumpyExpansionWrangler(object): box_child_ids=self.tree.box_child_ids, centers=self.tree.box_centers, + src_rscale=mpole_level_to_rscale(self.tree, source_level), + tgt_rscale=mpole_level_to_rscale(self.tree, target_level), + **self.kernel_extra_kwargs) assert mpoles_res is target_mpoles_view @@ -441,6 +454,9 @@ class SumpyExpansionWrangler(object): src_box_lists=src_box_lists, centers=self.tree.box_centers, + src_rscale=mpole_level_to_rscale(self.tree, lev), + tgt_rscale=local_level_to_rscale(self.tree, lev), + wait_for=mpole_exps.events, **self.kernel_extra_kwargs) @@ -484,6 +500,8 @@ class SumpyExpansionWrangler(object): centers=self.tree.box_centers, result=pot, + rscale=mpole_level_to_rscale(self.tree, isrc_level), + wait_for=wait_for, **kwargs) @@ -530,6 +548,8 @@ class SumpyExpansionWrangler(object): tgt_expansions=target_local_exps_view, tgt_base_ibox=target_level_start_ibox, + rscale=local_level_to_rscale(self.tree, lev), + wait_for=src_weights.events, **kwargs) @@ -572,6 +592,9 @@ class SumpyExpansionWrangler(object): box_parent_ids=self.tree.box_parent_ids, centers=self.tree.box_centers, + src_rscale=local_level_to_rscale(self.tree, source_lev), + tgt_rscale=local_level_to_rscale(self.tree, target_lev), + **self.kernel_extra_kwargs) assert local_exps_res is target_local_exps_view @@ -607,6 +630,8 @@ class SumpyExpansionWrangler(object): centers=self.tree.box_centers, result=pot, + rscale=local_level_to_rscale(self.tree, lev), + wait_for=local_exps.events, **kwargs) -- GitLab From b53acc463c22b5282d9461c1f564e795bc6bbf9f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 18:06:38 -0500 Subject: [PATCH 38/60] Add use_rscale to {Y,H}2DLocalExpansion --- sumpy/expansion/local.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 296e6c43..8e4b7364 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -326,12 +326,12 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): class H2DLocalExpansion(_FourierBesselLocalExpansion): - def __init__(self, kernel, order): + def __init__(self, kernel, order, use_rscale=None): from sumpy.kernel import HelmholtzKernel assert (isinstance(kernel.get_base_kernel(), HelmholtzKernel) and kernel.dim == 2) - super(H2DLocalExpansion, self).__init__(kernel, order) + super(H2DLocalExpansion, self).__init__(kernel, order, use_rscale) from sumpy.expansion.multipole import H2DMultipoleExpansion self.mpole_expn_class = H2DMultipoleExpansion @@ -341,12 +341,12 @@ class H2DLocalExpansion(_FourierBesselLocalExpansion): class Y2DLocalExpansion(_FourierBesselLocalExpansion): - def __init__(self, kernel, order): + def __init__(self, kernel, order, use_rscale=None): from sumpy.kernel import YukawaKernel assert (isinstance(kernel.get_base_kernel(), YukawaKernel) and kernel.dim == 2) - super(Y2DLocalExpansion, self).__init__(kernel, order) + super(Y2DLocalExpansion, self).__init__(kernel, order, use_rscale) from sumpy.expansion.multipole import Y2DMultipoleExpansion self.mpole_expn_class = Y2DMultipoleExpansion -- GitLab From 0f6c95dd82e249c9c8fbbd189a5c9e941300ee3e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 18:09:09 -0500 Subject: [PATCH 39/60] Fix Py2 syntax error in sumpy.toys --- sumpy/toys.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/sumpy/toys.py b/sumpy/toys.py index 251c0803..c5f544f6 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -318,8 +318,7 @@ class PointSources(PotentialSource): evt, (potential,) = self.toy_ctx.get_p2p()( self.toy_ctx.queue, targets, self.points, [self.weights], out_host=True, - **self.toy_ctx.extra_source_and_kernel_kwargs, - ) + **self.toy_ctx.extra_source_and_kernel_kwargs) return potential -- GitLab From c17c463d68347e777f8ac50d174ab2f8c6f24fbb Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 25 Aug 2017 18:31:00 -0500 Subject: [PATCH 40/60] Add back get_global_scaling_const to LineTaylorLocalExpansion to allow sumpy.qbx to function with it --- sumpy/expansion/local.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 8e4b7364..495eed00 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -52,6 +52,11 @@ __doc__ = """ class LineTaylorLocalExpansion(LocalExpansionBase): + # This is here to conform this to enough of the kernel interface + # to make it fit into sumpy.qbx.LayerPotential. + def get_global_scaling_const(self): + return self.kernel.get_global_scaling_const() + def get_storage_index(self, k): return k -- GitLab From 0c10b084116b97e0946c55547feafdbbfee69acc Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 26 Aug 2017 21:37:19 -0500 Subject: [PATCH 41/60] Fix sumpy FMM integration with rscale --- sumpy/fmm.py | 50 ++++++++++++++++++++++++++------------------------ 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/sumpy/fmm.py b/sumpy/fmm.py index ee4ba146..77166e1b 100644 --- a/sumpy/fmm.py +++ b/sumpy/fmm.py @@ -43,11 +43,7 @@ from sumpy import ( E2EFromCSR, E2EFromChildren, E2EFromParent) -def mpole_level_to_rscale(tree, level): - return 1 # tree.root_extent * (2**-level) - - -def local_level_to_rscale(tree, level): +def level_to_rscale(tree, level): return tree.root_extent * (2**-level) @@ -333,7 +329,7 @@ class SumpyExpansionWrangler(object): tgt_base_ibox=level_start_ibox, wait_for=src_weights.events, - rscale=mpole_level_to_rscale(self.tree, lev), + rscale=level_to_rscale(self.tree, lev), **kwargs) @@ -353,18 +349,23 @@ class SumpyExpansionWrangler(object): tree = self.tree evt = None - # 2 is the last relevant source_level. - # 1 is the last relevant target_level. - # (Nobody needs a multipole on level 0, i.e. for the root box.) - for source_level in range(tree.nlevels-1, 1, -1): - start, stop = level_start_source_parent_box_nrs[ - source_level:source_level+2] - if start == stop: - continue + # nlevels-1 is the last valid level index + # nlevels-2 is the last valid level that could have children + # + # 3 is the last relevant source_level. + # 2 is the last relevant target_level. + # (because no level 1 box will be well-separated from another) + for source_level in range(tree.nlevels-1, 2, -1): target_level = source_level - 1 assert target_level > 0 + start, stop = level_start_source_parent_box_nrs[ + target_level:target_level+2] + if start == stop: + print("source", source_level, "empty") + continue + m2m = self.code.m2m( self.level_orders[source_level], self.level_orders[target_level]) @@ -385,16 +386,17 @@ class SumpyExpansionWrangler(object): box_child_ids=self.tree.box_child_ids, centers=self.tree.box_centers, - src_rscale=mpole_level_to_rscale(self.tree, source_level), - tgt_rscale=mpole_level_to_rscale(self.tree, target_level), + src_rscale=level_to_rscale(self.tree, source_level), + tgt_rscale=level_to_rscale(self.tree, target_level), **self.kernel_extra_kwargs) - assert mpoles_res is target_mpoles_view if evt is not None: mpoles.add_event(evt) + # }}} + return mpoles def eval_direct(self, target_boxes, source_box_starts, @@ -454,8 +456,8 @@ class SumpyExpansionWrangler(object): src_box_lists=src_box_lists, centers=self.tree.box_centers, - src_rscale=mpole_level_to_rscale(self.tree, lev), - tgt_rscale=local_level_to_rscale(self.tree, lev), + src_rscale=level_to_rscale(self.tree, lev), + tgt_rscale=level_to_rscale(self.tree, lev), wait_for=mpole_exps.events, @@ -500,7 +502,7 @@ class SumpyExpansionWrangler(object): centers=self.tree.box_centers, result=pot, - rscale=mpole_level_to_rscale(self.tree, isrc_level), + rscale=level_to_rscale(self.tree, isrc_level), wait_for=wait_for, @@ -548,7 +550,7 @@ class SumpyExpansionWrangler(object): tgt_expansions=target_local_exps_view, tgt_base_ibox=target_level_start_ibox, - rscale=local_level_to_rscale(self.tree, lev), + rscale=level_to_rscale(self.tree, lev), wait_for=src_weights.events, @@ -592,8 +594,8 @@ class SumpyExpansionWrangler(object): box_parent_ids=self.tree.box_parent_ids, centers=self.tree.box_centers, - src_rscale=local_level_to_rscale(self.tree, source_lev), - tgt_rscale=local_level_to_rscale(self.tree, target_lev), + src_rscale=level_to_rscale(self.tree, source_lev), + tgt_rscale=level_to_rscale(self.tree, target_lev), **self.kernel_extra_kwargs) @@ -630,7 +632,7 @@ class SumpyExpansionWrangler(object): centers=self.tree.box_centers, result=pot, - rscale=local_level_to_rscale(self.tree, lev), + rscale=level_to_rscale(self.tree, lev), wait_for=local_exps.events, -- GitLab From e1c32a7a0efcc3c8cccc2906a2658a0fedc7c303 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 26 Aug 2017 23:56:51 -0500 Subject: [PATCH 42/60] Fix equality comparison and persistent hashing of expansions for use_rscale --- sumpy/expansion/__init__.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index a74be289..d8f04e7d 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -119,12 +119,14 @@ class ExpansionBase(object): key_hash.update(type(self).__name__.encode("utf8")) key_builder.rec(key_hash, self.kernel) key_builder.rec(key_hash, self.order) + key_builder.rec(key_hash, self.use_rscale) def __eq__(self, other): return ( type(self) == type(other) and self.kernel == other.kernel - and self.order == other.order) + and self.order == other.order + and self.use_rscale == other.use_rscale) def __ne__(self, other): return not self.__eq__(other) -- GitLab From 9e9027d264c9151e646c5a1a43b71f3706194227 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 27 Aug 2017 00:01:38 -0500 Subject: [PATCH 43/60] Fix/move to the right place the expansion's propagation of the kernel interface --- sumpy/expansion/__init__.py | 7 +++++-- sumpy/expansion/local.py | 5 ----- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index d8f04e7d..325b4c34 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -66,6 +66,9 @@ class ExpansionBase(object): # {{{ propagate kernel interface + # This is here to conform this to enough of the kernel interface + # to make it fit into sumpy.qbx.LayerPotential. + @property def dim(self): return self.kernel.dim @@ -80,8 +83,8 @@ class ExpansionBase(object): def get_code_transformer(self): return self.kernel.get_code_transformer() - def get_scaling(self): - return self.kernel.get_scaling() + def get_global_scaling_const(self): + return self.kernel.get_global_scaling_const() def get_args(self): return self.kernel.get_args() diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 495eed00..8e4b7364 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -52,11 +52,6 @@ __doc__ = """ class LineTaylorLocalExpansion(LocalExpansionBase): - # This is here to conform this to enough of the kernel interface - # to make it fit into sumpy.qbx.LayerPotential. - def get_global_scaling_const(self): - return self.kernel.get_global_scaling_const() - def get_storage_index(self, k): return k -- GitLab From 7a576e553a1e3ec28da7796d944b950d63243cdf Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 27 Aug 2017 20:17:44 -0500 Subject: [PATCH 44/60] Simplify, generalize rscaling to all kernels, improve docs --- sumpy/expansion/local.py | 14 ++- sumpy/expansion/multipole.py | 23 ++++- sumpy/kernel.py | 192 +++++++++++++---------------------- sumpy/version.py | 2 +- 4 files changed, 95 insertions(+), 136 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 8e4b7364..e97addf8 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -24,7 +24,6 @@ THE SOFTWARE. from six.moves import range, zip import sumpy.symbolic as sym -from sumpy.tools import sympy_vec_subs from sumpy.expansion import ( ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, @@ -173,13 +172,12 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): src_coeff_exprs, src_expansion.get_coefficient_identifiers()): - kernel_deriv = self.kernel.adjust_proxy_expression( - sympy_vec_subs( - dvec, dvec/src_rscale, + kernel_deriv = ( + src_expansion.get_scaled_multipole( taker.diff(add_mi(deriv, term)), - ), - src_rscale, sum(deriv) + sum(term) - ) / src_rscale**sum(deriv) + dvec, src_rscale, + nderivatives=sum(deriv) + sum(term), + nderivatives_for_scaling=sum(term))) local_result.append( coeff * kernel_deriv * tgt_rscale**sum(deriv)) @@ -195,7 +193,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # canceling "rscales" closer to each other in the hope of helping # with that. result = [ - (taker.diff(mi) * tgt_rscale**sum(mi)).expand() + (taker.diff(mi) * tgt_rscale**sum(mi)).expand(deep=False) for mi in self.get_coefficient_identifiers()] logger.info("building translation operator: done") diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 40eab2dd..3b957dda 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -97,18 +97,31 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): self.derivative_wrangler.get_stored_mpole_coefficients_from_full( result, rscale)) + def get_scaled_multipole(self, expr, bvec, rscale, nderivatives, + nderivatives_for_scaling=None): + if nderivatives_for_scaling is None: + nderivatives_for_scaling = nderivatives + + if self.kernel.has_efficient_scale_adjustment: + return ( + self.kernel.adjust_for_kernel_scaling( + sympy_vec_subs( + bvec, bvec/rscale, + expr), + rscale, nderivatives) + / rscale ** (nderivatives - nderivatives_for_scaling)) + else: + return (rscale**nderivatives_for_scaling * expr) + def evaluate(self, coeffs, bvec, rscale): if not self.use_rscale: rscale = 1 taker = self.get_kernel_derivative_taker(bvec) + result = sym.Add(*tuple( coeff - * self.kernel.adjust_proxy_expression( - sympy_vec_subs( - bvec, bvec/rscale, - taker.diff(mi)), - rscale, sum(mi)) + * self.get_scaled_multipole(taker.diff(mi), bvec, rscale, sum(mi)) for coeff, mi in zip(coeffs, self.get_coefficient_identifiers()))) return result diff --git a/sumpy/kernel.py b/sumpy/kernel.py index d443c791..710dd5d5 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -106,7 +106,8 @@ class Kernel(object): .. automethod:: prepare_loopy_kernel .. automethod:: get_code_transformer .. automethod:: get_expression - .. automethod:: adjust_proxy_expression + .. attribute:: has_efficient_scale_adjustment + .. automethod:: adjust_for_kernel_scaling .. automethod:: postprocess_at_source .. automethod:: postprocess_at_target .. automethod:: get_global_scaling_const @@ -172,96 +173,75 @@ class Kernel(object): """ return lambda expr: expr - # How does rscaling work for symbolic expressions? - # ------------------------------------------------ - # - # The thing is that all this is just a big no-op. We are trying to - # differentiate - # - # $G(\frac{d}\alpha \alpha)$ - # - # while shuffling around powers of $\alpha$ to keep all terms as close as - # possible to unit-scale. - # - # We call $\delta = d/\alpha$ the unit-scaled geometry variable. - # - # We could just differentiate - # - # $ \delta \mapsto G( \delta \alpha),$ - # - # which would in principle tell us everything we need to know about the - # scaling behavior of the kernel. In this iteration I chose not to do that - # because I felt the extra $\alpha$ would put an unnecessary "clutter" - # burden on the symbolic machinery, making it scale even less well than it - # already does. I don't know whether that's true and/or conceptually - # mistaken. - # - # Instead, I decided that we'll differentiate "just $G$" and have a - # fix-up routine that provides the result of substituting $\alpha \delta$ - # into $\partial^\mu G$ after the fact by modifying the expression "from - # the outside". We then supply $\delta$ as input. The thought behind this - # is that any arithmetic in $\partial^\mu G$ is carried out on (close-to) - # unit-scale quantities. - # - # But since we've differentiated a $G$ "with the wrong scaling", there is - # also a chain-rule term $\alpha^{|\mu|}$ that needs to be multiplied in. - # This is deferred until the evaluation of the expansion. - # - # See also the MR that introduced this: - # https://gitlab.tiker.net/inducer/sumpy/merge_requests/46 + def get_expression(self, dist_vec): + r"""Return a :mod:`sympy` expression for the kernel.""" + raise NotImplementedError - def get_expression(self, scaled_dist_vec): - r"""Return a proxy expression :math:`\tilde G(\vec \delta)` where - :math:`\vec \delta` is *scaled_dist_vec*. Here, - :math:`\vec \delta =\vec d/\alpha`, - where :math:`\alpha` is a scaling factor (often written as *rscale*) - used to ensure that :math:`\delta` has roughly unit magnitude and - :math:`\vec d=\text{target} - \text{source}`. (In - an FMM, :math:`\alpha` would be chosen proportional to the - box size.) :math:`\tilde G` should be chosen so that (a) - no powers of :math:`\alpha` appear when symbolic derivatives - of it are taken with respect to :math:`\vec \delta` and - (b) the correct partial derivatives (with input scaling by :math:`\alpha` - can nonetheless be recovered) by :meth:`adjust_proxy_expression`. - - Let :math:`\hat G` be the result of calling - :meth:`adjust_proxy_expression` with - :math:`\partial_{\vec \delta}^\mu \tilde G` and with - *nderivatives* set equal to :math:`|\mu|=\sum_i \mu_i`, - where :math:`\mu` is a multi-index indicating - how many derivatives are being taken along each input axis. - - Then the following equality is required to hold: + has_efficient_scale_adjustment = False + + def adjust_for_kernel_scaling(self, expr, rscale, nderivatives): + r""" + Consider a Taylor multipole expansion: .. math:: - \hat G(\vec \delta) \alpha^{|\mu|} - = \partial_{\vec \delta}^\mu \tilde G(\vec \delta) \alpha^{|\mu|} - = \partial_{\vec d}^\mu G(\vec d). + f (x - y) = \sum_{i = 0}^{\infty} (\partial_y^i f) (x - y) \big|_{y = c} + \frac{(y - c)^i}{i!} . + + Now suppose we would like to use a scaled version :math:`g` of the + kernel :math:`f`: - In other words, :meth:`adjust_proxy_expression` is *not* - responsible for applying the scaling :math:`\alpha^{|\mu|}` - resulting from the chain rule (this is managed separately within - :mod:`sumpy` in an effort to control expansion coefficient scaling), - but it should compensate for any other effects of the - scaling of the input variable. + .. math:: - As a worked example, consider a kernel of :math:`\log(x)`. The - proxy expression would be :math:`\log(\delta)`. + \begin{eqnarray*} + f (x) & = & g (x / \alpha),\\ + f^{(i)} (x) & = & \frac{1}{\alpha^i} g^{(i)} (x / \alpha) . + \end{eqnarray*} - For ``nderivatives == 0``, :meth:`adjust_proxy_expression` - would return :math:`\tilde G + \log(\alpha)`. + where :math:`\alpha` is chosen to be on a length scale similar to + :math:`x` (for example by choosing :math:`\alpha` proporitional to the + size of the box for which the expansion is intended) so that :math:`x / + \alpha` is roughly of unit magnitude, to avoid arithmetic issues with + small arguments. This yields - For ``nderivatives > 0``, :meth:`adjust_proxy_expression` - would return :math:`\tilde G/\alpha^{\text{nderivatives}}`. + .. math:: - In each case, the added terms merely compensate for the missing - input scaling. - """ - raise NotImplementedError + f (x - y) = \sum_{i = 0}^{\infty} (\partial_y^i g) + \left( \frac{x - y}{\alpha} \right) \Bigg|_{y = c} + \cdot + \frac{(y - c)^i}{\alpha^i \cdot i!}. + + Observe that the :math:`(y - c)` term is now scaled to unit magnitude, + as is the argument of :math:`g`. - def adjust_proxy_expression(self, expr, rscale, nderivatives): - r"""See :meth:`get_expression`. + With :math:`\xi = x / \alpha`, we find + + .. math:: + + \begin{eqnarray*} + g (\xi) & = & f (\alpha \xi),\\ + g^{(i)} (\xi) & = & \alpha^i f^{(i)} (\alpha \xi) . + \end{eqnarray*} + + Generically for all kernels, :math:`f^{(i)} (\alpha \xi)` is computable + by taking a sufficient number of symbolic derivatives of :math:`f` and + providing :math:`\alpha \xi = x` as the argument. + + Now, for some kernels, like :math:`f (x) = \log x`, the powers of + :math:`\alpha^i` from the chain rule cancel with the ones from the + argument substituted into the kernel derivative: + + .. math:: + + g^{(i)} (\xi) = \alpha^i f^{(i)} (\alpha \xi) = C \cdot \alpha^i \cdot + \frac{1}{(\alpha x)^i} \quad (i > 0), + + making them what you might call *scale-invariant*. In this case, one + may set :attr:`has_efficient_scale_adjustment`. For these kernels only, + :meth:`adjust_for_kernel_scaling` provides a shortcut for scaled kernel + derivative computation. Given :math:`f^{(i)}` as *expr*, it directly + returns an expression for :math:`g^{(i)}`, where :math:`i` is given + as *nderivatives*. """ raise NotImplementedError @@ -421,7 +401,9 @@ class LaplaceKernel(ExpressionKernel): global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives): + has_efficient_scale_adjustment = True + + def adjust_for_kernel_scaling(self, expr, rscale, nderivatives): if self.dim == 2: if nderivatives == 0: import sympy as sp @@ -464,10 +446,6 @@ class BiharmonicKernel(ExpressionKernel): global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives): - assert rscale == 1 - return expr - def __getinitargs__(self): return (self.dim,) @@ -511,19 +489,6 @@ class HelmholtzKernel(ExpressionKernel): self.helmholtz_k_name = helmholtz_k_name self.allow_evanescent = allow_evanescent - def adjust_proxy_expression(self, expr, rscale, nderivatives): - import sympy as sp - result = expr.subs( - sp.Symbol(self.helmholtz_k_name), - sp.Symbol(self.helmholtz_k_name) * rscale) - - if self.dim == 2: - return result - elif self.dim == 3: - return result / rscale - else: - raise RuntimeError("unsupported dimensionality") - def __getinitargs__(self): return (self.dim, self.helmholtz_k_name, self.allow_evanescent) @@ -589,19 +554,6 @@ class YukawaKernel(ExpressionKernel): self.yukawa_lambda_name = yukawa_lambda_name - def adjust_proxy_expression(self, expr, rscale, nderivatives): - import sympy as sp - result = expr.subs( - sp.Symbol(self.helmholtz_k_name), - sp.Symbol(self.helmholtz_k_name) * rscale) - - if self.dim == 2: - return result - elif self.dim == 3: - return result / rscale - else: - raise RuntimeError("unsupported dimensionality") - def __getinitargs__(self): return (self.dim, self.yukawa_lambda_name) @@ -678,10 +630,6 @@ class StokesletKernel(ExpressionKernel): global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives): - assert rscale == 1 - return expr - def __getinitargs__(self): return (self.dim, self.icomp, self.jcomp, self.viscosity_mu_name) @@ -747,10 +695,6 @@ class StressletKernel(ExpressionKernel): global_scaling_const=scaling, is_complex_valued=False) - def adjust_proxy_expression(self, expr, rscale, nderivatives): - assert rscale == 1 - return expr - def __getinitargs__(self): return (self.dim, self.icomp, self.jcomp, self.kcomp, self.viscosity_mu_name) @@ -797,8 +741,12 @@ class KernelWrapper(Kernel): def get_expression(self, scaled_dist_vec): return self.inner_kernel.get_expression(scaled_dist_vec) - def adjust_proxy_expression(self, expr, rscale, nderivatives): - return self.inner_kernel.adjust_proxy_expression( + @property + def has_efficient_scale_adjustment(self): + return self.inner_kernel.has_efficient_scale_adjustment + + def adjust_for_kernel_scaling(self, expr, rscale, nderivatives): + return self.inner_kernel.adjust_for_kernel_scaling( expr, rscale, nderivatives) def postprocess_at_source(self, expr, avec): diff --git a/sumpy/version.py b/sumpy/version.py index 72c4bd46..e3df92b8 100644 --- a/sumpy/version.py +++ b/sumpy/version.py @@ -25,4 +25,4 @@ VERSION = (2016, 1) VERSION_STATUS = "beta1" VERSION_TEXT = ".".join(str(x) for x in VERSION) + VERSION_STATUS -KERNEL_VERSION = 16 +KERNEL_VERSION = 17 -- GitLab From b56d82b9533a68c2575d34d5582857c65b02dddf Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 27 Aug 2017 20:19:22 -0500 Subject: [PATCH 45/60] Get rid of coefficient scaling experiment --- examples/coeff-scaling-experiment.py | 65 ---------------------------- 1 file changed, 65 deletions(-) delete mode 100644 examples/coeff-scaling-experiment.py diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py deleted file mode 100644 index a9bbf32f..00000000 --- a/examples/coeff-scaling-experiment.py +++ /dev/null @@ -1,65 +0,0 @@ -import pyopencl as cl -import sumpy.toys as t -import numpy as np -from sumpy.visualization import FieldPlotter -import matplotlib.pyplot as plt - - -# FIXME: Get rid of this once everything is working - -def main(): - dim = 2 - order = 10 - - from sumpy.kernel import ( # noqa: F401 - YukawaKernel, HelmholtzKernel, LaplaceKernel, - BiharmonicKernel, StokesletKernel, StressletKernel, - AxisTargetDerivative) - tctx = t.ToyContext( - cl.create_some_context(), - LaplaceKernel(dim), - #AxisTargetDerivative(0, LaplaceKernel(dim)), - #YukawaKernel(dim), extra_source_kwargs={"lam": 5}, - #HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3}, - #BiharmonicKernel(dim), - #StokesletKernel(dim, 1, 1), extra_source_kwargs={"mu": 0.3}, - #StressletKernel(dim, 1, 1, 0), extra_source_kwargs={"mu": 0.3}, - ) - - np.random.seed(12) - scale = 2**(-14) - #scale = 1 - pts = np.random.rand(dim, 50) - 0.5 - pt_src = t.PointSources( - tctx, - scale * pts, - np.ones(50)) - - mctr = scale*np.array([0., 0, 0])[:dim] - mctr1 = scale*np.array([0.2, 0, 0])[:dim] - mexp1 = t.multipole_expand(pt_src, mctr1, order=order, rscale=scale) - mexp = t.multipole_expand(mexp1, mctr, order=order, rscale=2*scale) - - lctr1 = scale*np.array([2.8, 0, 0])[:dim] - lctr = scale*np.array([2.5, 0, 0])[:dim] - lexp1 = t.local_expand(mexp, lctr1, order=order, rscale=scale) - lexp = t.local_expand(lexp1, lctr, order=order, rscale=2*scale) - - #print(mexp.coeffs) - print(lexp.coeffs) - - diff = lexp - pt_src - - diag = np.sqrt(dim) - print(t.l_inf(diff, scale*0.5*diag, center=lctr) - / t.l_inf(pt_src, scale*0.5*diag, center=lctr)) - if 0: - fp = FieldPlotter(lexp.center, extent=scale*0.5) - - t.logplot(fp, diff, cmap="jet", vmin=-16, vmax=-7) - plt.colorbar() - plt.show() - - -if __name__ == "__main__": - main() -- GitLab From b25c9daa225633f8d4dd1cd4b0c765f6cc422ec3 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 28 Aug 2017 10:54:40 -0500 Subject: [PATCH 46/60] Add use_rscale to ExpansionBase.with_kernel --- sumpy/expansion/__init__.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 2a370c33..79826664 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -95,8 +95,7 @@ class ExpansionBase(object): # }}} def with_kernel(self, kernel): - return type(self)(kernel, self.order) - # FIXME: add self.use_rscale once the rscale MR is in + return type(self)(kernel, self.order, self.use_rscale) def __len__(self): return len(self.get_coefficient_identifiers()) -- GitLab From 5a15a2c91f91ab0f0ca69af1a004566e1ce98dbb Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 28 Aug 2017 17:16:06 -0500 Subject: [PATCH 47/60] (Slightly) better document base expansion interface --- sumpy/expansion/__init__.py | 23 +++++++++++++++++++++-- sumpy/expansion/local.py | 12 ++++++------ sumpy/expansion/multipole.py | 12 ++++++------ 3 files changed, 33 insertions(+), 14 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 79826664..fd734697 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -50,6 +50,15 @@ logger = logging.getLogger(__name__) # {{{ base class class ExpansionBase(object): + """ + .. automethod:: with_kernel + .. automethod:: __len__ + .. automethod:: get_coefficient_identifiers + .. automethod:: coefficients_from_source + .. automethod:: translate_from + .. automethod:: __eq__ + .. automethod:: __ne__ + """ def __init__(self, kernel, order, use_rscale=None): # Don't be tempted to remove target derivatives here. @@ -100,7 +109,13 @@ class ExpansionBase(object): def __len__(self): return len(self.get_coefficient_identifiers()) - def coefficients_from_source(self, avec, bvec): + def get_coefficient_identifiers(self): + """ + Returns the identifiers of the coefficients that actually get stored. + """ + raise NotImplementedError + + def coefficients_from_source(self, avec, bvec, rscale): """Form an expansion from a source point. :arg avec: vector from source to center. @@ -112,7 +127,7 @@ class ExpansionBase(object): """ raise NotImplementedError - def evaluate(self, coeffs, bvec): + def evaluate(self, coeffs, bvec, rscale): """ :return: a :mod:`sympy` expression corresponding to the evaluated expansion with the coefficients @@ -121,6 +136,10 @@ class ExpansionBase(object): raise NotImplementedError + def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, + dvec, tgt_rscale): + raise NotImplementedError + def update_persistent_hash(self, key_hash, key_builder): key_hash.update(type(self).__name__.encode("utf8")) key_builder.rec(key_hash, self.kernel) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index e97addf8..317900fb 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -201,8 +201,8 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): class VolumeTaylorLocalExpansion( - VolumeTaylorLocalExpansionBase, - VolumeTaylorExpansion): + VolumeTaylorExpansion, + VolumeTaylorLocalExpansionBase): def __init__(self, kernel, order, use_rscale=None): VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale) @@ -210,8 +210,8 @@ class VolumeTaylorLocalExpansion( class LaplaceConformingVolumeTaylorLocalExpansion( - VolumeTaylorLocalExpansionBase, - LaplaceConformingVolumeTaylorExpansion): + LaplaceConformingVolumeTaylorExpansion, + VolumeTaylorLocalExpansionBase): def __init__(self, kernel, order, use_rscale=None): VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale) @@ -220,8 +220,8 @@ class LaplaceConformingVolumeTaylorLocalExpansion( class HelmholtzConformingVolumeTaylorLocalExpansion( - VolumeTaylorLocalExpansionBase, - HelmholtzConformingVolumeTaylorExpansion): + HelmholtzConformingVolumeTaylorExpansion, + VolumeTaylorLocalExpansionBase): def __init__(self, kernel, order, use_rscale=None): VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 3b957dda..da9baf1d 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -193,8 +193,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): class VolumeTaylorMultipoleExpansion( - VolumeTaylorMultipoleExpansionBase, - VolumeTaylorExpansion): + VolumeTaylorExpansion, + VolumeTaylorMultipoleExpansionBase): def __init__(self, kernel, order, use_rscale=None): VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale) @@ -202,8 +202,8 @@ class VolumeTaylorMultipoleExpansion( class LaplaceConformingVolumeTaylorMultipoleExpansion( - VolumeTaylorMultipoleExpansionBase, - LaplaceConformingVolumeTaylorExpansion): + LaplaceConformingVolumeTaylorExpansion, + VolumeTaylorMultipoleExpansionBase): def __init__(self, kernel, order, use_rscale=None): VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale) @@ -212,8 +212,8 @@ class LaplaceConformingVolumeTaylorMultipoleExpansion( class HelmholtzConformingVolumeTaylorMultipoleExpansion( - VolumeTaylorMultipoleExpansionBase, - HelmholtzConformingVolumeTaylorExpansion): + HelmholtzConformingVolumeTaylorExpansion, + VolumeTaylorMultipoleExpansionBase): def __init__(self, kernel, order, use_rscale=None): VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale) -- GitLab From 52a97236f7c854529be8d6341a6a68a6083a4cf2 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 28 Aug 2017 17:16:24 -0500 Subject: [PATCH 48/60] Fix direct-QBX for rscale --- sumpy/qbx.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 8892739e..a452273e 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -61,7 +61,9 @@ def stringify_expn_index(i): def expand(expansion_nr, sac, expansion, avec, bvec): - coefficients = expansion.coefficients_from_source(avec, bvec) + rscale = sym.Symbol("rscale") + + coefficients = expansion.coefficients_from_source(avec, bvec, rscale) assigned_coeffs = [ sym.Symbol( @@ -71,7 +73,7 @@ def expand(expansion_nr, sac, expansion, avec, bvec): for i in expansion.get_coefficient_identifiers()] return sac.assign_unique("expn%d_result" % expansion_nr, - expansion.evaluate(assigned_coeffs, bvec)) + expansion.evaluate(assigned_coeffs, bvec, rscale)) # {{{ layer potential computation @@ -101,6 +103,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): return """ <> a[idim] = center[idim,itgt] - src[idim,isrc] {dup=idim} <> b[idim] = tgt[idim,itgt] - center[idim,itgt] {dup=idim} + <> rscale = expansion_radii[itgt] """ def get_src_tgt_arguments(self): @@ -111,6 +114,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): shape=(self.dim, "ntargets"), order="C"), lp.GlobalArg("center", None, shape=(self.dim, "ntargets"), order="C"), + lp.GlobalArg("expansion_radii", None, shape="ntargets"), lp.ValueArg("nsources", None), lp.ValueArg("ntargets", None), ] @@ -242,7 +246,8 @@ class LayerPotential(LayerPotentialBase): for iknl in range(len(self.expansions)) ] - def __call__(self, queue, targets, sources, centers, strengths, **kwargs): + def __call__(self, queue, targets, sources, centers, strengths, expansion_radii, + **kwargs): """ :arg strengths: are required to have area elements and quadrature weights already multiplied in. @@ -253,7 +258,8 @@ class LayerPotential(LayerPotentialBase): for i, dens in enumerate(strengths): kwargs["strength_%d" % i] = dens - return knl(queue, src=sources, tgt=targets, center=centers, **kwargs) + return knl(queue, src=sources, tgt=targets, center=centers, + expansion_radii=expansion_radii, **kwargs) # }}} @@ -283,7 +289,7 @@ class LayerPotentialMatrixGenerator(LayerPotentialBase): for iknl in range(len(self.expansions)) ] - def __call__(self, queue, targets, sources, centers, **kwargs): + def __call__(self, queue, targets, sources, centers, expansion_radii, **kwargs): knl = self.get_optimized_kernel() return knl(queue, src=sources, tgt=targets, center=centers, -- GitLab From 53eab946e8fc50b8d2a57a3d30ae9236753def44 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 28 Aug 2017 19:28:29 -0500 Subject: [PATCH 49/60] Fix line-Taylor QBX tests for rscale --- sumpy/expansion/local.py | 6 ++++-- test/test_qbx.py | 5 ++++- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 317900fb..9e891a2d 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -57,7 +57,8 @@ class LineTaylorLocalExpansion(LocalExpansionBase): def get_coefficient_identifiers(self): return list(range(self.order+1)) - def coefficients_from_source(self, avec, bvec): + def coefficients_from_source(self, avec, bvec, rscale): + # no point in heeding rscale here--just ignore it if bvec is None: raise RuntimeError("cannot use line-Taylor expansions in a setting " "where the center-target vector is not known at coefficient " @@ -97,7 +98,8 @@ class LineTaylorLocalExpansion(LocalExpansionBase): .subs("tau", 0) for i in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec): + def evaluate(self, coeffs, bvec, rscale): + # no point in heeding rscale here--just ignore it from pytools import factorial return sym.Add(*( coeffs[self.get_storage_index(i)] / factorial(i) diff --git a/test/test_qbx.py b/test/test_qbx.py index aa86edf8..3758a5d7 100644 --- a/test/test_qbx.py +++ b/test/test_qbx.py @@ -82,8 +82,11 @@ def test_direct(ctx_getter): radius = 7 * h centers = unit_circle * (1 - radius) + expansion_radii = np.ones(n) * radius + strengths = (sigma * h,) - evt, (result_qbx,) = lpot(queue, targets, sources, centers, strengths) + evt, (result_qbx,) = lpot(queue, targets, sources, centers, strengths, + expansion_radii=expansion_radii) eocrec.add_data_point(h, np.max(np.abs(result_ref - result_qbx))) -- GitLab From 08bb556390130f80b15d93978d18bdcc4723d306 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 28 Aug 2017 19:39:34 -0500 Subject: [PATCH 50/60] Bump kernel version for direct-QBX kernel changes --- sumpy/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/version.py b/sumpy/version.py index e3df92b8..97d6e94e 100644 --- a/sumpy/version.py +++ b/sumpy/version.py @@ -25,4 +25,4 @@ VERSION = (2016, 1) VERSION_STATUS = "beta1" VERSION_TEXT = ".".join(str(x) for x in VERSION) + VERSION_STATUS -KERNEL_VERSION = 17 +KERNEL_VERSION = 18 -- GitLab From 490e32eb9ac3d212dcbabcead1245b3c7fee8f6a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 28 Aug 2017 20:23:37 -0500 Subject: [PATCH 51/60] Fix line Taylor coefficient complexity test after rscale --- test/test_codegen.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_codegen.py b/test/test_codegen.py index fb165598..0a02b6b4 100644 --- a/test/test_codegen.py +++ b/test/test_codegen.py @@ -91,7 +91,7 @@ def test_line_taylor_coeff_growth(): expn = LineTaylorLocalExpansion(LaplaceKernel(2), order) avec = make_sym_vector("a", 2) bvec = make_sym_vector("b", 2) - coeffs = expn.coefficients_from_source(avec, bvec) + coeffs = expn.coefficients_from_source(avec, bvec, rscale=1) sym2pymbolic = SympyToPymbolicMapper() coeffs_pymbolic = [sym2pymbolic(c) for c in coeffs] -- GitLab From 0d67f0bdc55faad802357dbefcc8c95262a4c239 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 29 Aug 2017 10:58:06 -0500 Subject: [PATCH 52/60] Change sympy references to sumpy.symbolic --- sumpy/e2p.py | 2 +- sumpy/expansion/__init__.py | 2 +- sumpy/kernel.py | 2 +- sumpy/p2e.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index 79f8c982..5a5c1154 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -81,7 +81,7 @@ class E2PBase(KernelCacheWrapper): from sumpy.symbolic import make_sym_vector bvec = make_sym_vector("b", self.dim) - import sympy as sp + import sumpy.symbolic as sp rscale = sp.Symbol("rscale") from sumpy.assignment_collection import SymbolicAssignmentCollection diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index fd734697..bbcd869c 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -29,7 +29,7 @@ import logging from pytools import memoize_method import sumpy.symbolic as sym from sumpy.tools import MiDerivativeTaker -import sympy as sp +import sumpy.symbolic as sp from collections import defaultdict diff --git a/sumpy/kernel.py b/sumpy/kernel.py index a17d0c00..c655c6f3 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -406,7 +406,7 @@ class LaplaceKernel(ExpressionKernel): def adjust_for_kernel_scaling(self, expr, rscale, nderivatives): if self.dim == 2: if nderivatives == 0: - import sympy as sp + import sumpy.symbolic as sp return (expr + sp.log(rscale)) else: return expr diff --git a/sumpy/p2e.py b/sumpy/p2e.py index 66327fa5..d9354da7 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -77,7 +77,7 @@ class P2EBase(KernelCacheWrapper): from sumpy.symbolic import make_sym_vector avec = make_sym_vector("a", self.dim) - import sympy as sp + import sumpy.symbolic as sp rscale = sp.Symbol("rscale") from sumpy.assignment_collection import SymbolicAssignmentCollection -- GitLab From f5148f91ee337bbe9532f63efbf89b8016046588 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 29 Aug 2017 12:05:43 -0500 Subject: [PATCH 53/60] Unify the two vector-subst implementations --- sumpy/expansion/multipole.py | 8 ++++---- sumpy/symbolic.py | 13 +++++++------ sumpy/tools.py | 10 ---------- 3 files changed, 11 insertions(+), 20 deletions(-) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index da9baf1d..4db9250b 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -25,7 +25,7 @@ THE SOFTWARE. from six.moves import range, zip import sumpy.symbolic as sym # noqa -from sumpy.tools import sympy_vec_subs +from sumpy.symbolic import vector_subs from sumpy.expansion import ( ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, HelmholtzConformingVolumeTaylorExpansion) @@ -105,9 +105,9 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): if self.kernel.has_efficient_scale_adjustment: return ( self.kernel.adjust_for_kernel_scaling( - sympy_vec_subs( - bvec, bvec/rscale, - expr), + vector_subs( + expr, + bvec, bvec/rscale), rscale, nderivatives) / rscale ** (nderivatives - nderivatives_for_scaling)) else: diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py index 93295bca..f0c16a0b 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -25,7 +25,6 @@ THE SOFTWARE. import six from six.moves import range -from six.moves import zip import numpy as np from pymbolic.mapper import IdentityMapper as IdentityMapperBase @@ -193,12 +192,14 @@ def make_sym_vector(name, components): [sym.Symbol("%s%d" % (name, i)) for i in range(components)]) -def vector_subs(expr, old, new): - result = expr - for old_i, new_i in zip(old, new): - result = result.subs(old_i, new_i) +def vector_subs(expr, from_vec, to_vec): + subs_pairs = [] + assert (from_vec.rows, from_vec.cols) == (to_vec.rows, to_vec.cols) + for irow in range(from_vec.rows): + for icol in range(from_vec.cols): + subs_pairs.append((from_vec[irow, icol], to_vec[irow, icol])) - return result + return expr.subs(subs_pairs) def find_power_of(base, prod): diff --git a/sumpy/tools.py b/sumpy/tools.py index 380a8d98..65cd711e 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -439,14 +439,4 @@ def my_syntactic_subs(expr, subst_dict): return expr -def sympy_vec_subs(from_vec, to_vec, expr): - subs_pairs = [] - assert (from_vec.rows, from_vec.cols) == (to_vec.rows, to_vec.cols) - for irow in range(from_vec.rows): - for icol in range(from_vec.cols): - subs_pairs.append((from_vec[irow, icol], to_vec[irow, icol])) - - return expr.subs(subs_pairs) - - # vim: fdm=marker -- GitLab From 3e0b27434671a1481e5e537c9bf0f271a5ed64b1 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 29 Aug 2017 12:06:01 -0500 Subject: [PATCH 54/60] Implement sumpy.symbolic.log --- sumpy/symbolic.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py index f0c16a0b..2effca15 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -178,6 +178,9 @@ def checked_cse(exprs, symbols=None): # }}} +log = sym.log + + def sym_real_norm_2(x): return sym.sqrt((x.T*x)[0, 0]) -- GitLab From bd66344b3c790571b0cd0c3f7e387d3c9af028fc Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 29 Aug 2017 12:11:39 -0500 Subject: [PATCH 55/60] Add comment about use_rscale defaulting --- sumpy/expansion/__init__.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index bbcd869c..9fa81ce5 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -476,7 +476,8 @@ class VolumeTaylorExpansion(VolumeTaylorExpansionBase): derivative_wrangler_class = FullDerivativeWrangler derivative_wrangler_cache = {} - def __init__(self, kernel, order, use_rscale=None): + # not user-facing, be strict about having to pass use_rscale + def __init__(self, kernel, order, use_rscale): self.derivative_wrangler_key = (order, kernel.dim) @@ -485,6 +486,7 @@ class LaplaceConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): derivative_wrangler_class = LaplaceDerivativeWrangler derivative_wrangler_cache = {} + # not user-facing, be strict about having to pass use_rscale def __init__(self, kernel, order, use_rscale): self.derivative_wrangler_key = (order, kernel.dim) @@ -494,7 +496,8 @@ class HelmholtzConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): derivative_wrangler_class = HelmholtzDerivativeWrangler derivative_wrangler_cache = {} - def __init__(self, kernel, order, use_rscale=None): + # not user-facing, be strict about having to pass use_rscale + def __init__(self, kernel, order, use_rscale): helmholtz_k_name = kernel.get_base_kernel().helmholtz_k_name self.derivative_wrangler_key = (order, kernel.dim, helmholtz_k_name) -- GitLab From e286b81909b58cdd6a5e7c482a3e3a07569af0d0 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 29 Aug 2017 12:23:42 -0500 Subject: [PATCH 56/60] Tweak usage of global kernel scaling constant in description of adjust_for_kernel_scaling --- sumpy/kernel.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index c655c6f3..cbfdfaeb 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -227,13 +227,13 @@ class Kernel(object): by taking a sufficient number of symbolic derivatives of :math:`f` and providing :math:`\alpha \xi = x` as the argument. - Now, for some kernels, like :math:`f (x) = \log x`, the powers of + Now, for some kernels, like :math:`f (x) = C \log x`, the powers of :math:`\alpha^i` from the chain rule cancel with the ones from the argument substituted into the kernel derivative: .. math:: - g^{(i)} (\xi) = \alpha^i f^{(i)} (\alpha \xi) = C \cdot \alpha^i \cdot + g^{(i)} (\xi) = \alpha^i f^{(i)} (\alpha \xi) = C' \cdot \alpha^i \cdot \frac{1}{(\alpha x)^i} \quad (i > 0), making them what you might call *scale-invariant*. In this case, one -- GitLab From 39e933da9b4587d4f5ac3ce296631d831254e8e4 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 29 Aug 2017 12:25:01 -0500 Subject: [PATCH 57/60] Document rscale in description of adjust_for_kernel_scaling --- sumpy/kernel.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index cbfdfaeb..ddfba682 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -242,6 +242,8 @@ class Kernel(object): derivative computation. Given :math:`f^{(i)}` as *expr*, it directly returns an expression for :math:`g^{(i)}`, where :math:`i` is given as *nderivatives*. + + :arg rscale: The scaling parameter :math:`\alpha` above. """ raise NotImplementedError -- GitLab From ac45874a396bde1635d020d23f18f5b8a255a325 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 29 Aug 2017 12:28:15 -0500 Subject: [PATCH 58/60] Update DerivativeWrangler abstract methods to include rscale parameters --- sumpy/expansion/__init__.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 9fa81ce5..6b986233 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -170,10 +170,12 @@ class DerivativeWrangler(object): def get_coefficient_identifiers(self): raise NotImplementedError - def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives): + def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives, + rscale): raise NotImplementedError - def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients): + def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients, + rscale): raise NotImplementedError def get_derivative_taker(self, expr, var_list): -- GitLab From 9f2e5b3e880e92b0919e867fcb6d33bf53186356 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 29 Aug 2017 12:56:22 -0500 Subject: [PATCH 59/60] Fix LayerPotentialMatrixGenerator for rscale --- sumpy/qbx.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/qbx.py b/sumpy/qbx.py index a452273e..b5806ef2 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -293,7 +293,7 @@ class LayerPotentialMatrixGenerator(LayerPotentialBase): knl = self.get_optimized_kernel() return knl(queue, src=sources, tgt=targets, center=centers, - **kwargs) + expansion_radii=expansion_radii, **kwargs) # }}} -- GitLab From d5900e49690e397ac52d7d885537d22d03ffe9bd Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 29 Aug 2017 16:56:38 -0500 Subject: [PATCH 60/60] Implement sumpy.symbolic.log using SYMBOLIC_API --- sumpy/symbolic.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py index 2effca15..55f91284 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -76,7 +76,7 @@ _find_symbolic_backend() # Symbolic API common to SymEngine and sympy. # Before adding a function here, make sure it's present in both modules. SYMBOLIC_API = """ -Add Basic Mul Pow exp sqrt symbols sympify cos sin atan2 Function Symbol +Add Basic Mul Pow exp sqrt log symbols sympify cos sin atan2 Function Symbol Derivative Integer Matrix Subs I pi functions""".split() if USE_SYMENGINE: @@ -178,9 +178,6 @@ def checked_cse(exprs, symbols=None): # }}} -log = sym.log - - def sym_real_norm_2(x): return sym.sqrt((x.T*x)[0, 0]) -- GitLab