diff --git a/benchmarks/bench_translations.py b/benchmarks/bench_translations.py index ed97ed5456d00fcdc99daeb8e14841db8155b4dd..1c66ba747e2c61998e3106ddd9ff5becea39d7ca 100644 --- a/benchmarks/bench_translations.py +++ b/benchmarks/bench_translations.py @@ -5,12 +5,10 @@ from pyopencl.tools import ( # noqa from sumpy.expansion.multipole import ( VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion, - LaplaceConformingVolumeTaylorMultipoleExpansion, - HelmholtzConformingVolumeTaylorMultipoleExpansion) + LinearPDEConformingVolumeTaylorMultipoleExpansion) from sumpy.expansion.local import ( VolumeTaylorLocalExpansion, H2DLocalExpansion, - LaplaceConformingVolumeTaylorLocalExpansion, - HelmholtzConformingVolumeTaylorLocalExpansion) + LinearPDEConformingVolumeTaylorLocalExpansion) from sumpy.kernel import LaplaceKernel, HelmholtzKernel @@ -91,8 +89,8 @@ class LaplaceVolumeTaylorTranslation(TranslationBenchmarkSuite): class LaplaceConformingVolumeTaylorTranslation(TranslationBenchmarkSuite): knl = LaplaceKernel - local_expn_class = LaplaceConformingVolumeTaylorLocalExpansion - mpole_expn_class = LaplaceConformingVolumeTaylorMultipoleExpansion + local_expn_class = LinearPDEConformingVolumeTaylorLocalExpansion + mpole_expn_class = LinearPDEConformingVolumeTaylorMultipoleExpansion class HelmholtzVolumeTaylorTranslation(TranslationBenchmarkSuite): @@ -107,8 +105,8 @@ class HelmholtzVolumeTaylorTranslation(TranslationBenchmarkSuite): class HelmholtzConformingVolumeTaylorTranslation(TranslationBenchmarkSuite): knl = HelmholtzKernel - local_expn_class = HelmholtzConformingVolumeTaylorLocalExpansion - mpole_expn_class = HelmholtzConformingVolumeTaylorMultipoleExpansion + local_expn_class = LinearPDEConformingVolumeTaylorLocalExpansion + mpole_expn_class = LinearPDEConformingVolumeTaylorMultipoleExpansion class Helmholtz2DTranslation(TranslationBenchmarkSuite): diff --git a/examples/sym-exp-complexity.py b/examples/sym-exp-complexity.py index 49e20ead986b880b86529ee8d436843d10ef3cbb..ae91a932cc4c646a6d6dcadce7ce72ed1acb40fa 100644 --- a/examples/sym-exp-complexity.py +++ b/examples/sym-exp-complexity.py @@ -3,12 +3,10 @@ import pyopencl as cl import loopy as lp from sumpy.kernel import LaplaceKernel, HelmholtzKernel from sumpy.expansion.local import ( - LaplaceConformingVolumeTaylorLocalExpansion, - HelmholtzConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorLocalExpansion, ) from sumpy.expansion.multipole import ( - LaplaceConformingVolumeTaylorMultipoleExpansion, - HelmholtzConformingVolumeTaylorMultipoleExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, ) from sumpy.e2e import E2EFromCSR try: @@ -22,13 +20,13 @@ def find_flops(): if 0: knl = LaplaceKernel(2) - m_expn_cls = LaplaceConformingVolumeTaylorMultipoleExpansion - l_expn_cls = LaplaceConformingVolumeTaylorLocalExpansion + m_expn_cls = LinearPDEConformingVolumeTaylorMultipoleExpansion + l_expn_cls = LinearPDEConformingVolumeTaylorLocalExpansion flop_type = np.float64 else: knl = HelmholtzKernel(2) - m_expn_cls = HelmholtzConformingVolumeTaylorMultipoleExpansion - l_expn_cls = HelmholtzConformingVolumeTaylorLocalExpansion + m_expn_cls = LinearPDEConformingVolumeTaylorMultipoleExpansion + l_expn_cls = LinearPDEConformingVolumeTaylorLocalExpansion flop_type = np.complex128 orders = list(range(1, 11, 1)) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 79bc9cb8e79139f4b2470d4b88c69018992ca7f7..8bfb85e9670375a352a3ee48255b1ddc5a9e1f1b 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -24,7 +24,6 @@ import logging from pytools import memoize_method import sumpy.symbolic as sym from sumpy.tools import add_mi -from .diff_op import make_identity_diff_op, laplacian from typing import List, Tuple __doc__ = """ @@ -275,7 +274,7 @@ class ExpansionTermsWrangler: # to which it is orthogonal to and the constant `c` described above hyperplanes = [] if isinstance(self, LinearPDEBasedExpansionTermsWrangler): - pde_dict, = self.get_pde_as_diff_op().eqs + pde_dict, = self.knl.get_pde_as_diff_op().eqs if not all(ident.mi in mi_to_index for ident in pde_dict): # The order of the expansion is less than the order of the PDE. @@ -406,18 +405,18 @@ class CSEMatVecOperator: class LinearPDEBasedExpansionTermsWrangler(ExpansionTermsWrangler): """ .. automethod:: __init__ - .. automethod:: get_pde_as_diff_op """ - init_arg_names = ("order", "dim", "max_mi") + init_arg_names = ("order", "dim", "knl", "max_mi") - def __init__(self, order, dim, max_mi=None): + def __init__(self, order, dim, knl, max_mi=None): r""" :param order: order of the expansion :param dim: number of dimensions + :param knl: kernel for the PDE """ - super().__init__(order, dim, - max_mi) + super().__init__(order, dim, max_mi) + self.knl = knl def get_coefficient_identifiers(self): return self.stored_identifiers @@ -443,14 +442,6 @@ class LinearPDEBasedExpansionTermsWrangler(ExpansionTermsWrangler): stored_identifiers, _ = self.get_stored_ids_and_unscaled_projection_matrix() return stored_identifiers - def get_pde_as_diff_op(self): - r""" - Returns the PDE as a :class:`sumpy.expansion.diff_op.LinearPDESystemOperator` - object `L` where `L(u) = 0` is the PDE. - """ - - raise NotImplementedError - @memoize_method def get_stored_ids_and_unscaled_projection_matrix(self): from pytools import ProcessLogger @@ -460,7 +451,7 @@ class LinearPDEBasedExpansionTermsWrangler(ExpansionTermsWrangler): coeff_ident_enumerate_dict = {tuple(mi): i for (i, mi) in enumerate(mis)} - diff_op = self.get_pde_as_diff_op() + diff_op = self.knl.get_pde_as_diff_op() assert len(diff_op.eqs) == 1 pde_dict = {k.mi: v for k, v in diff_op.eqs[0].items()} for ident in pde_dict.keys(): @@ -582,46 +573,6 @@ class LinearPDEBasedExpansionTermsWrangler(ExpansionTermsWrangler): projection_with_rscale, shape) -class LaplaceExpansionTermsWrangler(LinearPDEBasedExpansionTermsWrangler): - - init_arg_names = ("order", "dim", "max_mi") - - def __init__(self, order, dim, max_mi=None): - super().__init__(order=order, dim=dim, - max_mi=max_mi) - - def get_pde_as_diff_op(self): - w = make_identity_diff_op(self.dim) - return laplacian(w) - - -class HelmholtzExpansionTermsWrangler(LinearPDEBasedExpansionTermsWrangler): - - init_arg_names = ("order", "dim", "helmholtz_k_name", "max_mi") - - def __init__(self, order, dim, helmholtz_k_name, max_mi=None): - self.helmholtz_k_name = helmholtz_k_name - super().__init__(order=order, dim=dim, - max_mi=max_mi) - - def get_pde_as_diff_op(self, **kwargs): - w = make_identity_diff_op(self.dim) - k = sym.Symbol(self.helmholtz_k_name) - return (laplacian(w) + k**2 * w) - - -class BiharmonicExpansionTermsWrangler(LinearPDEBasedExpansionTermsWrangler): - - init_arg_names = ("order", "dim", "max_mi") - - def __init__(self, order, dim, max_mi=None): - super().__init__(order=order, dim=dim, - max_mi=max_mi) - - def get_pde_as_diff_op(self, **kwargs): - w = make_identity_diff_op(self.dim) - return laplacian(laplacian(w)) - # }}} @@ -678,36 +629,47 @@ class VolumeTaylorExpansion(VolumeTaylorExpansionBase): self.expansion_terms_wrangler_key = (order, kernel.dim) -class LaplaceConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): +class LinearPDEConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): - expansion_terms_wrangler_class = LaplaceExpansionTermsWrangler + expansion_terms_wrangler_class = LinearPDEBasedExpansionTermsWrangler expansion_terms_wrangler_cache = {} # not user-facing, be strict about having to pass use_rscale def __init__(self, kernel, order, use_rscale): - self.expansion_terms_wrangler_key = (order, kernel.dim) + self.expansion_terms_wrangler_key = (order, kernel.dim, kernel) -class HelmholtzConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): +class LaplaceConformingVolumeTaylorExpansion( + LinearPDEConformingVolumeTaylorExpansion): - expansion_terms_wrangler_class = HelmholtzExpansionTermsWrangler - expansion_terms_wrangler_cache = {} + def __init__(self, *args, **kwargs): + from warnings import warn + warn("LaplaceConformingVolumeTaylorExpansion is deprecated. " + "Use LinearPDEConformingVolumeTaylorExpansion instead.", + DeprecationWarning, stacklevel=2) + super().__init__(*args, **kwargs) - # 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.expansion_terms_wrangler_key = (order, kernel.dim, helmholtz_k_name) +class HelmholtzConformingVolumeTaylorExpansion( + LinearPDEConformingVolumeTaylorExpansion): -class BiharmonicConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): + def __init__(self, *args, **kwargs): + from warnings import warn + warn("HelmholtzConformingVolumeTaylorExpansion is deprecated. " + "Use LinearPDEConformingVolumeTaylorExpansion instead.", + DeprecationWarning, stacklevel=2) + super().__init__(*args, **kwargs) - expansion_terms_wrangler_class = BiharmonicExpansionTermsWrangler - expansion_terms_wrangler_cache = {} - # not user-facing, be strict about having to pass use_rscale - def __init__(self, kernel, order, use_rscale): - self.expansion_terms_wrangler_key = (order, kernel.dim) +class BiharmonicConformingVolumeTaylorExpansion( + LinearPDEConformingVolumeTaylorExpansion): + def __init__(self, *args, **kwargs): + from warnings import warn + warn("BiharmonicConformingVolumeTaylorExpansion is deprecated. " + "Use LinearPDEConformingVolumeTaylorExpansion instead.", + DeprecationWarning, stacklevel=2) + super().__init__(*args, **kwargs) # }}} @@ -757,13 +719,10 @@ class DefaultExpansionFactory(ExpansionFactoryBase): def get_local_expansion_class(self, base_kernel): """Returns a subclass of :class:`ExpansionBase` suitable for *base_kernel*. """ - from sumpy.kernel import (HelmholtzKernel, LaplaceKernel, YukawaKernel, - BiharmonicKernel, StokesletKernel, StressletKernel) + from sumpy.kernel import (HelmholtzKernel, YukawaKernel) from sumpy.expansion.local import (H2DLocalExpansion, Y2DLocalExpansion, - HelmholtzConformingVolumeTaylorLocalExpansion, - LaplaceConformingVolumeTaylorLocalExpansion, - BiharmonicConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorLocalExpansion, VolumeTaylorLocalExpansion) if (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel) @@ -772,27 +731,20 @@ class DefaultExpansionFactory(ExpansionFactoryBase): elif (isinstance(base_kernel.get_base_kernel(), YukawaKernel) and base_kernel.dim == 2): return Y2DLocalExpansion - elif isinstance(base_kernel.get_base_kernel(), HelmholtzKernel): - return HelmholtzConformingVolumeTaylorLocalExpansion - elif isinstance(base_kernel.get_base_kernel(), LaplaceKernel): - return LaplaceConformingVolumeTaylorLocalExpansion - elif isinstance(base_kernel.get_base_kernel(), - (BiharmonicKernel, StokesletKernel, StressletKernel)): - return BiharmonicConformingVolumeTaylorLocalExpansion - else: + try: + base_kernel.get_base_kernel().get_pde_as_diff_op() + return LinearPDEConformingVolumeTaylorLocalExpansion + except NotImplementedError: return VolumeTaylorLocalExpansion def get_multipole_expansion_class(self, base_kernel): """Returns a subclass of :class:`ExpansionBase` suitable for *base_kernel*. """ - from sumpy.kernel import (HelmholtzKernel, LaplaceKernel, YukawaKernel, - BiharmonicKernel, StokesletKernel, StressletKernel) + from sumpy.kernel import (HelmholtzKernel, YukawaKernel) from sumpy.expansion.multipole import (H2DMultipoleExpansion, Y2DMultipoleExpansion, - LaplaceConformingVolumeTaylorMultipoleExpansion, - HelmholtzConformingVolumeTaylorMultipoleExpansion, - BiharmonicConformingVolumeTaylorMultipoleExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, VolumeTaylorMultipoleExpansion) if (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel) @@ -801,14 +753,10 @@ class DefaultExpansionFactory(ExpansionFactoryBase): elif (isinstance(base_kernel.get_base_kernel(), YukawaKernel) and base_kernel.dim == 2): return Y2DMultipoleExpansion - elif isinstance(base_kernel.get_base_kernel(), LaplaceKernel): - return LaplaceConformingVolumeTaylorMultipoleExpansion - elif isinstance(base_kernel.get_base_kernel(), HelmholtzKernel): - return HelmholtzConformingVolumeTaylorMultipoleExpansion - elif isinstance(base_kernel.get_base_kernel(), - (BiharmonicKernel, StokesletKernel, StressletKernel)): - return BiharmonicConformingVolumeTaylorMultipoleExpansion - else: + try: + base_kernel.get_base_kernel().get_pde_as_diff_op() + return LinearPDEConformingVolumeTaylorMultipoleExpansion + except NotImplementedError: return VolumeTaylorMultipoleExpansion # }}} diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 4d9aca2d9b1cdec15303dabfa2230812646c50d1..44994f6cac3a71345089e2cf4bff73ccb7184536 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -24,7 +24,8 @@ import sumpy.symbolic as sym from sumpy.tools import add_to_sac from sumpy.expansion import ( - ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, + ExpansionBase, VolumeTaylorExpansion, LinearPDEConformingVolumeTaylorExpansion, + LaplaceConformingVolumeTaylorExpansion, HelmholtzConformingVolumeTaylorExpansion, BiharmonicConformingVolumeTaylorExpansion) @@ -118,7 +119,7 @@ class LineTaylorLocalExpansion(LocalExpansionBase): .subs("tau", 0) for i in self.get_coefficient_identifiers()] - def evaluate(self, kernel, coeffs, bvec, rscale, sac=None): + def evaluate(self, tgt_kernel, coeffs, bvec, rscale, sac=None): # no point in heeding rscale here--just ignore it from pytools import factorial return sym.Add(*( @@ -493,34 +494,47 @@ class VolumeTaylorLocalExpansion( VolumeTaylorExpansion.__init__(self, kernel, order, use_rscale) -class LaplaceConformingVolumeTaylorLocalExpansion( - LaplaceConformingVolumeTaylorExpansion, +class LinearPDEConformingVolumeTaylorLocalExpansion( + LinearPDEConformingVolumeTaylorExpansion, VolumeTaylorLocalExpansionBase): def __init__(self, kernel, order, use_rscale=None): VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale) - LaplaceConformingVolumeTaylorExpansion.__init__( + LinearPDEConformingVolumeTaylorExpansion.__init__( self, kernel, order, use_rscale) +class LaplaceConformingVolumeTaylorLocalExpansion( + LaplaceConformingVolumeTaylorExpansion): + + def __init__(self, *args, **kwargs): + from warnings import warn + warn("LaplaceConformingVolumeTaylorLocalExpansion is deprecated. " + "Use LinearPDEConformingVolumeTaylorLocalExpansion instead.", + DeprecationWarning, stacklevel=2) + super().__init__(*args, **kwargs) + + class HelmholtzConformingVolumeTaylorLocalExpansion( - HelmholtzConformingVolumeTaylorExpansion, - VolumeTaylorLocalExpansionBase): + HelmholtzConformingVolumeTaylorExpansion): - def __init__(self, kernel, order, use_rscale=None): - VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale) - HelmholtzConformingVolumeTaylorExpansion.__init__( - self, kernel, order, use_rscale) + def __init__(self, *args, **kwargs): + from warnings import warn + warn("HelmholtzConformingVolumeTaylorLocalExpansion is deprecated. " + "Use LinearPDEConformingVolumeTaylorLocalExpansion instead.", + DeprecationWarning, stacklevel=2) + super().__init__(*args, **kwargs) class BiharmonicConformingVolumeTaylorLocalExpansion( - BiharmonicConformingVolumeTaylorExpansion, - VolumeTaylorLocalExpansionBase): - - def __init__(self, kernel, order, use_rscale=None): - VolumeTaylorLocalExpansionBase.__init__(self, kernel, order, use_rscale) - BiharmonicConformingVolumeTaylorExpansion.__init__( - self, kernel, order, use_rscale) + BiharmonicConformingVolumeTaylorExpansion): + + def __init__(self, *args, **kwargs): + from warnings import warn + warn("BiharmonicConformingVolumeTaylorLocalExpansion is deprecated. " + "Use LinearPDEConformingVolumeTaylorLocalExpansion instead.", + DeprecationWarning, stacklevel=2) + super().__init__(*args, **kwargs) # }}} diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index b9b1a07700b1019d18ea06b2613cd9acca4f251c..54e8b6e3c689a43497f5513fef9411b4ed252beb 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -23,7 +23,8 @@ THE SOFTWARE. import sumpy.symbolic as sym # noqa from sumpy.expansion import ( - ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, + ExpansionBase, VolumeTaylorExpansion, LinearPDEConformingVolumeTaylorExpansion, + LaplaceConformingVolumeTaylorExpansion, HelmholtzConformingVolumeTaylorExpansion, BiharmonicConformingVolumeTaylorExpansion) from pytools import factorial @@ -346,34 +347,47 @@ class VolumeTaylorMultipoleExpansion( VolumeTaylorExpansion.__init__(self, kernel, order, use_rscale) -class LaplaceConformingVolumeTaylorMultipoleExpansion( - LaplaceConformingVolumeTaylorExpansion, +class LinearPDEConformingVolumeTaylorMultipoleExpansion( + LinearPDEConformingVolumeTaylorExpansion, VolumeTaylorMultipoleExpansionBase): def __init__(self, kernel, order, use_rscale=None): VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale) - LaplaceConformingVolumeTaylorExpansion.__init__( + LinearPDEConformingVolumeTaylorExpansion.__init__( self, kernel, order, use_rscale) +class LaplaceConformingVolumeTaylorMultipoleExpansion( + LaplaceConformingVolumeTaylorExpansion): + + def __init__(self, *args, **kwargs): + from warnings import warn + warn("LaplaceConformingVolumeTaylorMultipoleExpansion is deprecated. " + "Use LinearPDEConformingVolumeTaylorMultipoleExpansion instead.", + DeprecationWarning, stacklevel=2) + super().__init__(*args, **kwargs) + + class HelmholtzConformingVolumeTaylorMultipoleExpansion( - HelmholtzConformingVolumeTaylorExpansion, - VolumeTaylorMultipoleExpansionBase): + HelmholtzConformingVolumeTaylorExpansion): - def __init__(self, kernel, order, use_rscale=None): - VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale) - HelmholtzConformingVolumeTaylorExpansion.__init__( - self, kernel, order, use_rscale) + def __init__(self, *args, **kwargs): + from warnings import warn + warn("HelmholtzConformingVolumeTaylorMultipoleExpansion is deprecated. " + "Use LinearPDEConformingVolumeTaylorMultipoleExpansion instead.", + DeprecationWarning, stacklevel=2) + super().__init__(*args, **kwargs) class BiharmonicConformingVolumeTaylorMultipoleExpansion( - BiharmonicConformingVolumeTaylorExpansion, - VolumeTaylorMultipoleExpansionBase): - - def __init__(self, kernel, order, use_rscale=None): - VolumeTaylorMultipoleExpansionBase.__init__(self, kernel, order, use_rscale) - BiharmonicConformingVolumeTaylorExpansion.__init__( - self, kernel, order, use_rscale) + BiharmonicConformingVolumeTaylorExpansion): + + def __init__(self, *args, **kwargs): + from warnings import warn + warn("BiharmonicConformingVolumeTaylorMultipoleExpansion is deprecated. " + "Use LinearPDEConformingVolumeTaylorMultipoleExpansion instead.", + DeprecationWarning, stacklevel=2) + super().__init__(*args, **kwargs) # }}} diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 494a6107248004883cd70c9e79dc58474c8d74ea..a74025e8c490ec0b94e8f69b34eefce2ac8294b9 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -25,6 +25,7 @@ import loopy as lp import numpy as np from pymbolic.mapper import IdentityMapper, CSECachingMapperMixin from sumpy.symbolic import pymbolic_real_norm_2 +import sumpy.symbolic as sym from pymbolic.primitives import make_sym_vector from pymbolic import var from collections import defaultdict @@ -379,6 +380,14 @@ class ExpressionKernel(Kernel): from sumpy.tools import ExprDerivativeTaker return ExprDerivativeTaker(self.get_expression(dvec), dvec, rscale, sac) + def get_pde_as_diff_op(self): + r""" + Returns the PDE for the kernel as a + :class:`sumpy.expansion.diff_op.LinearPDESystemOperator` object `L` + where `L(u) = 0` is the PDE. + """ + raise NotImplementedError + one_kernel_2d = ExpressionKernel( dim=2, @@ -431,6 +440,11 @@ class LaplaceKernel(ExpressionKernel): from sumpy.tools import LaplaceDerivativeTaker return LaplaceDerivativeTaker(self.get_expression(dvec), dvec, rscale, sac) + def get_pde_as_diff_op(self): + from sumpy.expansion.diff_op import make_identity_diff_op, laplacian + w = make_identity_diff_op(self.dim) + return laplacian(w) + class BiharmonicKernel(ExpressionKernel): init_arg_names = ("dim",) @@ -475,6 +489,11 @@ class BiharmonicKernel(ExpressionKernel): return RadialDerivativeTaker(self.get_expression(dvec), dvec, rscale, sac) + def get_pde_as_diff_op(self): + from sumpy.expansion.diff_op import make_identity_diff_op, laplacian + w = make_identity_diff_op(self.dim) + return laplacian(laplacian(w)) + class HelmholtzKernel(ExpressionKernel): init_arg_names = ("dim", "helmholtz_k_name", "allow_evanescent") @@ -552,6 +571,13 @@ class HelmholtzKernel(ExpressionKernel): from sumpy.tools import HelmholtzDerivativeTaker return HelmholtzDerivativeTaker(self.get_expression(dvec), dvec, rscale, sac) + def get_pde_as_diff_op(self): + from sumpy.expansion.diff_op import make_identity_diff_op, laplacian + + w = make_identity_diff_op(self.dim) + k = sym.Symbol(self.helmholtz_k_name) + return (laplacian(w) + k**2 * w) + class YukawaKernel(ExpressionKernel): init_arg_names = ("dim", "yukawa_lambda_name") @@ -633,6 +659,12 @@ class YukawaKernel(ExpressionKernel): from sumpy.tools import HelmholtzDerivativeTaker return HelmholtzDerivativeTaker(self.get_expression(dvec), dvec, rscale, sac) + def get_pde_as_diff_op(self): + from sumpy.expansion.diff_op import make_identity_diff_op, laplacian + w = make_identity_diff_op(self.dim) + lam = sym.Symbol(self.yukawa_lambda_name) + return (laplacian(w) - lam**2 * w) + class StokesletKernel(ExpressionKernel): init_arg_names = ("dim", "icomp", "jcomp", "viscosity_mu_name") @@ -700,6 +732,11 @@ class StokesletKernel(ExpressionKernel): mapper_method = "map_stokeslet_kernel" + def get_pde_as_diff_op(self): + from sumpy.expansion.diff_op import make_identity_diff_op, laplacian + w = make_identity_diff_op(self.dim) + return laplacian(laplacian(w)) + class StressletKernel(ExpressionKernel): init_arg_names = ("dim", "icomp", "jcomp", "kcomp", "viscosity_mu_name") @@ -769,6 +806,11 @@ class StressletKernel(ExpressionKernel): mapper_method = "map_stresslet_kernel" + def get_pde_as_diff_op(self): + from sumpy.expansion.diff_op import make_identity_diff_op, laplacian + w = make_identity_diff_op(self.dim) + return laplacian(laplacian(w)) + # }}} diff --git a/sumpy/qbx.py b/sumpy/qbx.py index ff447c315a059af6144f0e19e9bc9be90a3a2e48..f530cd2b9f484efc2e9b0f434e7b41c0cd95aa6c 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -98,17 +98,25 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): return coefficients def _evaluate(self, sac, avec, bvec, rscale, expansion_nr, coefficients): - expn = self.target_kernels[expansion_nr] + from sumpy.expansion.local import LineTaylorLocalExpansion + tgt_knl = self.target_kernels[expansion_nr] + if isinstance(tgt_knl, LineTaylorLocalExpansion): + # In LineTaylorLocalExpansion.evaluate, we can't run + # postprocess_at_target because the coefficients are assigned + # symbols and postprocess with a derivative will make them zero. + # Instead run postprocess here before the coeffients are assigned. + coefficients = [tgt_knl.postprocess_at_target(coeff, bvec) for + coeff in coefficients] + assigned_coeffs = [ sym.Symbol( sac.assign_unique("expn%dcoeff%s" % ( expansion_nr, stringify_expn_index(i)), - expn.postprocess_at_target( - coefficients[self.expansion.get_storage_index(i)], bvec))) + coefficients[self.expansion.get_storage_index(i)])) for i in self.expansion.get_coefficient_identifiers()] return sac.assign_unique("expn%d_result" % expansion_nr, - self.expansion.evaluate(expn, assigned_coeffs, bvec, rscale)) + self.expansion.evaluate(tgt_knl, assigned_coeffs, bvec, rscale)) def get_loopy_insns_and_result_names(self): from sumpy.symbolic import make_sym_vector @@ -348,7 +356,7 @@ class LayerPotentialMatrixGenerator(LayerPotentialBase): # }}} -# {{{ +# {{{ matrix block generator class LayerPotentialMatrixBlockGenerator(LayerPotentialBase): """Generator for a subset of the layer potential matrix entries. diff --git a/test/test_fmm.py b/test/test_fmm.py index fe159b48c4ef268bd2800d641f3c8de52457aefa..cb5465b1d352bd92e4c1a52d201eb967d31af75e 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -31,13 +31,11 @@ from sumpy.kernel import LaplaceKernel, HelmholtzKernel, YukawaKernel from sumpy.expansion.multipole import ( VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion, Y2DMultipoleExpansion, - LaplaceConformingVolumeTaylorMultipoleExpansion, - HelmholtzConformingVolumeTaylorMultipoleExpansion) + LinearPDEConformingVolumeTaylorMultipoleExpansion) from sumpy.expansion.local import ( VolumeTaylorLocalExpansion, H2DLocalExpansion, Y2DLocalExpansion, - LaplaceConformingVolumeTaylorLocalExpansion, - HelmholtzConformingVolumeTaylorLocalExpansion) + LinearPDEConformingVolumeTaylorLocalExpansion) import pytest @@ -58,27 +56,27 @@ else: False), (LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion, True), - (LaplaceKernel(2), LaplaceConformingVolumeTaylorLocalExpansion, - LaplaceConformingVolumeTaylorMultipoleExpansion, True), - (LaplaceKernel(2), LaplaceConformingVolumeTaylorLocalExpansion, - LaplaceConformingVolumeTaylorMultipoleExpansion, False), + (LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, True), + (LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, False), (LaplaceKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion, False), (LaplaceKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion, True), - (LaplaceKernel(3), LaplaceConformingVolumeTaylorLocalExpansion, - LaplaceConformingVolumeTaylorMultipoleExpansion, True), - (LaplaceKernel(3), LaplaceConformingVolumeTaylorLocalExpansion, - LaplaceConformingVolumeTaylorMultipoleExpansion, False), + (LaplaceKernel(3), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, True), + (LaplaceKernel(3), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, False), (HelmholtzKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion, False), - (HelmholtzKernel(2), HelmholtzConformingVolumeTaylorLocalExpansion, - HelmholtzConformingVolumeTaylorMultipoleExpansion, False), + (HelmholtzKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, False), (HelmholtzKernel(2), H2DLocalExpansion, H2DMultipoleExpansion, False), (HelmholtzKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion, False), - (HelmholtzKernel(3), HelmholtzConformingVolumeTaylorLocalExpansion, - HelmholtzConformingVolumeTaylorMultipoleExpansion, False), + (HelmholtzKernel(3), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion, False), (YukawaKernel(2), Y2DLocalExpansion, Y2DMultipoleExpansion, False), ]) def test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class, @@ -232,8 +230,8 @@ def test_unified_single_and_double(ctx_factory): queue = cl.CommandQueue(ctx) knl = LaplaceKernel(2) - local_expn_class = LaplaceConformingVolumeTaylorLocalExpansion - mpole_expn_class = LaplaceConformingVolumeTaylorMultipoleExpansion + local_expn_class = LinearPDEConformingVolumeTaylorLocalExpansion + mpole_expn_class = LinearPDEConformingVolumeTaylorMultipoleExpansion nsources = 1000 ntargets = 300 diff --git a/test/test_kernels.py b/test/test_kernels.py index 7efa44f4ebcc038151b752b529f6518042c1fa6f..cd34b84a3bf544f1ce9a2287cfd7c583a3862356 100644 --- a/test/test_kernels.py +++ b/test/test_kernels.py @@ -32,14 +32,10 @@ from pyopencl.tools import ( # noqa from sumpy.expansion.multipole import ( VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion, VolumeTaylorMultipoleExpansionBase, - LaplaceConformingVolumeTaylorMultipoleExpansion, - HelmholtzConformingVolumeTaylorMultipoleExpansion, - BiharmonicConformingVolumeTaylorMultipoleExpansion) + LinearPDEConformingVolumeTaylorMultipoleExpansion) from sumpy.expansion.local import ( VolumeTaylorLocalExpansion, H2DLocalExpansion, - LaplaceConformingVolumeTaylorLocalExpansion, - HelmholtzConformingVolumeTaylorLocalExpansion, - BiharmonicConformingVolumeTaylorLocalExpansion) + LinearPDEConformingVolumeTaylorLocalExpansion) from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative, DirectionalSourceDerivative, BiharmonicKernel, StokesletKernel) import sumpy.symbolic as sym @@ -107,8 +103,8 @@ def test_p2p(ctx_factory, exclude_self): @pytest.mark.parametrize(("base_knl", "expn_class"), [ - (LaplaceKernel(2), LaplaceConformingVolumeTaylorLocalExpansion), - (LaplaceKernel(2), LaplaceConformingVolumeTaylorMultipoleExpansion), + (LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion), + (LaplaceKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion), ]) def test_p2e_multiple(ctx_factory, base_knl, expn_class): @@ -220,13 +216,13 @@ def test_p2e_multiple(ctx_factory, base_knl, expn_class): @pytest.mark.parametrize(("base_knl", "expn_class"), [ (LaplaceKernel(2), VolumeTaylorLocalExpansion), (LaplaceKernel(2), VolumeTaylorMultipoleExpansion), - (LaplaceKernel(2), LaplaceConformingVolumeTaylorLocalExpansion), - (LaplaceKernel(2), LaplaceConformingVolumeTaylorMultipoleExpansion), + (LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion), + (LaplaceKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), VolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), VolumeTaylorLocalExpansion), - (HelmholtzKernel(2), HelmholtzConformingVolumeTaylorLocalExpansion), - (HelmholtzKernel(2), HelmholtzConformingVolumeTaylorMultipoleExpansion), + (HelmholtzKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion), + (HelmholtzKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), H2DLocalExpansion), (HelmholtzKernel(2), H2DMultipoleExpansion), @@ -238,9 +234,9 @@ def test_p2e_multiple(ctx_factory, base_knl, expn_class): (HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorMultipoleExpansion), (HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorLocalExpansion), (HelmholtzKernel(2, allow_evanescent=True), - HelmholtzConformingVolumeTaylorLocalExpansion), + LinearPDEConformingVolumeTaylorLocalExpansion), (HelmholtzKernel(2, allow_evanescent=True), - HelmholtzConformingVolumeTaylorMultipoleExpansion), + LinearPDEConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2, allow_evanescent=True), H2DLocalExpansion), (HelmholtzKernel(2, allow_evanescent=True), H2DMultipoleExpansion), ]) @@ -457,16 +453,16 @@ def test_p2e2p(ctx_factory, 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), + (LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), - (HelmholtzKernel(2), HelmholtzConformingVolumeTaylorLocalExpansion, - HelmholtzConformingVolumeTaylorMultipoleExpansion), + (HelmholtzKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion), (HelmholtzKernel(2), H2DLocalExpansion, H2DMultipoleExpansion), (StokesletKernel(2, 0, 0), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), - (StokesletKernel(2, 0, 0), BiharmonicConformingVolumeTaylorLocalExpansion, - BiharmonicConformingVolumeTaylorMultipoleExpansion), + (StokesletKernel(2, 0, 0), LinearPDEConformingVolumeTaylorLocalExpansion, + LinearPDEConformingVolumeTaylorMultipoleExpansion), ]) def test_translations(ctx_factory, knl, local_expn_class, mpole_expn_class): logging.basicConfig(level=logging.INFO) diff --git a/test/test_qbx.py b/test/test_qbx.py index 11f7da7b65d3dde2e63da2cf1f29e91718cddc95..379ef0e221327f13d5452c34e8f1601de9fad421 100644 --- a/test/test_qbx.py +++ b/test/test_qbx.py @@ -29,18 +29,20 @@ from pyopencl.tools import ( # noqa import logging logger = logging.getLogger(__name__) +from sumpy.expansion.local import ( + LineTaylorLocalExpansion, VolumeTaylorLocalExpansion) +import pytest -try: - import faulthandler -except ImportError: - pass -else: - faulthandler.enable() +@pytest.mark.parametrize("expn_class", [ + LineTaylorLocalExpansion, + VolumeTaylorLocalExpansion, + ]) +def test_direct_qbx_vs_eigval(ctx_factory, expn_class): + """This evaluates a single layer potential on a circle using a known + eigenvalue/eigenvector combination. + """ - -def test_direct(ctx_factory): - # This evaluates a single layer potential on a circle. logging.basicConfig(level=logging.INFO) ctx = ctx_factory() @@ -52,9 +54,8 @@ def test_direct(ctx_factory): order = 12 from sumpy.qbx import LayerPotential - from sumpy.expansion.local import LineTaylorLocalExpansion - lpot = LayerPotential(ctx, expansion=LineTaylorLocalExpansion(lknl, order), + lpot = LayerPotential(ctx, expansion=expn_class(lknl, order), target_kernels=(lknl,), source_kernels=(lknl,)) mode_nr = 25 @@ -95,8 +96,77 @@ def test_direct(ctx_factory): assert eocrec.order_estimate() > order - slack -# You can test individual routines by typing -# $ python test_kernels.py 'test_p2p(cl.create_some_context)' +@pytest.mark.parametrize("expn_class", [ + LineTaylorLocalExpansion, + VolumeTaylorLocalExpansion, + ]) +def test_direct_qbx_vs_eigval_with_tgt_deriv(ctx_factory, expn_class): + """This evaluates a single layer potential on a circle using a known + eigenvalue/eigenvector combination. + """ + + logging.basicConfig(level=logging.INFO) + + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + from sumpy.kernel import LaplaceKernel, AxisTargetDerivative + lknl = LaplaceKernel(2) + + order = 8 + + from sumpy.qbx import LayerPotential + + lpot_dx = LayerPotential(ctx, expansion=expn_class(lknl, order), + target_kernels=(AxisTargetDerivative(0, lknl),), source_kernels=(lknl,)) + lpot_dy = LayerPotential(ctx, expansion=expn_class(lknl, order), + target_kernels=(AxisTargetDerivative(1, lknl),), source_kernels=(lknl,)) + + mode_nr = 15 + + from pytools.convergence import EOCRecorder + + eocrec = EOCRecorder() + + for n in [200, 300, 400]: + t = np.linspace(0, 2 * np.pi, n, endpoint=False) + unit_circle = np.exp(1j * t) + unit_circle = np.array([unit_circle.real, unit_circle.imag]) + + sigma = np.cos(mode_nr * t) + #eigval = 1/(2*mode_nr) + eigval = 0.5 + + result_ref = eigval * sigma + + h = 2 * np.pi / n + + targets = unit_circle + sources = unit_circle + + radius = 7 * h + centers = unit_circle * (1 - radius) + + expansion_radii = np.ones(n) * radius + + strengths = (sigma * h,) + + evt, (result_qbx_dx,) = lpot_dx(queue, targets, sources, centers, strengths, + expansion_radii=expansion_radii) + evt, (result_qbx_dy,) = lpot_dy(queue, targets, sources, centers, strengths, + expansion_radii=expansion_radii) + + normals = unit_circle + result_qbx = normals[0] * result_qbx_dx + normals[1] * result_qbx_dy + + eocrec.add_data_point(h, np.max(np.abs(result_ref - result_qbx))) + + if expn_class is not LineTaylorLocalExpansion: + print(eocrec) + + slack = 1.5 + assert eocrec.order_estimate() > order - slack + if __name__ == "__main__": if len(sys.argv) > 1: