From 71571ba877f1a05d43c1280767a0bb52c38d2872 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Mon, 30 Apr 2018 16:47:57 -0400 Subject: [PATCH 1/9] Various readability/doc/convention improvements --- sumpy/expansion/__init__.py | 75 +++++++++++++++++++------------------ sumpy/expansion/local.py | 3 ++ 2 files changed, 41 insertions(+), 37 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index b335917c..7ab2125c 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -217,6 +217,11 @@ class FullDerivativeWrangler(DerivativeWrangler): def copy(self, **kwargs): order = kwargs.pop('order', self.order) dim = kwargs.pop('dim', self.dim) + + if kwargs: + raise TypeError("unexpected keyword arguments '%s'" + % ", ".join(kwargs)) + return type(self)(order, dim) @@ -263,14 +268,14 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): """ 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:`~LinearRecurrenceBasedDerivativeWrangler.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) self.deriv_multiplier = deriv_multiplier @@ -306,19 +311,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 +352,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 +375,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 +406,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 +417,18 @@ 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, - Note that coeff should be numeric instead of symbolic to enable use of - numeric linear algebra routines. `deriv_multiplier` can be symbolic + where :math:`\mathbf\alpha` is a coefficient vector. + + 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. @@ -444,16 +455,6 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): 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 LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index b9c5d0a6..0336ef43 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -169,8 +169,10 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): src_max_sum = max(sum(mi) for mi in src_expansion.get_coefficient_identifiers()) + assert src_max_sum == src_expansion.order tgt_max_sum = max(sum(mi) for mi in self.get_coefficient_identifiers()) + assert tgt_max_sum == tgt_expansion.order max_sum = src_max_sum + tgt_max_sum new_deriv_wrangler = \ @@ -209,6 +211,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): local_result.append( coeff * kernel_deriv * tgt_rscale**sum(deriv)) result.append(sym.Add(*local_result)) + else: from sumpy.tools import MiDerivativeTaker expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=src_rscale) -- GitLab From edd50e7336e30474e31cfa1b6f289f2c83aace1b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 30 Apr 2018 18:09:54 -0500 Subject: [PATCH 2/9] Rename derivative wrangler -> expansion terms wrangler, improve understandability of wrangled M2L --- sumpy/expansion/__init__.py | 105 +++++++++++++++-------------------- sumpy/expansion/local.py | 58 +++++++++---------- sumpy/expansion/multipole.py | 8 +-- 3 files changed, 80 insertions(+), 91 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 7ab2125c..22223b47 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,9 @@ class ExpansionBase(object): # }}} -# {{{ derivative wrangler +# {{{ expansion terms wrangler -class DerivativeWrangler(object): +class ExpansionTermsWrangler(object): def __init__(self, order, dim): self.order = order @@ -178,9 +177,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 +189,21 @@ 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): + order = kwargs.pop('order', self.order) + dim = kwargs.pop('dim', self.dim) + + if kwargs: + raise TypeError("unexpected keyword arguments '%s'" + % ", ".join(kwargs)) + return type(self)(order, dim) -class FullDerivativeWrangler(DerivativeWrangler): - def get_derivative_taker(self, expr, dvec): - return MiDerivativeTaker(expr, dvec) +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,17 +212,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) - - if kwargs: - raise TypeError("unexpected keyword arguments '%s'" - % ", ".join(kwargs)) - - return type(self)(order, dim) - # {{{ sparse matrix-vector multiplication @@ -260,7 +248,7 @@ def _spmv(spmat, x, sparse_vectors): # }}} -class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): +class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler): """ .. automethod:: __init__ .. automethod:: get_pde_dict @@ -271,13 +259,14 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): r""" :param order: order of the expansion :param dim: number of dimensions - :param deriv_multiplier: a symbolic expression that is used to 'normalize out' - constant coefficients in the PDE in - :func:`~LinearRecurrenceBasedDerivativeWrangler.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|}}`. + :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): @@ -451,15 +440,11 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler): """ raise NotImplementedError - def get_derivative_taker(self, expr, var_list): - from sumpy.tools import MiDerivativeTaker - return MiDerivativeTaker(expr, var_list) - -class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): +class LaplaceExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler): 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 = {} @@ -479,11 +464,11 @@ class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): return idx -class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): +class HelmholtzExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler): def __init__(self, order, dim, 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 = {} @@ -511,31 +496,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 @@ -549,33 +536,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 0336ef43..716fa93f 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,50 +167,52 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): from sumpy.tools import add_mi - src_max_sum = max(sum(mi) for mi in - src_expansion.get_coefficient_identifiers()) - assert src_max_sum == src_expansion.order - tgt_max_sum = max(sum(mi) for mi in - self.get_coefficient_identifiers()) - assert tgt_max_sum == tgt_expansion.order + # 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) + stored_coeffs = tgtplusderiv_exp_terms_wrangler \ + .get_stored_mpole_coefficients_from_full( + embedded_coeffs, src_rscale) 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 diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 92ae80f3..25488952 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)) -- GitLab From 11b7d29b99367eb4a0e2ae1fa73ce04c183e75a8 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 30 Apr 2018 18:27:34 -0500 Subject: [PATCH 3/9] More M2L summary comments --- sumpy/expansion/local.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 716fa93f..2f494a66 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -194,10 +194,12 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): tgtplusderiv_ident_to_index[add_mi(lexp_mi, term)]] \ = coeff + # Compress the embedded coefficient set stored_coeffs = tgtplusderiv_exp_terms_wrangler \ .get_stored_mpole_coefficients_from_full( embedded_coeffs, src_rscale) + # Sum the embeeded coefficient set for i, coeff in enumerate(stored_coeffs): if coeff == 0: continue -- GitLab From f4394132f3324c1ff3d8727ee6ad1a53c63a87bd Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 30 Apr 2018 19:09:48 -0500 Subject: [PATCH 4/9] Do away with embedding construction in M2L --- sumpy/expansion/local.py | 70 ++++++++++++++-------------------------- 1 file changed, 25 insertions(+), 45 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 2f494a66..ad29ea05 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -157,64 +157,44 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # # coeff0 * diff(kernel, mi0) + coeff1 * diff(kernel, mi1) + ... # - # To get the local expansion coefficients, we take derivatives of - # the multipole expansion. + # The local is formed by differentiating the mpole, i.e. the + # target variable is in the mpole basis--so that is what the + # derivatives for the local will hit. # - # This code speeds up derivative taking by caching all kernel - # derivatives. + # We group the calculation here so that we take all derivatives + # in one big gulp, which ends up being more efficient than + # evaluating the mpole (as a giant sum) and then differentiating + # that. + # + # In particular, this speeds up derivative taking by allowing all + # kernel derivatives to be cached. taker = src_expansion.get_kernel_derivative_taker(dvec) from sumpy.tools import add_mi - # 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() - - tgtplusderiv_ident_to_index = dict((ident, i) for i, ident in - enumerate(tgtplusderiv_full_coeff_ids)) + src_exp_terms_wrangler = src_expansion.expansion_terms_wrangler + src_exp_stored_coeff_ids = \ + src_exp_terms_wrangler.get_coefficient_identifiers() result = [] - for lexp_mi in self.get_coefficient_identifiers(): - lexp_mi_terms = [] - - # 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()): - embedded_coeffs[ - tgtplusderiv_ident_to_index[add_mi(lexp_mi, term)]] \ - = coeff - - # Compress the embedded coefficient set - stored_coeffs = tgtplusderiv_exp_terms_wrangler \ - .get_stored_mpole_coefficients_from_full( - embedded_coeffs, src_rscale) - - # Sum the embeeded coefficient set - for i, coeff in enumerate(stored_coeffs): - if coeff == 0: - continue - nderivatives_for_scaling = \ - sum(tgtplusderiv_coeff_ids[i])-sum(lexp_mi) + for tgt_coeff_id in self.get_coefficient_identifiers(): + tgt_coeff_terms = [] + + for src_coeff_id, coeff in zip( + src_exp_stored_coeff_ids, src_coeff_exprs): + nderivatives_for_scaling = sum(src_coeff_id) + term_mi = add_mi(src_coeff_id, tgt_coeff_id) kernel_deriv = ( src_expansion.get_scaled_multipole( - taker.diff(tgtplusderiv_coeff_ids[i]), + taker.diff(term_mi), dvec, src_rscale, - nderivatives=sum(tgtplusderiv_coeff_ids[i]), + nderivatives=sum(term_mi), nderivatives_for_scaling=nderivatives_for_scaling)) - lexp_mi_terms.append( - coeff * kernel_deriv * tgt_rscale**sum(lexp_mi)) - result.append(sym.Add(*lexp_mi_terms)) + tgt_coeff_terms.append( + coeff * kernel_deriv * tgt_rscale**sum(tgt_coeff_id)) + result.append(sym.Add(*tgt_coeff_terms)) else: from sumpy.tools import MiDerivativeTaker -- GitLab From e6729cedfe0910f9a5341e5dcd7269319743c7cf Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 30 Apr 2018 19:14:53 -0500 Subject: [PATCH 5/9] Remove ExpansionTermsWrangler.copy, fix docstring formatting --- sumpy/expansion/__init__.py | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 22223b47..5b6ac7c2 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -189,16 +189,6 @@ class ExpansionTermsWrangler(object): res = sorted(gnitstam(self.order, self.dim), key=sum) return res - def copy(self, **kwargs): - order = kwargs.pop('order', self.order) - dim = kwargs.pop('dim', self.dim) - - if kwargs: - raise TypeError("unexpected keyword arguments '%s'" - % ", ".join(kwargs)) - - return type(self)(order, dim) - class FullExpansionTermsWrangler(ExpansionTermsWrangler): @@ -261,9 +251,9 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler): :param dim: number of dimensions :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} / + :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(LinearRecurrenceBasedExpansionTermsWrangler, self).__init__(order, dim) -- GitLab From b2400584b2908475c4502569d1fe4418d646acc1 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 30 Apr 2018 22:59:36 -0500 Subject: [PATCH 6/9] Revert "Do away with embedding construction in M2L" This reverts commit f4394132f3324c1ff3d8727ee6ad1a53c63a87bd. --- sumpy/expansion/local.py | 70 ++++++++++++++++++++++++++-------------- 1 file changed, 45 insertions(+), 25 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index ad29ea05..2f494a66 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -157,44 +157,64 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # # coeff0 * diff(kernel, mi0) + coeff1 * diff(kernel, mi1) + ... # - # The local is formed by differentiating the mpole, i.e. the - # target variable is in the mpole basis--so that is what the - # derivatives for the local will hit. + # To get the local expansion coefficients, we take derivatives of + # the multipole expansion. # - # We group the calculation here so that we take all derivatives - # in one big gulp, which ends up being more efficient than - # evaluating the mpole (as a giant sum) and then differentiating - # that. - # - # In particular, this speeds up derivative taking by allowing all - # kernel derivatives to be cached. + # 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 - src_exp_terms_wrangler = src_expansion.expansion_terms_wrangler - src_exp_stored_coeff_ids = \ - src_exp_terms_wrangler.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() - result = [] - for tgt_coeff_id in self.get_coefficient_identifiers(): - tgt_coeff_terms = [] + tgtplusderiv_ident_to_index = dict((ident, i) for i, ident in + enumerate(tgtplusderiv_full_coeff_ids)) - for src_coeff_id, coeff in zip( - src_exp_stored_coeff_ids, src_coeff_exprs): - nderivatives_for_scaling = sum(src_coeff_id) - term_mi = add_mi(src_coeff_id, tgt_coeff_id) + result = [] + for lexp_mi in self.get_coefficient_identifiers(): + lexp_mi_terms = [] + + # 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()): + embedded_coeffs[ + tgtplusderiv_ident_to_index[add_mi(lexp_mi, term)]] \ + = coeff + + # Compress the embedded coefficient set + stored_coeffs = tgtplusderiv_exp_terms_wrangler \ + .get_stored_mpole_coefficients_from_full( + embedded_coeffs, src_rscale) + + # Sum the embeeded coefficient set + for i, coeff in enumerate(stored_coeffs): + if coeff == 0: + continue + nderivatives_for_scaling = \ + sum(tgtplusderiv_coeff_ids[i])-sum(lexp_mi) kernel_deriv = ( src_expansion.get_scaled_multipole( - taker.diff(term_mi), + taker.diff(tgtplusderiv_coeff_ids[i]), dvec, src_rscale, - nderivatives=sum(term_mi), + nderivatives=sum(tgtplusderiv_coeff_ids[i]), nderivatives_for_scaling=nderivatives_for_scaling)) - tgt_coeff_terms.append( - coeff * kernel_deriv * tgt_rscale**sum(tgt_coeff_id)) - result.append(sym.Add(*tgt_coeff_terms)) + 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 -- GitLab From 7e003fb3f22c72c025b98d1c4690009dfca85006 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 30 Apr 2018 22:59:39 -0500 Subject: [PATCH 7/9] Revert "Remove ExpansionTermsWrangler.copy, fix docstring formatting" This reverts commit e6729cedfe0910f9a5341e5dcd7269319743c7cf. --- sumpy/expansion/__init__.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 5b6ac7c2..22223b47 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -189,6 +189,16 @@ class ExpansionTermsWrangler(object): res = sorted(gnitstam(self.order, self.dim), key=sum) return res + def copy(self, **kwargs): + order = kwargs.pop('order', self.order) + dim = kwargs.pop('dim', self.dim) + + if kwargs: + raise TypeError("unexpected keyword arguments '%s'" + % ", ".join(kwargs)) + + return type(self)(order, dim) + class FullExpansionTermsWrangler(ExpansionTermsWrangler): @@ -251,9 +261,9 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler): :param dim: number of dimensions :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} / + :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(LinearRecurrenceBasedExpansionTermsWrangler, self).__init__(order, dim) -- GitLab From b471e3a4dcef9319b649f57203c3b19c69843390 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 30 Apr 2018 23:00:23 -0500 Subject: [PATCH 8/9] Comment/docstring fixes --- sumpy/expansion/__init__.py | 6 +++--- sumpy/expansion/local.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 22223b47..413d35ab 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -261,9 +261,9 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler): :param dim: number of dimensions :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} / + :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(LinearRecurrenceBasedExpansionTermsWrangler, self).__init__(order, dim) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 2f494a66..11b1368a 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -199,7 +199,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): .get_stored_mpole_coefficients_from_full( embedded_coeffs, src_rscale) - # Sum the embeeded coefficient set + # Sum the embedded coefficient set for i, coeff in enumerate(stored_coeffs): if coeff == 0: continue -- GitLab From 5988025708b67ab3b95cc649d85b6ef25de9c66c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 30 Apr 2018 23:14:00 -0500 Subject: [PATCH 9/9] Untangle ExpansionTermsWrangler.copy implementation --- sumpy/expansion/__init__.py | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 413d35ab..4f0a3a01 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -162,6 +162,8 @@ class ExpansionBase(object): class ExpansionTermsWrangler(object): + init_arg_names = ("order", "dim") + def __init__(self, order, dim): self.order = order self.dim = dim @@ -190,14 +192,18 @@ class ExpansionTermsWrangler(object): return res def copy(self, **kwargs): - order = kwargs.pop('order', self.order) - dim = kwargs.pop('dim', self.dim) + 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) if kwargs: raise TypeError("unexpected keyword arguments '%s'" % ", ".join(kwargs)) - return type(self)(order, dim) + return type(self)(**new_kwargs) class FullExpansionTermsWrangler(ExpansionTermsWrangler): @@ -255,6 +261,8 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler): .. 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 @@ -412,7 +420,9 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler): \sum_{\nu,c_\nu\in \text{pde\_dict}} \frac{c_\nu\cdot \alpha_\nu} - {\text{deriv\_multiplier}^{\sum \text{mi}}} = 0, + {\text{deriv\_multiplier}^{ + \sum \text{mi} + }} = 0, where :math:`\mathbf\alpha` is a coefficient vector. @@ -443,6 +453,8 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler): class LaplaceExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler): + init_arg_names = ("order", "dim") + def __init__(self, order, dim): super(LaplaceExpansionTermsWrangler, self).__init__(order, dim, 1) @@ -466,7 +478,11 @@ class LaplaceExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler) 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(HelmholtzExpansionTermsWrangler, self).__init__(order, dim, multiplier) -- GitLab