diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index b335917c3eb87b94d3eadc7bb2c30de8a82881a0..4f0a3a01482af78dacbfe0635ce45249a610f930 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -28,13 +28,12 @@ import numpy as np
import logging
from pytools import memoize_method
import sumpy.symbolic as sym
-from sumpy.tools import MiDerivativeTaker
from collections import defaultdict
__doc__ = """
.. autoclass:: ExpansionBase
-.. autoclass:: LinearRecurrenceBasedDerivativeWrangler
+.. autoclass:: LinearRecurrenceBasedExpansionTermsWrangler
Expansion Factories
^^^^^^^^^^^^^^^^^^^
@@ -159,9 +158,11 @@ class ExpansionBase(object):
# }}}
-# {{{ derivative wrangler
+# {{{ expansion terms wrangler
-class DerivativeWrangler(object):
+class ExpansionTermsWrangler(object):
+
+ init_arg_names = ("order", "dim")
def __init__(self, order, dim):
self.order = order
@@ -178,9 +179,6 @@ class DerivativeWrangler(object):
rscale):
raise NotImplementedError
- def get_derivative_taker(self, expr, var_list):
- raise NotImplementedError
-
@memoize_method
def get_full_coefficient_identifiers(self):
"""
@@ -193,18 +191,25 @@ class DerivativeWrangler(object):
res = sorted(gnitstam(self.order, self.dim), key=sum)
return res
- @memoize_method
- def copy(self, order):
- raise NotImplementedError
+ def copy(self, **kwargs):
+ new_kwargs = dict(
+ (name, getattr(self, name))
+ for name in self.init_arg_names)
+ new_kwargs["order"] = kwargs.pop("order", self.order)
+ new_kwargs["dim"] = kwargs.pop("dim", self.dim)
-class FullDerivativeWrangler(DerivativeWrangler):
+ if kwargs:
+ raise TypeError("unexpected keyword arguments '%s'"
+ % ", ".join(kwargs))
- def get_derivative_taker(self, expr, dvec):
- return MiDerivativeTaker(expr, dvec)
+ return type(self)(**new_kwargs)
+
+
+class FullExpansionTermsWrangler(ExpansionTermsWrangler):
get_coefficient_identifiers = (
- DerivativeWrangler.get_full_coefficient_identifiers)
+ ExpansionTermsWrangler.get_full_coefficient_identifiers)
def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives,
rscale):
@@ -213,12 +218,6 @@ class FullDerivativeWrangler(DerivativeWrangler):
get_stored_mpole_coefficients_from_full = (
get_full_kernel_derivatives_from_stored)
- @memoize_method
- def copy(self, **kwargs):
- order = kwargs.pop('order', self.order)
- dim = kwargs.pop('dim', self.dim)
- return type(self)(order, dim)
-
# {{{ sparse matrix-vector multiplication
@@ -255,24 +254,27 @@ def _spmv(spmat, x, sparse_vectors):
# }}}
-class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
+class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler):
"""
.. automethod:: __init__
.. automethod:: get_pde_dict
.. automethod:: get_reduced_coeffs
"""
+ init_arg_names = ("order", "dim", "deriv_multiplier")
+
def __init__(self, order, dim, deriv_multiplier):
- """
+ r"""
:param order: order of the expansion
:param dim: number of dimensions
- :param deriv_multiplier: symbolic expression that should be multiplied
- with each coefficient_identifier to recover the linear
- recurrence for the PDE.
-
- .. seealso:: :func:`~LinearRecurrenceBasedDerivativeWrangler.get_pde_dict`
+ :param deriv_multiplier: a symbolic expression that is used to
+ 'normalize out' constant coefficients in the PDE in
+ :func:`~LinearRecurrenceBasedExpansionTermsWrangler.get_pde_dict`, so
+ that the Taylor coefficient with multi-index :math:`\nu` as seen by
+ that representation of the PDE is :math:`\text{coeff} /
+ {\text{deriv\_multiplier}^{|\nu|}}`.
"""
- super(LinearRecurrenceBasedDerivativeWrangler, self).__init__(order, dim)
+ super(LinearRecurrenceBasedExpansionTermsWrangler, self).__init__(order, dim)
self.deriv_multiplier = deriv_multiplier
def get_coefficient_identifiers(self):
@@ -306,19 +308,19 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
Return a matrix that expresses every derivative in terms of a
set of "stored" derivatives.
- For example, for the recurrence
+ For example, for the recurrence::
u_xx + u_yy + u_zz = 0
- the coefficient matrix features the following entries:
+ the coefficient matrix features the following entries::
- ... u_xx u_yy ... <= cols = only stored derivatives
- ==================
- ...| ... ... ... ...
- |
- u_zz| ... -1 -1 ...
+ ... u_xx u_yy ... <= cols = only stored derivatives
+ ==================
+ ...| ... ... ... ...
+ |
+ u_zz| ... -1 -1 ...
- ^ rows = one for every derivative
+ ^ rows = one for every derivative
"""
stored_identifiers, coeff_matrix = self._get_stored_ids_and_coeff_mat()
@@ -347,9 +349,8 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
tol = 1e-13
stored_identifiers = []
- import time
- start_time = time.time()
- logger.debug("computing recurrence for Taylor coefficients: start")
+ from pytools import ProcessLogger
+ plog = ProcessLogger(logger, "compute recurrence for Taylor coefficients")
pde_dict = self.get_pde_dict()
pde_mat = []
@@ -371,13 +372,18 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
else:
pde_mat.append(eq)
if len(pde_mat) > 0:
- """
+ r"""
Find a matrix :math:`s` such that :math:`K = S^T K_{[r]}`
where :math:`r` is the indices of the reduced coefficients and
:math:`K` is a column vector of coefficients. Let :math:`P` be the
- PDE matrix and :math:`N` be the nullspace of :math:`P`.
+ PDE matrix, i.e. the matrix obtained by applying the PDE
+ as an identity to each of the Taylor coefficients.
+ (Realize that those, as functions of space, each still satisfy the PDE.)
+ As a result, if :math:`\mathbf\alpha` is a vector of Taylor coefficients,
+ one would expect :math:`P\mathbf\alpha = \mathbf 0`.
+ Further, let :math:`N` be the nullspace of :math:`P`.
Then, :math:`S^T = N * N_{[r, :]}^{-1}` which implies,
- :math:`S = N_{[r, :]}^{-T} N^T = solve(N_{[r, :]}^T, N^T)`
+ :math:`S = N_{[r, :]}^{-T} N^T = N_{[r, :]}^{-T} N^T`.
"""
pde_mat = np.array(pde_mat, dtype=np.float64)
n = nullspace(pde_mat, atol=tol)
@@ -397,9 +403,7 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
if np.abs(s[i][j]) > tol:
coeff_matrix[j].append((i, s[i][j]))
- logger.debug("computing recurrence for Taylor coefficients: "
- "done after {dur:.2f} seconds"
- .format(dur=time.time() - start_time))
+ plog.done()
logger.debug("number of Taylor coefficients was reduced from {orig} to {red}"
.format(orig=len(self.get_full_coefficient_identifiers()),
@@ -410,14 +414,20 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
def get_pde_dict(self):
r"""
Returns the PDE as a dictionary of (mi, coeff) such that mi
- is the multi index of the derivative and the PDE is represented by,
+ is the multi-index of the derivative and the PDE is represented by,
.. math::
- \sum(\frac{mi \times coeff}{deriv\_multiplier^{sum(mi)}}) = 0
+ \sum_{\nu,c_\nu\in \text{pde\_dict}}
+ \frac{c_\nu\cdot \alpha_\nu}
+ {\text{deriv\_multiplier}^{
+ \sum \text{mi}
+ }} = 0,
+
+ where :math:`\mathbf\alpha` is a coefficient vector.
- Note that coeff should be numeric instead of symbolic to enable use of
- numeric linear algebra routines. `deriv_multiplier` can be symbolic
+ Note that *coeff* should be a number (not symbolic) to enable use of
+ numeric linear algebra routines. *deriv_multiplier* can be symbolic
and should be used when the PDE has symbolic values as coefficients
for the derivatives.
@@ -440,25 +450,13 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
"""
raise NotImplementedError
- def get_derivative_taker(self, expr, var_list):
- from sumpy.tools import MiDerivativeTaker
- return MiDerivativeTaker(expr, var_list)
-
- @memoize_method
- def copy(self, **kwargs):
- obj = type(self).__new__(type(self))
- order = kwargs.pop('order', self.order)
- dim = kwargs.pop('dim', self.dim)
- deriv_multiplier = kwargs.pop('deriv_multiplier', self.deriv_multiplier)
- LinearRecurrenceBasedDerivativeWrangler.__init__(obj, order,
- dim, deriv_multiplier)
- return obj
+class LaplaceExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler):
-class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
+ init_arg_names = ("order", "dim")
def __init__(self, order, dim):
- super(LaplaceDerivativeWrangler, self).__init__(order, dim, 1)
+ super(LaplaceExpansionTermsWrangler, self).__init__(order, dim, 1)
def get_pde_dict(self):
pde_dict = {}
@@ -478,11 +476,15 @@ class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
return idx
-class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
+class HelmholtzExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler):
+
+ init_arg_names = ("order", "dim", "helmholtz_k_name")
def __init__(self, order, dim, helmholtz_k_name):
+ self.helmholtz_k_name = helmholtz_k_name
+
multiplier = sym.Symbol(helmholtz_k_name)
- super(HelmholtzDerivativeWrangler, self).__init__(order, dim, multiplier)
+ super(HelmholtzExpansionTermsWrangler, self).__init__(order, dim, multiplier)
def get_pde_dict(self, **kwargs):
pde_dict = {}
@@ -510,31 +512,33 @@ class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
class VolumeTaylorExpansionBase(object):
@classmethod
- def get_or_make_derivative_wrangler(cls, *key):
+ def get_or_make_expansion_terms_wrangler(cls, *key):
"""
- This stores the derivative wrangler at the class attribute level because
- precomputing the derivative wrangler may become expensive.
+ This stores the expansion terms wrangler at the class attribute level
+ because recreating the expansion terms wrangler implicitly empties its
+ caches.
"""
try:
- wrangler = cls.derivative_wrangler_cache[key]
+ wrangler = cls.expansion_terms_wrangler_cache[key]
except KeyError:
- wrangler = cls.derivative_wrangler_class(*key)
- cls.derivative_wrangler_cache[key] = wrangler
+ wrangler = cls.expansion_terms_wrangler_class(*key)
+ cls.expansion_terms_wrangler_cache[key] = wrangler
return wrangler
@property
- def derivative_wrangler(self):
- return self.get_or_make_derivative_wrangler(*self.derivative_wrangler_key)
+ def expansion_terms_wrangler(self):
+ return self.get_or_make_expansion_terms_wrangler(
+ *self.expansion_terms_wrangler_key)
def get_coefficient_identifiers(self):
"""
Returns the identifiers of the coefficients that actually get stored.
"""
- return self.derivative_wrangler.get_coefficient_identifiers()
+ return self.expansion_terms_wrangler.get_coefficient_identifiers()
def get_full_coefficient_identifiers(self):
- return self.derivative_wrangler.get_full_coefficient_identifiers()
+ return self.expansion_terms_wrangler.get_full_coefficient_identifiers()
@property
@memoize_method
@@ -548,33 +552,33 @@ class VolumeTaylorExpansionBase(object):
class VolumeTaylorExpansion(VolumeTaylorExpansionBase):
- derivative_wrangler_class = FullDerivativeWrangler
- derivative_wrangler_cache = {}
+ expansion_terms_wrangler_class = FullExpansionTermsWrangler
+ expansion_terms_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)
+ self.expansion_terms_wrangler_key = (order, kernel.dim)
class LaplaceConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase):
- derivative_wrangler_class = LaplaceDerivativeWrangler
- derivative_wrangler_cache = {}
+ expansion_terms_wrangler_class = LaplaceExpansionTermsWrangler
+ expansion_terms_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)
+ self.expansion_terms_wrangler_key = (order, kernel.dim)
class HelmholtzConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase):
- derivative_wrangler_class = HelmholtzDerivativeWrangler
- derivative_wrangler_cache = {}
+ expansion_terms_wrangler_class = HelmholtzExpansionTermsWrangler
+ expansion_terms_wrangler_cache = {}
# 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)
+ self.expansion_terms_wrangler_key = (order, kernel.dim, helmholtz_k_name)
# }}}
diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py
index b9c5d0a673b1c7be2a42509a01a32d5ceb1fc7d3..11b1368a43bd8c57ad51065f84bea3b7dc2a5ca5 100644
--- a/sumpy/expansion/local.py
+++ b/sumpy/expansion/local.py
@@ -128,7 +128,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
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(
+ self.expansion_terms_wrangler.get_full_kernel_derivatives_from_stored(
coeffs, rscale))
bvec = bvec * rscale**-1
result = sum(
@@ -167,48 +167,55 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
from sumpy.tools import add_mi
- src_max_sum = max(sum(mi) for mi in
- src_expansion.get_coefficient_identifiers())
- tgt_max_sum = max(sum(mi) for mi in
- self.get_coefficient_identifiers())
+ # Create a expansion terms wrangler for derivatives up to order
+ # (tgt order)+(src order).
+ tgtplusderiv_exp_terms_wrangler = \
+ src_expansion.expansion_terms_wrangler.copy(
+ order=self.order + src_expansion.order)
+ tgtplusderiv_coeff_ids = \
+ tgtplusderiv_exp_terms_wrangler.get_coefficient_identifiers()
+ tgtplusderiv_full_coeff_ids = \
+ tgtplusderiv_exp_terms_wrangler.get_full_coefficient_identifiers()
- max_sum = src_max_sum + tgt_max_sum
- new_deriv_wrangler = \
- src_expansion.derivative_wrangler.copy(order=max_sum)
- new_coeffs = new_deriv_wrangler.get_coefficient_identifiers()
- new_full_coeffs = new_deriv_wrangler.get_full_coefficient_identifiers()
-
- ident_to_index = dict((ident, i) for i, ident in
- enumerate(new_full_coeffs))
+ tgtplusderiv_ident_to_index = dict((ident, i) for i, ident in
+ enumerate(tgtplusderiv_full_coeff_ids))
result = []
- for deriv in self.get_coefficient_identifiers():
- local_result = []
+ for lexp_mi in self.get_coefficient_identifiers():
+ lexp_mi_terms = []
- full_coeffs = [0] * len(new_full_coeffs)
+ # Embed the source coefficients in the coefficient array
+ # for the (tgt order)+(src order) wrangler, offset by lexp_mi.
+ embedded_coeffs = [0] * len(tgtplusderiv_full_coeff_ids)
for coeff, term in zip(
src_coeff_exprs,
src_expansion.get_coefficient_identifiers()):
- full_coeffs[ident_to_index[add_mi(deriv, term)]] = coeff
+ embedded_coeffs[
+ tgtplusderiv_ident_to_index[add_mi(lexp_mi, term)]] \
+ = coeff
- stored_coeffs = \
- new_deriv_wrangler.get_stored_mpole_coefficients_from_full(
- full_coeffs, src_rscale)
+ # Compress the embedded coefficient set
+ stored_coeffs = tgtplusderiv_exp_terms_wrangler \
+ .get_stored_mpole_coefficients_from_full(
+ embedded_coeffs, src_rscale)
+ # Sum the embedded coefficient set
for i, coeff in enumerate(stored_coeffs):
if coeff == 0:
continue
- nderivatives_for_scaling = sum(new_coeffs[i])-sum(deriv)
+ nderivatives_for_scaling = \
+ sum(tgtplusderiv_coeff_ids[i])-sum(lexp_mi)
kernel_deriv = (
src_expansion.get_scaled_multipole(
- taker.diff(new_coeffs[i]),
+ taker.diff(tgtplusderiv_coeff_ids[i]),
dvec, src_rscale,
- nderivatives=sum(new_coeffs[i]),
+ nderivatives=sum(tgtplusderiv_coeff_ids[i]),
nderivatives_for_scaling=nderivatives_for_scaling))
- local_result.append(
- coeff * kernel_deriv * tgt_rscale**sum(deriv))
- result.append(sym.Add(*local_result))
+ lexp_mi_terms.append(
+ coeff * kernel_deriv * tgt_rscale**sum(lexp_mi))
+ result.append(sym.Add(*lexp_mi_terms))
+
else:
from sumpy.tools import MiDerivativeTaker
expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=src_rscale)
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index 92ae80f3308b4bf8820fbc00da9229725058e3e1..25488952e4b7425371ecf23aefa261c0c70900f5 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -94,7 +94,7 @@ 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(
+ self.expansion_terms_wrangler.get_stored_mpole_coefficients_from_full(
result, rscale))
def get_scaled_multipole(self, expr, bvec, rscale, nderivatives,
@@ -127,8 +127,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
return result
def get_kernel_derivative_taker(self, bvec):
- return (self.derivative_wrangler.get_derivative_taker(
- self.kernel.get_expression(bvec), bvec))
+ from sumpy.tools import MiDerivativeTaker
+ return MiDerivativeTaker(self.kernel.get_expression(bvec), bvec)
def translate_from(self, src_expansion, src_coeff_exprs, src_rscale,
dvec, tgt_rscale):
@@ -188,7 +188,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
logger.info("building translation operator: done")
return (
- self.derivative_wrangler.get_stored_mpole_coefficients_from_full(
+ self.expansion_terms_wrangler.get_stored_mpole_coefficients_from_full(
result, tgt_rscale))