diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index b67fae02b5069f2a5a92c6b63bea405000a379d6..c1281f9691613ffcb5c3e8ab4a50e7ce7e3bd3da 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -23,6 +23,8 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +import numpy as np +import sympy as sp from pytools import memoize_method __doc__ = """ @@ -33,6 +35,7 @@ __doc__ = """ # {{{ base class class ExpansionBase(object): + def __init__(self, kernel, order): # Don't be tempted to remove target derivatives here. # Line Taylor QBX can't do without them, because it can't @@ -106,30 +109,161 @@ class ExpansionBase(object): def __ne__(self, other): return not self.__eq__(other) - # }}} # {{{ volume taylor class VolumeTaylorExpansionBase(object): - @memoize_method - def _storage_loc_dict(self): - return dict( - (idx, i) - for i, idx in enumerate(self.get_coefficient_identifiers())) - def get_storage_index(self, k): - return self._storage_loc_dict()[k] + def get_coefficient_identifiers(self): + """ + Returns the identifiers of the coefficients that actually get stored. + """ + raise NotImplementedError @memoize_method - def get_coefficient_identifiers(self): + def get_full_coefficient_identifiers(self): + """ + Returns identifiers for every coefficient in the complete expansion. + """ from pytools import ( generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam) return sorted(gnitstam(self.order, self.kernel.dim), key=sum) + def stored_to_full(self, coeff_idx, stored_coeffs): + raise NotImplementedError + + def full_to_stored(self, coeff_idx, full_coeffs): + raise NotImplementedError + + +class VolumeTaylorExpansion(VolumeTaylorExpansionBase): + + get_coefficient_identifiers = ( + VolumeTaylorExpansionBase.get_full_coefficient_identifiers) + + def stored_to_full(self, stored_coeffs): + return stored_coeffs + + full_to_stored = stored_to_full + + +class LinearRecurrenceBasedVolumeTaylorExpansion(VolumeTaylorExpansionBase): + + def __init__(self): + self.precompute_coeff_matrix() + + def get_coefficient_identifiers(self): + return self.stored_identifiers + + def stored_to_full(self, stored_coeffs): + return self.coeff_matrix.dot(stored_coeffs) + + def full_to_stored(self, full_coeffs): + return self.coeff_matrix.T.dot(full_coeffs) + + def precompute_coeff_matrix(self): + stored_identifiers = [] + identifiers_so_far = [] + + ncoeffs = len(self.get_full_coefficient_identifiers()) + coeff_matrix = [] + + # Build up the matrix by row. + for identifier in self.get_full_coefficient_identifiers(): + expr = self.try_get_recurrence_for_derivative( + identifier, identifiers_so_far, ncoeffs) + + if expr is None: + # Identifier should be stored + row = np.zeros(ncoeffs, dtype=object) + row[len(stored_identifiers)] = 1 + stored_identifiers.append(identifier) + else: + # Rewrite in terms of the stored identifiers + ncoeffs_so_far = len(coeff_matrix) + row = np.dot(np.transpose(coeff_matrix), expr[:ncoeffs_so_far]) + + coeff_matrix.append(row) + identifiers_so_far.append(identifier) + + self.stored_identifiers = stored_identifiers + ncols = len(stored_identifiers) + self.coeff_matrix = np.vstack(coeff_matrix)[:, :ncols] + + +class LaplaceConformingVolumeTaylorExpansion( + LinearRecurrenceBasedVolumeTaylorExpansion): + + def try_get_recurrence_for_derivative(self, deriv, in_terms_of, ncoeffs): + deriv = np.array(deriv) + + for dim in np.nonzero(2 <= deriv)[0]: + # Check if we can reduce this dimension in terms of the other + # dimensions. + + reduced_deriv = deriv.copy() + reduced_deriv[dim] -= 2 + + needed_derivs = [] + for other_dim in range(self.kernel.dim): + if other_dim == dim: + continue + needed_deriv = reduced_deriv.copy() + needed_deriv[other_dim] += 2 + + needed_derivs.append(tuple(needed_deriv)) + + expr = np.zeros(ncoeffs, dtype=object) + try: + for needed_deriv in needed_derivs: + deriv_idx = in_terms_of.index(needed_deriv) + expr[deriv_idx] = -1 + except ValueError: + continue + + return expr + + +class HelmholtzConformingVolumeTaylorExpansion( + LinearRecurrenceBasedVolumeTaylorExpansion): + + def try_get_recurrence_for_derivative(self, deriv, in_terms_of, ncoeffs): + deriv = np.array(deriv) + + for dim in np.nonzero(2 <= deriv)[0]: + # Check if we can reduce this dimension in terms of the other + # dimensions. + + reduced_deriv = deriv.copy() + reduced_deriv[dim] -= 2 + + needed_derivs = [] + for other_dim in range(self.kernel.dim): + if other_dim == dim: + continue + needed_deriv = reduced_deriv.copy() + needed_deriv[other_dim] += 2 + + needed_derivs.append((-1, tuple(needed_deriv))) + + k = sp.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + + needed_derivs.append((-k*k, tuple(reduced_deriv))) + + expr = np.zeros(ncoeffs, dtype=object) + try: + for coeff, needed_deriv in needed_derivs: + deriv_idx = in_terms_of.index(needed_deriv) + expr[deriv_idx] = coeff + except ValueError: + continue + + return expr + # }}} diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 92cae781c81832124cc925cdeaaadfa1ab9bfbb8..2bb87b46f4f7de7d224b6cb8584951965d65e09f 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -25,7 +25,9 @@ THE SOFTWARE. import sympy as sp -from sumpy.expansion import ExpansionBase, VolumeTaylorExpansionBase +from sumpy.expansion import ( + ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, + HelmholtzConformingVolumeTaylorExpansion) class LocalExpansionBase(ExpansionBase): @@ -46,6 +48,7 @@ __doc__ = """ # {{{ line taylor class LineTaylorLocalExpansion(LocalExpansionBase): + def get_storage_index(self, k): return k @@ -81,7 +84,11 @@ class LineTaylorLocalExpansion(LocalExpansionBase): # {{{ volume taylor -class VolumeTaylorLocalExpansion(LocalExpansionBase, VolumeTaylorExpansionBase): +class VolumeTaylorLocalExpansionBase(LocalExpansionBase): + """ + Coefficients represent derivative values of the kernel. + """ + def coefficients_from_source(self, avec, bvec): from sumpy.tools import mi_derivative ppkernel = self.kernel.postprocess_at_source( @@ -91,11 +98,14 @@ class VolumeTaylorLocalExpansion(LocalExpansionBase, VolumeTaylorExpansionBase): def evaluate(self, coeffs, bvec): from sumpy.tools import mi_power, mi_factorial - return sum( + evaluated_coeffs = self.stored_to_full(coeffs) + result = sum( coeff * self.kernel.postprocess_at_target(mi_power(bvec, mi), bvec) / mi_factorial(mi) - for coeff, mi in zip(coeffs, self.get_coefficient_identifiers())) + for coeff, mi in zip( + evaluated_coeffs, self.get_full_coefficient_identifiers())) + return result def translate_from(self, src_expansion, src_coeff_exprs, dvec): logger.info("building translation operator: %s(%d) -> %s(%d): start" @@ -112,6 +122,33 @@ class VolumeTaylorLocalExpansion(LocalExpansionBase, VolumeTaylorExpansionBase): logger.info("building translation operator: done") return result + +class VolumeTaylorLocalExpansion( + VolumeTaylorLocalExpansionBase, + VolumeTaylorExpansion): + + def __init__(self, kernel, order): + VolumeTaylorLocalExpansionBase.__init__(self, kernel, order) + VolumeTaylorExpansion.__init__(self) + + +class LaplaceConformingVolumeTaylorLocalExpansion( + VolumeTaylorLocalExpansionBase, + LaplaceConformingVolumeTaylorExpansion): + + def __init__(self, kernel, order): + VolumeTaylorLocalExpansionBase.__init__(self, kernel, order) + LaplaceConformingVolumeTaylorExpansion.__init__(self) + + +class HelmholtzConformingVolumeTaylorLocalExpansion( + VolumeTaylorLocalExpansionBase, + HelmholtzConformingVolumeTaylorExpansion): + + def __init__(self, kernel, order): + VolumeTaylorLocalExpansionBase.__init__(self, kernel, order) + HelmholtzConformingVolumeTaylorExpansion.__init__(self) + # }}} diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 4d6594017d664c91fd238776ace4bdf6d38ca091..f579bd304d404a862196aa46eb136c5d49194ebc 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -25,7 +25,9 @@ THE SOFTWARE. from six.moves import range, zip import sympy as sp # noqa -from sumpy.expansion import ExpansionBase, VolumeTaylorExpansionBase +from sumpy.expansion import ( + ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, + HelmholtzConformingVolumeTaylorExpansion) import logging logger = logging.getLogger(__name__) @@ -45,13 +47,16 @@ class MultipoleExpansionBase(ExpansionBase): # {{{ volume taylor -class VolumeTaylorMultipoleExpansion( - MultipoleExpansionBase, VolumeTaylorExpansionBase): +class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): + """ + Coefficients represent the terms in front of the kernel derivatives. + """ + def coefficients_from_source(self, avec, bvec): from sumpy.kernel import DirectionalSourceDerivative kernel = self.kernel - from sumpy.tools import mi_power + from sumpy.tools import mi_power, mi_factorial if isinstance(kernel, DirectionalSourceDerivative): if kernel.get_base_kernel() is not kernel.inner_kernel: @@ -61,7 +66,7 @@ class VolumeTaylorMultipoleExpansion( from sumpy.symbolic import make_sympy_vector dir_vec = make_sympy_vector(kernel.dir_vec_name, kernel.dim) - coeff_identifiers = self.get_coefficient_identifiers() + coeff_identifiers = self.get_full_coefficient_identifiers() result = [0] * len(coeff_identifiers) for idim in range(kernel.dim): @@ -75,22 +80,23 @@ class VolumeTaylorMultipoleExpansion( result[i] += ( - mi_power(avec, derivative_mi) * mi[idim] * dir_vec[idim]) - - return result + for i, mi in enumerate(coeff_identifiers): + result[i] /= mi_factorial(mi) else: - return [ - mi_power(avec, mi) - for mi in self.get_coefficient_identifiers()] + result = [ + mi_power(avec, mi) / mi_factorial(mi) + for mi in self.get_full_coefficient_identifiers()] + return self.full_to_stored(result) def evaluate(self, coeffs, bvec): ppkernel = self.kernel.postprocess_at_target( self.kernel.get_expression(bvec), bvec) - from sumpy.tools import mi_derivative, mi_factorial - return sum( - coeff / mi_factorial(mi) - * mi_derivative(ppkernel, bvec, mi) + from sumpy.tools import mi_derivative + result = sum( + coeff * mi_derivative(ppkernel, bvec, mi) for coeff, mi in zip(coeffs, self.get_coefficient_identifiers())) + return result def translate_from(self, src_expansion, src_coeff_exprs, dvec): if not isinstance(src_expansion, type(self)): @@ -104,14 +110,19 @@ class VolumeTaylorMultipoleExpansion( type(self).__name__, self.order)) + from sumpy.tools import mi_factorial + src_mi_to_index = dict((mi, i) for i, mi in enumerate( src_expansion.get_coefficient_identifiers())) - result = [0] * len(self.get_coefficient_identifiers()) + for i, mi in enumerate(src_expansion.get_coefficient_identifiers()): + src_coeff_exprs[i] *= mi_factorial(mi) + + result = [0] * len(self.get_full_coefficient_identifiers()) from pytools import generate_nonnegative_integer_tuples_below as gnitb for i, tgt_mi in enumerate( - self.get_coefficient_identifiers()): + self.get_full_coefficient_identifiers()): tgt_mi_plus_one = tuple(mi_i + 1 for mi_i in tgt_mi) @@ -119,7 +130,7 @@ class VolumeTaylorMultipoleExpansion( try: src_index = src_mi_to_index[src_mi] except KeyError: - # tgt higher-order than source--dumb, but not life-threatening + # Omitted coefficients: not life-threatening continue contrib = src_coeff_exprs[src_index] @@ -133,8 +144,37 @@ class VolumeTaylorMultipoleExpansion( result[i] += contrib + result[i] /= mi_factorial(tgt_mi) + logger.info("building translation operator: done") - return result + return self.full_to_stored(result) + + +class VolumeTaylorMultipoleExpansion( + VolumeTaylorMultipoleExpansionBase, + VolumeTaylorExpansion): + + def __init__(self, kernel, order): + VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order) + VolumeTaylorExpansion.__init__(self) + + +class LaplaceConformingVolumeTaylorMultipoleExpansion( + VolumeTaylorMultipoleExpansionBase, + LaplaceConformingVolumeTaylorExpansion): + + def __init__(self, kernel, order): + VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order) + LaplaceConformingVolumeTaylorExpansion.__init__(self) + + +class HelmholtzConformingVolumeTaylorMultipoleExpansion( + VolumeTaylorMultipoleExpansionBase, + HelmholtzConformingVolumeTaylorExpansion): + + def __init__(self, kernel, order): + VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order) + HelmholtzConformingVolumeTaylorExpansion.__init__(self) # }}} diff --git a/test/test_fmm.py b/test/test_fmm.py index 6e292fbae08abcdbc874790e9c522f386d0c109b..2e2bfea4803b6a5f5e17a6d05379620adff09855 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -32,9 +32,13 @@ from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) from sumpy.kernel import LaplaceKernel, HelmholtzKernel from sumpy.expansion.multipole import ( - VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion) + VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion, + LaplaceConformingVolumeTaylorMultipoleExpansion, + HelmholtzConformingVolumeTaylorMultipoleExpansion) from sumpy.expansion.local import ( - VolumeTaylorLocalExpansion, H2DLocalExpansion) + VolumeTaylorLocalExpansion, H2DLocalExpansion, + LaplaceConformingVolumeTaylorLocalExpansion, + HelmholtzConformingVolumeTaylorLocalExpansion) from sumpy.expansion.level_to_order import ( h2d_level_to_order_lookup, @@ -84,10 +88,18 @@ def test_level_to_order_lookup(ctx_getter, lookup_func, extra_args): @pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class", [ (LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), + (LaplaceKernel(2), LaplaceConformingVolumeTaylorLocalExpansion, + LaplaceConformingVolumeTaylorMultipoleExpansion), (LaplaceKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), + (LaplaceKernel(3), LaplaceConformingVolumeTaylorLocalExpansion, + LaplaceConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), + (HelmholtzKernel(2), HelmholtzConformingVolumeTaylorLocalExpansion, + HelmholtzConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), H2DLocalExpansion, H2DMultipoleExpansion), - (HelmholtzKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion) + (HelmholtzKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), + (HelmholtzKernel(3), HelmholtzConformingVolumeTaylorLocalExpansion, + HelmholtzConformingVolumeTaylorMultipoleExpansion), ]) def test_sumpy_fmm(ctx_getter, knl, local_expn_class, mpole_expn_class): logging.basicConfig(level=logging.INFO) diff --git a/test/test_kernels.py b/test/test_kernels.py index 465ce178443dd321c4860bc5ec153aa492e2e9db..b93008925936589be276f98fb847e20652489075 100644 --- a/test/test_kernels.py +++ b/test/test_kernels.py @@ -35,8 +35,14 @@ from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) from sumpy.expansion.multipole import ( - VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion) -from sumpy.expansion.local import VolumeTaylorLocalExpansion, H2DLocalExpansion + VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion, + VolumeTaylorMultipoleExpansionBase, + LaplaceConformingVolumeTaylorMultipoleExpansion, + HelmholtzConformingVolumeTaylorMultipoleExpansion) +from sumpy.expansion.local import ( + VolumeTaylorLocalExpansion, H2DLocalExpansion, + LaplaceConformingVolumeTaylorLocalExpansion, + HelmholtzConformingVolumeTaylorLocalExpansion) from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative, DirectionalSourceDerivative) from pytools.convergence import PConvergenceVerifier @@ -95,14 +101,22 @@ def test_p2p(ctx_getter): @pytest.mark.parametrize(("base_knl", "expn_class"), [ (LaplaceKernel(2), VolumeTaylorLocalExpansion), (LaplaceKernel(2), VolumeTaylorMultipoleExpansion), + (LaplaceKernel(2), LaplaceConformingVolumeTaylorLocalExpansion), + (LaplaceKernel(2), LaplaceConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), VolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), VolumeTaylorLocalExpansion), + (HelmholtzKernel(2), HelmholtzConformingVolumeTaylorLocalExpansion), + (HelmholtzKernel(2), HelmholtzConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), H2DLocalExpansion), (HelmholtzKernel(2), H2DMultipoleExpansion), (HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorMultipoleExpansion), (HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorLocalExpansion), + (HelmholtzKernel(2, allow_evanescent=True), + HelmholtzConformingVolumeTaylorLocalExpansion), + (HelmholtzKernel(2, allow_evanescent=True), + HelmholtzConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2, allow_evanescent=True), H2DLocalExpansion), (HelmholtzKernel(2, allow_evanescent=True), H2DMultipoleExpansion), ]) @@ -292,11 +306,12 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative): slack += 1 grad_slack += 2 - if isinstance(base_knl, HelmholtzKernel) and base_knl.allow_evanescent: - slack += 0.5 - grad_slack += 0.5 + if isinstance(base_knl, HelmholtzKernel): + if base_knl.allow_evanescent: + slack += 0.5 + grad_slack += 0.5 - if expn_class is VolumeTaylorMultipoleExpansion: + if issubclass(expn_class, VolumeTaylorMultipoleExpansionBase): slack += 0.3 grad_slack += 0.3 @@ -306,7 +321,11 @@ def test_p2e2p(ctx_getter, base_knl, expn_class, order, with_source_derivative): @pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class", [ (LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), + (LaplaceKernel(2), LaplaceConformingVolumeTaylorLocalExpansion, + LaplaceConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), + (HelmholtzKernel(2), HelmholtzConformingVolumeTaylorLocalExpansion, + HelmholtzConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), H2DLocalExpansion, H2DMultipoleExpansion) ]) def test_translations(ctx_getter, knl, local_expn_class, mpole_expn_class):