diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index aa8ae80adf6740131933528c34ecb98d8d7fb0c3..2264e926055b28bf7bde2c3ddaa54f38bed42221 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -12,8 +12,7 @@ dependencies: - islpy - pyopencl - python=3 -- symengine=0.3.0 -- python-symengine=0.3.0 +- python-symengine=0.6.0 - pyfmmlib - pip diff --git a/asv.conf.json b/asv.conf.json index 9632c4fb70310c528b43e3d055a96d96875e629f..6d6fd675c48cbdf02e0bc8448f22890dde2a3e26 100644 --- a/asv.conf.json +++ b/asv.conf.json @@ -68,7 +68,7 @@ // }, "matrix": { "numpy" : [""], - "sympy" : ["1.0"], + "sympy" : ["1.4"], "pyopencl" : [""], "islpy" : [""], "pocl" : [""], diff --git a/sumpy/e2p.py b/sumpy/e2p.py index 408b736408d9247db72d56da75e27e0b7dd6c857..7b0072ad55708ba205ceed3448914ad0d8eda281 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -91,7 +91,7 @@ class E2PBase(KernelCacheWrapper): coeff_exprs = [sym.Symbol("coeff%d" % i) for i in range(len(self.expansion.get_coefficient_identifiers()))] - value = self.expansion.evaluate(coeff_exprs, bvec, rscale, sac=sac) + value = self.expansion.evaluate(coeff_exprs, bvec, rscale) result_names = [ sac.assign_unique("result_%d_p" % i, diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 85c93850795318acd8abda04755b93b343cddcbc..9d3a917a9a16516bf6968a0446eccc0886554ed2 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -114,21 +114,19 @@ class ExpansionBase(object): """ raise NotImplementedError - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + def coefficients_from_source(self, avec, bvec, rscale): """Form an expansion from a source point. :arg avec: vector from source to center. :arg bvec: vector from center to target. Not usually necessary, except for line-Taylor expansion. - :arg sac: a symbolic assignment collection where temporary - expressions are stored. :returns: a list of :mod:`sympy` expressions representing the coefficients of the expansion. """ raise NotImplementedError - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): """ :return: a :mod:`sympy` expression corresponding to the evaluated expansion with the coefficients @@ -175,11 +173,11 @@ class ExpansionTermsWrangler(object): raise NotImplementedError def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives, - rscale, sac=None): + rscale): raise NotImplementedError def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients, - rscale, sac=None): + rscale): raise NotImplementedError @memoize_method @@ -218,6 +216,48 @@ class ExpansionTermsWrangler(object): return type(self)(**new_kwargs) + @memoize_method + def _get_coeff_identifier_split(self): + """ + This splits the coefficients into a O(p) number of disjoint sets + so that for each set, all the identifiers have the form, + (m_1, m_2, ..., m_{j-1}, c, m_{j+1}, ... , m_d) + where c is a constant. + + If this is an instance of LinearPDEBasedExpansionTermsWrangler, + then the number of sets will be O(1). + + Returns a List[Tuple[j, List[identifiers]]] + """ + res = [] + mis = self.get_full_coefficient_identifiers() + coeff_ident_enumerate_dict = dict((tuple(mi), i) for + (i, mi) in enumerate(mis)) + + max_mi = None + if isinstance(self, LinearPDEBasedExpansionTermsWrangler): + pde_dict = self.get_pde().eq + for ident in pde_dict.keys(): + if ident not in coeff_ident_enumerate_dict: + break + else: + max_mi_idx = max(coeff_ident_enumerate_dict[ident] for + ident in pde_dict.keys()) + max_mi = mis[max_mi_idx] + + if max_mi is None: + max_mi = [max(mi[d] for mi in mis) for d in range(self.dim)] + + seen_mis = set() + for d in range(self.dim): + filtered_mis = [mi for mi in mis if mi[d] < max_mi[d]] + for i in range(max_mi[d]): + new_mis = [mi for mi in filtered_mis if mi[d] == i + and mi not in seen_mis] + seen_mis.update(new_mis) + res.append((d, new_mis)) + return res + class FullExpansionTermsWrangler(ExpansionTermsWrangler): @@ -225,7 +265,7 @@ class FullExpansionTermsWrangler(ExpansionTermsWrangler): ExpansionTermsWrangler.get_full_coefficient_identifiers) def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives, - rscale, sac=None): + rscale): return stored_kernel_derivatives get_stored_mpole_coefficients_from_full = ( @@ -267,7 +307,7 @@ def _spmv(spmat, x, sparse_vectors): # }}} -def _fast_spmv(reconstruct_matrix, vec, sac, transpose=False): +def _fast_spmv(reconstruct_matrix, vec, transpose=False): if not transpose: res = [0] * len(reconstruct_matrix) stored_idx = 0 @@ -278,8 +318,7 @@ def _fast_spmv(reconstruct_matrix, vec, sac, transpose=False): else: for k, v in deps: res[row] += res[k] * v - new_sym = sym.Symbol(sac.assign_unique("expr", res[row])) - res[row] = new_sym + res[row] = sym.UnevaluatedExpr(res[row]) return res else: res = [] @@ -288,7 +327,7 @@ def _fast_spmv(reconstruct_matrix, vec, sac, transpose=False): if len(deps) == 0: res.append(expr_all[row]) continue - new_sym = sym.Symbol(sac.assign_unique("expr", expr_all[row])) + new_sym = sym.UnevaluatedExpr(expr_all[row]) for k, v in deps: expr_all[k] += new_sym * v res.reverse() @@ -316,26 +355,26 @@ class LinearPDEBasedExpansionTermsWrangler(ExpansionTermsWrangler): return self.stored_identifiers def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives, - rscale, sac=None): + rscale): coeff_matrix, reconstruct_matrix, use_reconstruct = \ self.get_coefficient_matrix(rscale) - if not use_reconstruct or sac is None: + if not use_reconstruct: return _spmv(coeff_matrix, stored_kernel_derivatives, sparse_vectors=False) - return _fast_spmv(reconstruct_matrix, stored_kernel_derivatives, sac) + return _fast_spmv(reconstruct_matrix, stored_kernel_derivatives) def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients, - rscale, sac=None): + rscale): # = M^T x, where M = coeff matrix coeff_matrix, reconstruct_matrix, use_reconstruct = \ self.get_coefficient_matrix(rscale) - if not use_reconstruct or sac is None: + if not use_reconstruct: result = [0] * len(self.stored_identifiers) for row, coeff in enumerate(full_mpole_coefficients): for col, val in coeff_matrix[row]: result[col] += coeff * val return result - return _fast_spmv(reconstruct_matrix, full_mpole_coefficients, sac, + return _fast_spmv(reconstruct_matrix, full_mpole_coefficients, transpose=True) @property diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index a58c598173a352afaa92d6a6d06855ca07275457..c32778ae14c525c2a7acdffbc2850ead0a160691 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -58,7 +58,7 @@ class LineTaylorLocalExpansion(LocalExpansionBase): def get_coefficient_identifiers(self): return list(range(self.order+1)) - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + 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 " @@ -99,7 +99,7 @@ class LineTaylorLocalExpansion(LocalExpansionBase): .subs("tau", 0) for i in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): # no point in heeding rscale here--just ignore it from pytools import factorial return sym.Add(*( @@ -116,7 +116,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): Coefficients represent derivative values of the kernel. """ - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + 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) @@ -126,91 +126,20 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): taker.diff(mi) * rscale ** sum(mi) for mi in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec, rscale, sac=None): - from pytools import factorial - + def evaluate(self, coeffs, bvec, rscale): evaluated_coeffs = ( self.expansion_terms_wrangler.get_full_kernel_derivatives_from_stored( coeffs, rscale)) bvec = [b*rscale**-1 for b in bvec] - mi_to_index = dict((mi, i) for i, mi in - enumerate(self.get_full_coefficient_identifiers())) - - # Sort multi-indices so that last dimension varies fastest - sorted_target_mis = sorted(self.get_full_coefficient_identifiers()) - dim = self.dim - - # Start with an invalid "seen" multi-index - seen_mi = [-1]*dim - # Local sum keep the sum of the terms that differ by each dimension - local_sum = [0]*dim - # Local multiplier keep the scalar that the sum has to be multiplied by - # when adding to the sum of the preceding dimension. - local_multiplier = [0]*dim - - # For the multi-indices in 3D, local_sum looks like this: - # - # Multi-index | coef | local_sum | local_mult - # (0, 0, 0) | c0 | 0, 0, c0 | 0, 1, 1 - # (0, 0, 1) | c1 | 0, 0, c0+c1*dz | 0, 1, 1 - # (0, 0, 2) | c2 | 0, 0, c0+c1*dz+c2*dz^2 | 0, 1, 1 - # (0, 1, 0) | c3 | 0, c0+c1*dz+c2*dz^2, c3 | 0, 1, dy - # (0, 1, 1) | c4 | 0, c0+c1*dz+c2*dz^2, c3+c4*dz | 0, 1, dy - # (0, 1, 2) | c5 | 0, c0+c1*dz+c2*dz^2, c3+c4*dz+c5*dz^2 | 0, 1, dy - # (0, 2, 0) | c6 | 0, c0+c1*dz+c2*dz^2, c6 | 0, 1, dy^2 - # | | +dy*(c3+c4*dz+c5*dz^2) | - # (0, 2, 1) | c7 | 0, c0+c1*dz+c2*dz^2, c6+c7*dz | 0, 1, dy^2 - # | | +dy*(c3+c4*dz+c5*dz^2) | - # (0, 2, 2) | c8 | 0, c0+c1*dz+c2*dz^2, c6+c7*dz+x8*dz^2 | 0, 1, dy^2 - # | | +dy*(c3+c4*dz+c5*dz^2) | - # (1, 0, 0) | c8 | c0+c1*dz+c2*dz^2, 0, 0 | 0, dx, 1 - # | | +dy*(c3+c4*dz+c5*dz^2) | - # | | +dy^2*(c6+c7*dz+c8*dz^2) | - - for mi in sorted_target_mis: - - # {{{ handle the case where a not-last dimension "clicked over" - - # (where d will be that not-last dimension) - - # Iterate in reverse order of dimensions to properly handle a - # "double click-over". - - for d in reversed(range(dim-1)): - if seen_mi[d] != mi[d]: - # If the dimension d of mi changed from the previous value - # then the sum for dimension d+1 is complete, add it to - # dimension d after multiplying and restart. - - local_sum[d] += local_sum[d+1]*local_multiplier[d+1] - local_sum[d+1] = 0 - local_multiplier[d+1] = bvec[d]**mi[d] / factorial(mi[d]) - - # }}} - - local_sum[dim-1] += evaluated_coeffs[mi_to_index[mi]] * \ - bvec[dim-1]**mi[dim-1] / factorial(mi[dim-1]) - seen_mi = mi - - for d in reversed(range(dim-1)): - local_sum[d] += local_sum[d+1]*local_multiplier[d+1] - - # {{{ simpler, functionally equivalent code - - if 0: - from sumpy.tools import mi_power, mi_factorial - - return sum( - coeff - * mi_power(bvec, mi, evaluate=False) - / mi_factorial(mi) - for coeff, mi in zip( - evaluated_coeffs, self.get_full_coefficient_identifiers())) - - # }}} + from sumpy.tools import mi_power, mi_factorial - return local_sum[0] + return sum( + coeff + * mi_power(bvec, mi, evaluate=False) + / mi_factorial(mi) + for coeff, mi in zip( + evaluated_coeffs, self.get_full_coefficient_identifiers())) def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, dvec, tgt_rscale): @@ -295,9 +224,68 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): lexp_mi_terms.append( coeff * kernel_deriv * tgt_rscale**sum(lexp_mi)) result.append(sym.Add(*lexp_mi_terms)) + logger.info("building translation operator: done") + return result - else: - from sumpy.tools import MiDerivativeTaker + rscale_ratio = sym.UnevaluatedExpr(tgt_rscale/src_rscale) + + from sumpy.tools import MiDerivativeTaker + from math import factorial + src_wrangler = src_expansion.expansion_terms_wrangler + src_coeffs = ( + src_wrangler.get_full_kernel_derivatives_from_stored( + src_coeff_exprs, src_rscale)) + src_mis = \ + src_expansion.expansion_terms_wrangler.get_full_coefficient_identifiers() + + src_mi_to_index = dict((mi, i) for i, mi in enumerate(src_mis)) + + tgt_mis_all = \ + self.expansion_terms_wrangler.get_coefficient_identifiers() + tgt_mi_to_index = dict((mi, i) for i, mi in enumerate(tgt_mis_all)) + + tgt_split = self.expansion_terms_wrangler._get_coeff_identifier_split() + + p = max(sum(mi) for mi in src_mis) + result = [0] * len(tgt_mis_all) + + # O(1) iterations + for const_dim in set(d for d, _ in tgt_split): + # Use the const_dim as the first dimension to vary so that the below + # algorithm is O(p^{d+1}) for full and O(p^{d}) for compressed + dims = [const_dim] + list(range(const_dim)) + \ + list(range(const_dim+1, self.dim)) + # Start with source coefficients + Y = src_coeffs # noqa: N806 + # O(1) iterations + for d in dims: + C = Y # noqa: N806 + Y = [0] * len(src_mis) # noqa: N806 + # Only O(p^{d-1}) operations are used in compressed + # All of them are used in full. + for i, s in enumerate(src_mis): + # O(p) iterations + for q in range(p+1-sum(s)): + src_mi = list(s) + src_mi[d] += q + src_mi = tuple(src_mi) + if src_mi in src_mi_to_index: + Y[i] += (dvec[d]/src_rscale) ** q * \ + C[src_mi_to_index[src_mi]] / factorial(q) + + # This is O(p) in full and O(1) in compressed + for d, tgt_mis in tgt_split: + if d != const_dim: + continue + # O(p^{d-1}) iterations + for mi in tgt_mis: + if mi not in src_mi_to_index: + continue + result[tgt_mi_to_index[mi]] = Y[src_mi_to_index[mi]] \ + * rscale_ratio ** sum(mi) + + # {{{ simpler, functionally equivalent code + if 0: # 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. Scaling dvec, then differentiating, @@ -314,7 +302,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): result = [ (taker.diff(mi).xreplace(replace_dict) * rscale_ratio**sum(mi)) for mi in self.get_coefficient_identifiers()] - + # }}} logger.info("building translation operator: done") return result @@ -369,7 +357,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): def get_coefficient_identifiers(self): return list(range(-self.order, self.order+1)) - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + def coefficients_from_source(self, avec, bvec, rscale): if not self.use_rscale: rscale = 1 @@ -387,7 +375,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): * sym.exp(sym.I * l * source_angle_rel_center), avec) for l in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): if not self.use_rscale: rscale = 1 diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 1a53aa40f15f78fb3d94c1b2f8bb9932d6daba90..b7e2769ba21b470cb163515da6e70c55eb360d2b 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -56,7 +56,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): Coefficients represent the terms in front of the kernel derivatives. """ - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + def coefficients_from_source(self, avec, bvec, rscale): from sumpy.kernel import DirectionalSourceDerivative kernel = self.kernel @@ -109,7 +109,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): for mi in self.get_full_coefficient_identifiers()] return ( self.expansion_terms_wrangler.get_stored_mpole_coefficients_from_full( - result, rscale, sac=sac)) + result, rscale)) def get_scaled_multipole(self, expr, bvec, rscale, nderivatives, nderivatives_for_scaling=None): @@ -127,7 +127,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): else: return (rscale**nderivatives_for_scaling * expr) - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): if not self.use_rscale: rscale = 1 @@ -164,18 +164,18 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): from sumpy.tools import mi_factorial src_mi_to_index = dict((mi, i) for i, mi in enumerate( - src_expansion.get_coefficient_identifiers())) + src_expansion.get_full_coefficient_identifiers())) tgt_mi_to_index = dict((mi, i) for i, mi in enumerate( self.get_full_coefficient_identifiers())) - src_coeff_exprs = list(src_coeff_exprs) + src_coeff_exprs_full = \ + [0]*len(src_expansion.get_full_coefficient_identifiers()) + for i, mi in enumerate(src_expansion.get_coefficient_identifiers()): - src_coeff_exprs[i] *= mi_factorial(mi) * \ + src_coeff_exprs_full[src_mi_to_index[mi]] = src_coeff_exprs[i] * \ sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(mi) - result = [0] * len(self.get_full_coefficient_identifiers()) - # This algorithm uses the observation that M2M coefficients # have the following form in 2D # @@ -183,43 +183,99 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): # d_x^i d_y^j \binom{m}{i} \binom{n}{j}$ # and can be rewritten as follows. # - # Let $T_{m, n} = \sum_{i\le m} A_{i, n} d_x^i \binom{m}{i}$. + # Let $Y_{m, n} = \sum_{i\le m} A_{i, n} d_x^i \binom{m}{i}$. # - # Then, $B_{m, n} = \sum_{j\le n} T_{m, j} d_y^j \binom{n}{j}$. + # Then, $B_{m, n} = \sum_{j\le n} Y_{m, j} d_y^j \binom{n}{j}$. # - # $T_{m, n}$ are $p^2$ number of temporary variables that are + # $Y_{m, n}$ are $p^2$ number of temporary variables that are # reused for different M2M coefficients and costs $p$ per variable. - # Total cost for calculating $T_{m, n}$ is $p^3$ and similar + # Total cost for calculating $Y_{m, n}$ is $p^3$ and similar # for $B_{m, n}$ # In other words, we're better off computing the translation - # one dimension at a time. This is realized here by using - # the output from one dimension as the input to the next. - # per_dim_coeffs_to_translate serves as the array of input - # coefficients for each dimension's translation. - - dim_coeffs_to_translate = src_coeff_exprs + # one dimension at a time. If the coefficients in the source + # expansion have the form (0, m) and (m, 0), then we calculate + # the output from (0, m) with the second dimension as the fastest + # varying dimension and then calculate the output from (m, 0) + # with the first dimension as the fastest varying dimension. + # However, the contribution from (0, 0) is counted twice and we + # need to subtract that. + + src_split = \ + src_expansion.expansion_terms_wrangler._get_coeff_identifier_split() + result = [0] * len(self.get_full_coefficient_identifiers()) - mi_to_index = src_mi_to_index - for d in range(self.dim): + non_zero_coeffs_per_dim = [[] for d in range(self.dim)] + for d, src_mi in src_split: + non_zero_coeffs_per_dim[d] += src_mi + + for const_dim in set(d for d, _ in src_split): + dim_coeffs_to_translate = \ + [0] * len(src_expansion.get_full_coefficient_identifiers()) + for mi in non_zero_coeffs_per_dim[const_dim]: + idx = src_mi_to_index[mi] + dim_coeffs_to_translate[idx] = src_coeff_exprs_full[idx] + + # Use the const_dim as the last dimension to vary + dims = list(range(const_dim)) + \ + list(range(const_dim+1, self.dim)) + [const_dim] + for d in dims: + temp = [0] * len(src_expansion.get_full_coefficient_identifiers()) + for i, tgt_mi in enumerate( + src_expansion.get_full_coefficient_identifiers()): + src_mis_per_dim = [] + for mi_i in range(tgt_mi[d]+1): + new_mi = list(tgt_mi) + new_mi[d] = mi_i + src_mis_per_dim.append(tuple(new_mi)) + + for src_mi in src_mis_per_dim: + try: + src_index = src_mi_to_index[src_mi] + except KeyError: + # Omitted coefficients: not life-threatening + continue + + contrib = dim_coeffs_to_translate[src_index] + for idim in range(self.dim): + n = tgt_mi[idim] + k = src_mi[idim] + assert n >= k + contrib /= mi_factorial((n-k,)) + contrib *= \ + sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k) + + temp[i] += contrib + + dim_coeffs_to_translate = temp[:] + + for i, mi in enumerate(src_expansion.get_full_coefficient_identifiers()): + result[tgt_mi_to_index[mi]] += dim_coeffs_to_translate[i] + + # {{{ simpler, functionally equivalent code + if 0: + src_mi_to_index = dict((mi, i) for i, mi in enumerate( + src_expansion.get_coefficient_identifiers())) result = [0] * len(self.get_full_coefficient_identifiers()) + + for i, mi in enumerate(src_expansion.get_coefficient_identifiers()): + src_coeff_exprs[i] *= mi_factorial(mi) + + from pytools import generate_nonnegative_integer_tuples_below as gnitb + for i, tgt_mi in enumerate( self.get_full_coefficient_identifiers()): - src_mis_per_dim = [] - for mi_i in range(tgt_mi[d]+1): - new_mi = list(tgt_mi) - new_mi[d] = mi_i - src_mis_per_dim.append(tuple(new_mi)) + tgt_mi_plus_one = tuple(mi_i + 1 for mi_i in tgt_mi) - for src_mi in src_mis_per_dim: + for src_mi in gnitb(tgt_mi_plus_one): try: - src_index = mi_to_index[src_mi] + src_index = src_mi_to_index[src_mi] except KeyError: # Omitted coefficients: not life-threatening continue - contrib = dim_coeffs_to_translate[src_index] + contrib = src_coeff_exprs[src_index] for idim in range(self.dim): n = tgt_mi[idim] @@ -229,14 +285,11 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): contrib *= (binomial(n, k) * sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k)) - result[i] += contrib - - # Defer division by target factorial until the very end - if d == self.dim-1: - result[i] /= mi_factorial(tgt_mi) + result[i] += (contrib + * sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(src_mi)) - dim_coeffs_to_translate = result[:] - mi_to_index = tgt_mi_to_index + result[i] /= mi_factorial(tgt_mi) + # }}} logger.info("building translation operator: done") return ( @@ -294,7 +347,7 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): def get_coefficient_identifiers(self): return list(range(-self.order, self.order+1)) - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + def coefficients_from_source(self, avec, bvec, rscale): if not self.use_rscale: rscale = 1 @@ -314,7 +367,7 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): avec) for l in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): if not self.use_rscale: rscale = 1 diff --git a/sumpy/p2e.py b/sumpy/p2e.py index f5472fc42e306c21eeec5091b26d548cdb66c71a..2aa648040a7b51645345160ac8288d877ea17c11 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -88,7 +88,7 @@ class P2EBase(KernelCacheWrapper): coeff_names = [ sac.assign_unique("coeff%d" % i, coeff_i) for i, coeff_i in enumerate( - self.expansion.coefficients_from_source(avec, None, rscale, sac))] + self.expansion.coefficients_from_source(avec, None, rscale))] sac.run_global_cse()