diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index a8f78cfef8fb56fb8cda90ca97e6b2ac91813a12..5551537a78451f3d9dc5851cd1b83c7f0f428f4e 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -161,12 +161,12 @@ class ExpansionBase(object): def __ne__(self, other): return not self.__eq__(other) - def get_kernel_derivative_taker(self, dvec): + def get_kernel_derivative_taker(self, dvec, rscale=1): """Return a MiDerivativeTaker instance that supports taking derivatives of the kernel with respect to dvec """ from sumpy.tools import MiDerivativeTaker - return MiDerivativeTaker(self.kernel.get_expression(dvec), dvec) + return MiDerivativeTaker(self.kernel.get_expression(dvec), dvec, rscale) # }}} @@ -648,9 +648,9 @@ class LaplaceConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): def __init__(self, kernel, order, use_rscale): self.expansion_terms_wrangler_key = (order, kernel.dim) - def get_kernel_derivative_taker(self, dvec): + def get_kernel_derivative_taker(self, dvec, rscale=1): from sumpy.tools import LaplaceDerivativeTaker - return LaplaceDerivativeTaker(self.kernel.get_expression(dvec), dvec) + return LaplaceDerivativeTaker(self.kernel.get_expression(dvec), dvec, rscale) class HelmholtzConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index c52d00160b5216c0840c772eed2bb3aaef8b3ada..032516efc3635a7d53c47bd4d0fdb62e2f3ef53f 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 pytools import single_valued from sumpy.expansion import ( ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, @@ -123,11 +124,22 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): from sumpy.tools import MiDerivativeTakerWrapper result = [] - taker = self.get_kernel_derivative_taker(avec) + taker = self.get_kernel_derivative_taker(avec, rscale) + expr_dict = {(0,)*self.dim: 1} + expr_dict = self.kernel.get_derivative_transformation_at_source(expr_dict) + pp_nderivatives = single_valued(sum(mi) for mi in expr_dict.keys()) + for mi in self.get_coefficient_identifiers(): wrapper = MiDerivativeTakerWrapper(taker, mi) mi_expr = self.kernel.postprocess_at_source(wrapper, avec) - expr = mi_expr * rscale ** sum(mi) + # By passing `rscale` to the derivative taker we are taking a scaled + # version of the derivative which is `expr.diff(mi)*rscale**sum(mi)` + # which might be implemented efficiently for kernels like Laplace. + # One caveat is that `postprocess_at_source` might take more + # derivatives which would multiply the expression by more `rscale`s + # than necessary as the derivative taker does not know about + # `postprocess_at_source`. This is corrected by dividing by `rscale`. + expr = mi_expr / rscale ** pp_nderivatives result.append(expr) return result @@ -175,8 +187,6 @@ 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, _fft_uneval_expr from pytools import generate_nonnegative_integer_tuples_below as gnitb @@ -205,17 +215,12 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # The vector has the kernel derivatives and depends only on the distance # between the two centers - taker = src_expansion.get_kernel_derivative_taker(dvec) + taker = src_expansion.get_kernel_derivative_taker(dvec, src_rscale) vector_stored = [] # Calculate the kernel derivatives for the compressed set for term in \ srcplusderiv_terms_wrangler.get_coefficient_identifiers(): - kernel_deriv = src_expansion.get_scaled_multipole( - taker.diff(term), - dvec, src_rscale, - nderivatives=sum(term), - nderivatives_for_scaling=sum(term), - ) + kernel_deriv = taker.diff(term) vector_stored.append(kernel_deriv) # Calculate the kernel derivatives for the full set vector_full = \ diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index c1514357a415554f693eb961b699de412bb773c4..de9aa500b0c19863ae0c539c0cf6ad4ae86e4ad6 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -80,22 +80,6 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): self.expansion_terms_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( - vector_xreplace( - expr, - bvec, bvec * rscale**-1), - rscale, nderivatives) - / rscale ** (nderivatives - nderivatives_for_scaling)) - else: - return (rscale**nderivatives_for_scaling * expr) - def evaluate(self, coeffs, bvec, rscale, sac, knl=None): from sumpy.tools import MiDerivativeTakerWrapper from pytools import single_valued @@ -104,7 +88,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): if knl is None: knl = self.kernel - taker = self.get_kernel_derivative_taker(bvec) + taker = self.get_kernel_derivative_taker(bvec, rscale) expr_dict = {(0,)*self.dim: 1} expr_dict = knl.get_derivative_transformation_at_target(expr_dict) pp_nderivatives = single_valued(sum(mi) for mi in expr_dict.keys()) @@ -113,8 +97,9 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): for coeff, mi in zip(coeffs, self.get_coefficient_identifiers()): wrapper = MiDerivativeTakerWrapper(taker, mi) mi_expr = knl.postprocess_at_target(wrapper, bvec) - expr = coeff * self.get_scaled_multipole(mi_expr, bvec, rscale, - sum(mi) + pp_nderivatives, sum(mi)) + # For details about this correction, see the explanation at + # VolumeTaylorLocalExpansionBase.coefficients_from_source + expr = coeff * mi_expr / rscale**pp_nderivatives result.append(expr) result = sym.Add(*tuple(result)) @@ -365,7 +350,7 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): for l in self.get_coefficient_identifiers()) def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, - dvec, tgt_rscale): + dvec, tgt_rscale, sac): if not isinstance(src_expansion, type(self)): raise RuntimeError("do not know how to translate %s to %s" % (type(src_expansion).__name__, diff --git a/sumpy/kernel.py b/sumpy/kernel.py index e470e48d99b1c601f9ef3c881d6058b8f39f5d0c..092cbb4910d900c6d13e96356746ab7df2867d85 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -123,8 +123,6 @@ class Kernel(object): .. automethod:: prepare_loopy_kernel .. automethod:: get_code_transformer .. automethod:: get_expression - .. attribute:: has_efficient_scale_adjustment - .. automethod:: adjust_for_kernel_scaling .. automethod:: postprocess_at_source .. automethod:: postprocess_at_target .. automethod:: get_global_scaling_const @@ -193,77 +191,6 @@ class Kernel(object): r"""Return a :mod:`sympy` expression for the kernel.""" raise NotImplementedError - has_efficient_scale_adjustment = False - - def adjust_for_kernel_scaling(self, expr, rscale, nderivatives): - r""" - Consider a Taylor multipole expansion: - - .. math:: - - 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`: - - .. math:: - - \begin{eqnarray*} - f (x) & = & g (x / \alpha),\\ - f^{(i)} (x) & = & \frac{1}{\alpha^i} g^{(i)} (x / \alpha) . - \end{eqnarray*} - - 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 - - .. math:: - - 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`. - - 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) = 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 - \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*. - - :arg rscale: The scaling parameter :math:`\alpha` above. - """ - - raise NotImplementedError - def _diff(self, expr, vec, mi): """Take the derivative of an expression or a MiDerivativeTakerWrapper """ @@ -467,22 +394,6 @@ class LaplaceKernel(ExpressionKernel): global_scaling_const=scaling, is_complex_valued=False) - has_efficient_scale_adjustment = True - - def adjust_for_kernel_scaling(self, expr, rscale, nderivatives): - if self.dim == 2: - if nderivatives == 0: - import sumpy.symbolic as sp - return (expr + sp.log(rscale)) - else: - return expr - - elif self.dim == 3: - return expr / rscale - - else: - raise NotImplementedError("unsupported dimensionality") - def __getinitargs__(self): return (self.dim,) @@ -814,14 +725,6 @@ class KernelWrapper(Kernel): def get_expression(self, scaled_dist_vec): return self.inner_kernel.get_expression(scaled_dist_vec) - @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 get_derivative_transformation_at_source(self, expr_dict): return self.inner_kernel.get_derivative_transformation_at_source(expr_dict) diff --git a/sumpy/tools.py b/sumpy/tools.py index 88c5aea2924910306322b63f842ddec31b6d9635..4cc529ed1d27bb5c53332fd7829fd3ac9b5c9819 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -104,11 +104,77 @@ def mi_power(vector, mi, evaluate=True): class MiDerivativeTaker(object): - def __init__(self, expr, var_list): + def __init__(self, expr, var_list, rscale=1): + r""" + A class to take scaled derivatives of the symbolic expression + expr w.r.t. variables var_list and the scaling parameter rscale. + + Consider a Taylor multipole expansion: + + .. math:: + + 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`: + + .. math:: + + \begin{eqnarray*} + f (x) & = & g (x / \alpha),\\ + f^{(i)} (x) & = & \frac{1}{\alpha^i} g^{(i)} (x / \alpha) . + \end{eqnarray*} + + 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 + + .. math:: + + 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`. + + 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) = C \log x`, the powers of + :math:`\alpha^i` from the chain rule cancel with the ones from the + argument substituted into the kernel derivatives: + + .. 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*. + + This derivative taker returns :math:`g^{(i)}(\xi) = \alpha^i f^{(i)}` + given :math:`f^{(0)}` as *expr* and :math:`\alpha` as *rscale*. + """ + assert isinstance(expr, sym.Basic) self.var_list = var_list empty_mi = (0,) * len(var_list) self.cache_by_mi = {empty_mi: expr} + self.rscale = rscale def mi_dist(self, a, b): return np.array(a, dtype=int) - np.array(b, dtype=int) @@ -122,7 +188,7 @@ class MiDerivativeTaker(object): for next_deriv, next_mi in self.get_derivative_taking_sequence( current_mi, mi): - expr = expr.diff(next_deriv) + expr = expr.diff(next_deriv) * self.rscale self.cache_by_mi[next_mi] = expr return expr @@ -144,9 +210,10 @@ class MiDerivativeTaker(object): class LaplaceDerivativeTaker(MiDerivativeTaker): - def __init__(self, expr, var_list): - super(LaplaceDerivativeTaker, self).__init__(expr, var_list) + def __init__(self, expr, var_list, rscale=1): + super(LaplaceDerivativeTaker, self).__init__(expr, var_list, rscale) self.r = sym.sqrt(sum(v**2 for v in var_list)) + self.scaled_r = sym.sqrt(sum((v/rscale)**2 for v in var_list)) def diff(self, mi): # Return zero for negative values. Makes the algorithm readable. @@ -174,17 +241,17 @@ class LaplaceDerivativeTaker(MiDerivativeTaker): n = mi[i] if i == d: if dim == 3: - expr -= (2*n-1)*x*self.diff(tuple(mi_minus_one)) + expr -= (2*n-1)*(x/self.rscale)*self.diff(tuple(mi_minus_one)) expr -= (n-1)**2*self.diff(tuple(mi_minus_two)) else: - expr -= 2*x*(n-1)*self.diff(tuple(mi_minus_one)) + expr -= 2*(x/self.rscale)*(n-1)*self.diff(tuple(mi_minus_one)) expr -= (n-1)*(n-2)*self.diff(tuple(mi_minus_two)) if n == 2 and sum(mi) == 2: expr += 1 else: - expr -= 2*n*x*self.diff(tuple(mi_minus_one)) + expr -= 2*n*(x/self.rscale)*self.diff(tuple(mi_minus_one)) expr -= n*(n-1)*self.diff(tuple(mi_minus_two)) - expr /= self.r**2 + expr /= self.scaled_r**2 self.cache_by_mi[mi] = expr return expr diff --git a/test/test_codegen.py b/test/test_codegen.py index 02884ebade0ee850e7da9b71734f0b5a1aa5aaee..fa15592796a88dff555396e2624f38da04f31756 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, rscale=1, None) + coeffs = expn.coefficients_from_source(avec, bvec, rscale=1, sac=None) sym2pymbolic = SympyToPymbolicMapper() coeffs_pymbolic = [sym2pymbolic(c) for c in coeffs] diff --git a/test/test_kernels.py b/test/test_kernels.py index c52427c49ab1140262f847f5e383a11f962717c3..f37e51c669285d7c0500ac82f89beb30f0ad47e8 100644 --- a/test/test_kernels.py +++ b/test/test_kernels.py @@ -366,7 +366,7 @@ def _m2l_translate_simple(tgt_expansion, src_expansion, src_coeff_exprs, src_rsc # # To get the local expansion coefficients, we take derivatives of # the multipole expansion. - taker = src_expansion.get_kernel_derivative_taker(dvec) + taker = src_expansion.get_kernel_derivative_taker(dvec, src_rscale) from sumpy.tools import add_mi @@ -377,12 +377,7 @@ def _m2l_translate_simple(tgt_expansion, src_expansion, src_coeff_exprs, src_rsc src_coeff_exprs, src_expansion.get_coefficient_identifiers()): - kernel_deriv = ( - src_expansion.get_scaled_multipole( - taker.diff(add_mi(deriv, term)), - dvec, src_rscale, - nderivatives=sum(deriv) + sum(term), - nderivatives_for_scaling=sum(term))) + kernel_deriv = taker.diff(add_mi(deriv, term)) / src_rscale**sum(deriv) local_result.append( coeff * kernel_deriv * tgt_rscale**sum(deriv))