diff --git a/examples/expansion-toys.py b/examples/expansion-toys.py index 782bbcf51f0b539531482d196eb0e55d06c0b4e3..78485ce89a875654f0226f5d3edb882251083934 100644 --- a/examples/expansion-toys.py +++ b/examples/expansion-toys.py @@ -6,8 +6,14 @@ import matplotlib.pyplot as plt def main(): - from sumpy.kernel import LaplaceKernel - tctx = t.ToyContext(cl.create_some_context(), LaplaceKernel(2)) + from sumpy.kernel import ( # noqa: F401 + YukawaKernel, HelmholtzKernel, LaplaceKernel) + tctx = t.ToyContext( + cl.create_some_context(), + #LaplaceKernel(2), + YukawaKernel(2), extra_source_kwargs={"lam": 5}, + #HelmholtzKernel(2), extra_source_kwargs={"k": 0.3}, + ) pt_src = t.PointSources( tctx, @@ -21,13 +27,13 @@ def main(): plt.colorbar() plt.show() - mexp = t.multipole_expand(pt_src, [0, 0], 9) + mexp = t.multipole_expand(pt_src, [0, 0], 5) mexp2 = t.multipole_expand(mexp, [0, 0.25]) lexp = t.local_expand(mexp, [3, 0]) - lexp2 = t.local_expand(lexp, [3, 1]) + lexp2 = t.local_expand(lexp, [3, 1], 3) - diff = mexp - pt_src - diff = mexp2 - pt_src + #diff = mexp - pt_src + #diff = mexp2 - pt_src diff = lexp2 - pt_src print(t.l_inf(diff, 1.2, center=lexp2.center)) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index f50bd5fbb9facd732c0c192c787a0a420d72b58c..b5415b1e60250a732bcc8af141960fd57b33b9c0 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -32,6 +32,13 @@ from sumpy.tools import MiDerivativeTaker __doc__ = """ .. autoclass:: ExpansionBase + +Expansion Factories +^^^^^^^^^^^^^^^^^^^ + +.. autoclass:: ExpansionFactoryBase +.. autoclass:: DefaultExpansionFactory +.. autoclass:: VolumeTaylorExpansionFactory """ logger = logging.getLogger(__name__) @@ -434,4 +441,89 @@ class HelmholtzConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase): # }}} +# {{{ expansion factory + +class ExpansionFactoryBase(object): + """An interface + .. automethod:: get_local_expansion_class + .. automethod:: get_multipole_expansion_class + """ + + def get_local_expansion_class(self, base_kernel): + """Returns a subclass of :class:`ExpansionBase` suitable for *base_kernel*. + """ + raise NotImplementedError() + + def get_multipole_expansion_class(self, base_kernel): + """Returns a subclass of :class:`ExpansionBase` suitable for *base_kernel*. + """ + raise NotImplementedError() + + +class VolumeTaylorExpansionFactory(ExpansionFactoryBase): + """An implementation of :class:`ExpansionFactoryBase` that uses Volume Taylor + expansions for each kernel. + """ + + def get_local_expansion_class(self, base_kernel): + """Returns a subclass of :class:`ExpansionBase` suitable for *base_kernel*. + """ + from sumpy.expansion.local import VolumeTaylorLocalExpansion + return VolumeTaylorLocalExpansion + + def get_multipole_expansion_class(self, base_kernel): + """Returns a subclass of :class:`ExpansionBase` suitable for *base_kernel*. + """ + from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion + return VolumeTaylorMultipoleExpansion + + +class DefaultExpansionFactory(ExpansionFactoryBase): + """An implementation of :class:`ExpansionFactoryBase` that gives the 'best known' + expansion for each kernel. + """ + + 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 + if (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel) + and base_kernel.dim == 2): + from sumpy.expansion.local import H2DLocalExpansion + return H2DLocalExpansion + elif (isinstance(base_kernel.get_base_kernel(), YukawaKernel) + and base_kernel.dim == 2): + from sumpy.expansion.local import Y2DLocalExpansion + return Y2DLocalExpansion + elif isinstance(base_kernel.get_base_kernel(), LaplaceKernel): + from sumpy.expansion.local import \ + LaplaceConformingVolumeTaylorLocalExpansion + return LaplaceConformingVolumeTaylorLocalExpansion + else: + from sumpy.expansion.local import VolumeTaylorLocalExpansion + 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 + if (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel) + and base_kernel.dim == 2): + from sumpy.expansion.multipole import H2DMultipoleExpansion + return H2DMultipoleExpansion + elif (isinstance(base_kernel.get_base_kernel(), YukawaKernel) + and base_kernel.dim == 2): + from sumpy.expansion.multipole import Y2DMultipoleExpansion + return Y2DMultipoleExpansion + elif isinstance(base_kernel.get_base_kernel(), LaplaceKernel): + from sumpy.expansion.multipole import ( + LaplaceConformingVolumeTaylorMultipoleExpansion) + return LaplaceConformingVolumeTaylorMultipoleExpansion + else: + from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion + return VolumeTaylorMultipoleExpansion + +# }}} + + # vim: fdm=marker diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 02a47fed6abdb87de3097e7c5261475693c936cd..76f3d1ebaaa6a12037d333559deb49ffc8c19e48 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -41,6 +41,7 @@ __doc__ = """ .. autoclass:: VolumeTaylorLocalExpansion .. autoclass:: H2DLocalExpansion +.. autoclass:: Y2DLocalExpansion .. autoclass:: LineTaylorLocalExpansion """ @@ -202,16 +203,9 @@ class HelmholtzConformingVolumeTaylorLocalExpansion( # }}} -# {{{ 2D J-expansion - -class H2DLocalExpansion(LocalExpansionBase): - def __init__(self, kernel, order): - from sumpy.kernel import HelmholtzKernel - assert (isinstance(kernel.get_base_kernel(), HelmholtzKernel) - and kernel.dim == 2) - - LocalExpansionBase.__init__(self, kernel, order) +# {{{ 2D Bessel-based-expansion +class _FourierBesselLocalExpansion(LocalExpansionBase): def get_storage_index(self, k): return self.order+k @@ -222,13 +216,13 @@ class H2DLocalExpansion(LocalExpansionBase): from sumpy.symbolic import sym_real_norm_2 hankel_1 = sym.Function("hankel_1") - k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + arg_scale = self.get_bessel_arg_scaling() # The coordinates are negated since avec points from source to center. source_angle_rel_center = sym.atan2(-avec[1], -avec[0]) avec_len = sym_real_norm_2(avec) return [self.kernel.postprocess_at_source( - hankel_1(l, k * avec_len) + hankel_1(l, arg_scale * avec_len) * sym.exp(sym.I * l * source_angle_rel_center), avec) for l in self.get_coefficient_identifiers()] @@ -238,20 +232,20 @@ class H2DLocalExpansion(LocalExpansionBase): bvec_len = sym_real_norm_2(bvec) target_angle_rel_center = sym.atan2(bvec[1], bvec[0]) - k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + arg_scale = self.get_bessel_arg_scaling() return sum(coeffs[self.get_storage_index(l)] * self.kernel.postprocess_at_target( - bessel_j(l, k * bvec_len) + bessel_j(l, arg_scale * bvec_len) * sym.exp(sym.I * l * -target_angle_rel_center), bvec) for l in self.get_coefficient_identifiers()) def translate_from(self, src_expansion, src_coeff_exprs, dvec): from sumpy.symbolic import sym_real_norm_2 - k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + arg_scale = self.get_bessel_arg_scaling() - if isinstance(src_expansion, H2DLocalExpansion): + if isinstance(src_expansion, type(self)): dvec_len = sym_real_norm_2(dvec) bessel_j = sym.Function("bessel_j") new_center_angle_rel_old_center = sym.atan2(dvec[1], dvec[0]) @@ -259,28 +253,57 @@ class H2DLocalExpansion(LocalExpansionBase): for l in self.get_coefficient_identifiers(): translated_coeffs.append( sum(src_coeff_exprs[src_expansion.get_storage_index(m)] - * bessel_j(m - l, k * dvec_len) + * bessel_j(m - l, arg_scale * dvec_len) * sym.exp(sym.I * (m - l) * -new_center_angle_rel_old_center) for m in src_expansion.get_coefficient_identifiers())) return translated_coeffs - from sumpy.expansion.multipole import H2DMultipoleExpansion - if isinstance(src_expansion, H2DMultipoleExpansion): + if isinstance(src_expansion, self.mpole_expn_class): dvec_len = sym_real_norm_2(dvec) hankel_1 = sym.Function("hankel_1") new_center_angle_rel_old_center = sym.atan2(dvec[1], dvec[0]) translated_coeffs = [] for l in self.get_coefficient_identifiers(): translated_coeffs.append( - sum((-1) ** l * hankel_1(m + l, k * dvec_len) + sum((-1) ** l * hankel_1(m + l, arg_scale * dvec_len) * sym.exp(sym.I * (m + l) * new_center_angle_rel_old_center) * src_coeff_exprs[src_expansion.get_storage_index(m)] for m in src_expansion.get_coefficient_identifiers())) return translated_coeffs - raise RuntimeError("do not know how to translate %s to " - "local 2D Helmholtz Bessel expansion" - % type(src_expansion).__name__) + raise RuntimeError("do not know how to translate %s to %s" + % (type(src_expansion).__name__, + type(self).__name__)) + + +class H2DLocalExpansion(_FourierBesselLocalExpansion): + def __init__(self, kernel, order): + from sumpy.kernel import HelmholtzKernel + assert (isinstance(kernel.get_base_kernel(), HelmholtzKernel) + and kernel.dim == 2) + + super(H2DLocalExpansion, self).__init__(kernel, order) + + from sumpy.expansion.multipole import H2DMultipoleExpansion + self.mpole_expn_class = H2DMultipoleExpansion + + def get_bessel_arg_scaling(self): + return sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + + +class Y2DLocalExpansion(_FourierBesselLocalExpansion): + def __init__(self, kernel, order): + from sumpy.kernel import YukawaKernel + assert (isinstance(kernel.get_base_kernel(), YukawaKernel) + and kernel.dim == 2) + + super(Y2DLocalExpansion, self).__init__(kernel, order) + + from sumpy.expansion.multipole import Y2DMultipoleExpansion + self.mpole_expn_class = Y2DMultipoleExpansion + + def get_bessel_arg_scaling(self): + return sym.I * sym.Symbol(self.kernel.get_base_kernel().yukawa_lambda_name) # }}} diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index ce3b78d8ae25049f48f1b1569ac7483bf8586923..958bca6dbec6b8fdcbbb77545f039ae92a65d3d1 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -37,6 +37,7 @@ __doc__ = """ .. autoclass:: VolumeTaylorMultipoleExpansion .. autoclass:: H2DMultipoleExpansion +.. autoclass:: Y2DMultipoleExpansion """ @@ -185,16 +186,9 @@ class HelmholtzConformingVolumeTaylorMultipoleExpansion( # }}} -# {{{ 2D H-expansion - -class H2DMultipoleExpansion(MultipoleExpansionBase): - def __init__(self, kernel, order): - from sumpy.kernel import HelmholtzKernel - assert (isinstance(kernel.get_base_kernel(), HelmholtzKernel) - and kernel.dim == 2) - - MultipoleExpansionBase.__init__(self, kernel, order) +# {{{ 2D Hankel-based expansions +class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): def get_storage_index(self, k): return self.order+k @@ -206,12 +200,12 @@ class H2DMultipoleExpansion(MultipoleExpansionBase): bessel_j = sym.Function("bessel_j") avec_len = sym_real_norm_2(avec) - k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + arg_scale = self.get_bessel_arg_scaling() # The coordinates are negated since avec points from source to center. source_angle_rel_center = sym.atan2(-avec[1], -avec[0]) return [self.kernel.postprocess_at_source( - bessel_j(l, k * avec_len) * + bessel_j(l, arg_scale * avec_len) * sym.exp(sym.I * l * -source_angle_rel_center), avec) for l in self.get_coefficient_identifiers()] @@ -221,35 +215,59 @@ class H2DMultipoleExpansion(MultipoleExpansionBase): bvec_len = sym_real_norm_2(bvec) target_angle_rel_center = sym.atan2(bvec[1], bvec[0]) - k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + arg_scale = self.get_bessel_arg_scaling() return sum(coeffs[self.get_storage_index(l)] * self.kernel.postprocess_at_target( - hankel_1(l, k * bvec_len) + hankel_1(l, arg_scale * bvec_len) * sym.exp(sym.I * l * target_angle_rel_center), bvec) for l in self.get_coefficient_identifiers()) def translate_from(self, src_expansion, src_coeff_exprs, dvec): if not isinstance(src_expansion, type(self)): - raise RuntimeError("do not know how to translate %s to " - "multipole 2D Helmholtz Bessel expansion" - % type(src_expansion).__name__) + raise RuntimeError("do not know how to translate %s to %s" + % (type(src_expansion).__name__, + type(self).__name__)) from sumpy.symbolic import sym_real_norm_2 dvec_len = sym_real_norm_2(dvec) bessel_j = sym.Function("bessel_j") new_center_angle_rel_old_center = sym.atan2(dvec[1], dvec[0]) - k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + arg_scale = self.get_bessel_arg_scaling() translated_coeffs = [] for l in self.get_coefficient_identifiers(): translated_coeffs.append( sum(src_coeff_exprs[src_expansion.get_storage_index(m)] - * bessel_j(m - l, k * dvec_len) + * bessel_j(m - l, arg_scale * dvec_len) * sym.exp(sym.I * (m - l) * new_center_angle_rel_old_center) for m in src_expansion.get_coefficient_identifiers())) return translated_coeffs + +class H2DMultipoleExpansion(_HankelBased2DMultipoleExpansion): + def __init__(self, kernel, order): + from sumpy.kernel import HelmholtzKernel + assert (isinstance(kernel.get_base_kernel(), HelmholtzKernel) + and kernel.dim == 2) + + super(H2DMultipoleExpansion, self).__init__(kernel, order) + + def get_bessel_arg_scaling(self): + return sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + + +class Y2DMultipoleExpansion(_HankelBased2DMultipoleExpansion): + def __init__(self, kernel, order): + from sumpy.kernel import YukawaKernel + assert (isinstance(kernel.get_base_kernel(), YukawaKernel) + and kernel.dim == 2) + + super(Y2DMultipoleExpansion, self).__init__(kernel, order) + + def get_bessel_arg_scaling(self): + return sym.I * sym.Symbol(self.kernel.get_base_kernel().yukawa_lambda_name) + # }}} # vim: fdm=marker diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 735ac0f6b0d7fb285ce81b8467f64b170b0419db..118e9f77367a184aa369ee59ce782693364aae85 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -49,6 +49,7 @@ PDE kernels .. autoclass:: LaplaceKernel .. autoclass:: HelmholtzKernel +.. autoclass:: YukawaKernel .. autoclass:: StokesletKernel .. autoclass:: StressletKernel @@ -72,7 +73,6 @@ Transforming kernels .. autoclass:: AxisTargetDerivativeRemover .. autoclass:: TargetDerivativeRemover .. autoclass:: DerivativeCounter -.. autoclass:: KernelDimensionSetter """ @@ -323,9 +323,6 @@ class LaplaceKernel(ExpressionKernel): r = pymbolic_real_norm_2(make_sym_vector("d", dim)) expr = 1/r scaling = 1/(4*var("pi")) - elif dim is None: - expr = None - scaling = None else: raise RuntimeError("unsupported dimensionality") @@ -403,9 +400,6 @@ class HelmholtzKernel(ExpressionKernel): r = pymbolic_real_norm_2(make_sym_vector("d", dim)) expr = var("exp")(var("I")*k*r)/r scaling = 1/(4*var("pi")) - elif dim is None: - expr = None - scaling = None else: raise RuntimeError("unsupported dimensionality") @@ -429,11 +423,8 @@ class HelmholtzKernel(ExpressionKernel): self.allow_evanescent)) def __repr__(self): - if self._dim is not None: - return "HelmKnl%dD(%s)" % ( - self.dim, self.helmholtz_k_name) - else: - return "HelmKnl(%s)" % (self.helmholtz_k_name) + return "HelmKnl%dD(%s)" % ( + self.dim, self.helmholtz_k_name) def prepare_loopy_kernel(self, loopy_knl): from sumpy.codegen import (bessel_preamble_generator, bessel_mangler) @@ -458,6 +449,65 @@ class HelmholtzKernel(ExpressionKernel): mapper_method = "map_helmholtz_kernel" +class YukawaKernel(ExpressionKernel): + init_arg_names = ("dim", "yukawa_lambda_name") + + def __init__(self, dim=None, yukawa_lambda_name="lam"): + """ + :arg yukawa_lambda_name: The argument name to use for the Yukawa + parameter when generating functions to evaluate this kernel. + """ + lam = var(yukawa_lambda_name) + + if dim == 2: + r = pymbolic_real_norm_2(make_sym_vector("d", dim)) + + # http://dlmf.nist.gov/10.27#E8 + expr = var("hankel_1")(0, var("I")*lam*r) + scaling_for_K0 = 1/2*var("pi") # noqa: N806 + + scaling = -1/(2*var("pi")) * scaling_for_K0 + else: + raise RuntimeError("unsupported dimensionality") + + ExpressionKernel.__init__( + self, + dim, + expression=expr, + scaling=scaling, + is_complex_valued=True) + + self.yukawa_lambda_name = yukawa_lambda_name + + def __getinitargs__(self): + return (self._dim, self.yukawa_lambda_name) + + def update_persistent_hash(self, key_hash, key_builder): + key_hash.update(type(self).__name__.encode("utf8")) + key_builder.rec(key_hash, (self.dim, self.yukawa_lambda_name)) + + def __repr__(self): + return "YukKnl%dD(%s)" % ( + self.dim, self.yukawa_lambda_name) + + def prepare_loopy_kernel(self, loopy_knl): + from sumpy.codegen import (bessel_preamble_generator, bessel_mangler) + loopy_knl = lp.register_function_manglers(loopy_knl, + [bessel_mangler]) + loopy_knl = lp.register_preamble_generators(loopy_knl, + [bessel_preamble_generator]) + + return loopy_knl + + def get_args(self): + return [ + KernelArgument( + loopy_arg=lp.ValueArg(self.yukawa_lambda_name, np.float64), + )] + + mapper_method = "map_yukawa_kernel" + + class StokesletKernel(ExpressionKernel): init_arg_names = ("dim", "icomp", "jcomp", "viscosity_mu_name") @@ -827,6 +877,7 @@ class KernelIdentityMapper(KernelMapper): map_laplace_kernel = map_expression_kernel map_biharmonic_kernel = map_expression_kernel map_helmholtz_kernel = map_expression_kernel + map_yukawa_kernel = map_expression_kernel map_stokeslet_kernel = map_expression_kernel map_stresslet_kernel = map_expression_kernel @@ -861,6 +912,7 @@ class DerivativeCounter(KernelCombineMapper): map_laplace_kernel = map_expression_kernel map_biharmonic_kernel = map_expression_kernel map_helmholtz_kernel = map_expression_kernel + map_yukawa_kernel = map_expression_kernel map_stokeslet_kernel = map_expression_kernel map_stresslet_kernel = map_expression_kernel @@ -870,63 +922,6 @@ class DerivativeCounter(KernelCombineMapper): map_directional_target_derivative = map_axis_target_derivative map_directional_source_derivative = map_axis_target_derivative - -class KernelDimensionSetter(KernelIdentityMapper): - """Deprecated: This is no longer used and will be removed in 2018. - """ - - def __init__(self, dim): - self.dim = dim - - def map_expression_kernel(self, kernel): - if kernel._dim is not None and self.dim != kernel.dim: - raise RuntimeError("cannot set kernel dimension to new value (%d)" - "different from existing one (%d)" - % (self.dim, kernel.dim)) - - return kernel - - def map_laplace_kernel(self, kernel): - if kernel._dim is not None and self.dim != kernel.dim: - raise RuntimeError("cannot set kernel dimension to new value (%d)" - "different from existing one (%d)" - % (self.dim, kernel.dim)) - - return LaplaceKernel(self.dim) - - def map_helmholtz_kernel(self, kernel): - if kernel._dim is not None and self.dim != kernel.dim: - raise RuntimeError("cannot set kernel dimension to new value (%d)" - "different from existing one (%d)" - % (self.dim, kernel.dim)) - - return HelmholtzKernel(self.dim, - helmholtz_k_name=kernel.helmholtz_k_name, - allow_evanescent=kernel.allow_evanescent) - - def map_stokeslet_kernel(self, kernel): - if kernel._dim is not None and self.dim != kernel.dim: - raise RuntimeError("cannot set kernel dimension to new value (%d)" - "different from existing one (%d)" - % (self.dim, kernel.dim)) - - return StokesletKernel(self.dim, - kernel.icomp, - kernel.jcomp, - viscosity_mu_name=kernel.viscosity_mu_name) - - def map_stresslet_kernel(self, kernel): - if kernel._dim is not None and self.dim != kernel.dim: - raise RuntimeError("cannot set kernel dimension to new value (%d)" - "different from existing one (%d)" - % (self.dim, kernel.dim)) - - return StressletKernel(self.dim, - kernel.icomp, - kernel.jcomp, - kernel.kcomp, - viscosity_mu_name=kernel.viscosity_mu_name) - # }}} diff --git a/sumpy/toys.py b/sumpy/toys.py index 08330406e6f00f4fa3a0a60dbb4f37cd8142424c..d61af4faa80b1e40a7015f3b0c4db9ed4b00e231 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -42,17 +42,21 @@ class ToyContext(object): def __init__(self, cl_context, kernel, mpole_expn_class=None, local_expn_class=None, + expansion_factory=None, extra_source_kwargs=None): self.cl_context = cl_context self.queue = cl.CommandQueue(self.cl_context) self.kernel = kernel + if expansion_factory is None: + from sumpy.expansion import DefaultExpansionFactory + expansion_factory = DefaultExpansionFactory() if mpole_expn_class is None: - from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion - mpole_expn_class = VolumeTaylorMultipoleExpansion + mpole_expn_class = \ + expansion_factory.get_multipole_expansion_class(kernel) if local_expn_class is None: - from sumpy.expansion.local import VolumeTaylorLocalExpansion - local_expn_class = VolumeTaylorLocalExpansion + local_expn_class = \ + expansion_factory.get_local_expansion_class(kernel) if extra_source_kwargs is None: extra_source_kwargs = {} @@ -294,7 +298,7 @@ class PointSources(PotentialSource): def eval(self, targets): evt, (potential,) = self.toy_ctx.get_p2p()( self.toy_ctx.queue, targets, self.points, [self.weights], - out_host=True) + out_host=True, **self.toy_ctx.extra_source_kwargs) return potential diff --git a/test/test_fmm.py b/test/test_fmm.py index ee97996896ec9e695d82126fefc7bbb3cffb71b0..c85759fef3f2c6bf6457cb55185ed3f15d32f369 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -30,13 +30,15 @@ import numpy.linalg as la import pyopencl as cl from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) -from sumpy.kernel import LaplaceKernel, HelmholtzKernel +from sumpy.kernel import LaplaceKernel, HelmholtzKernel, YukawaKernel from sumpy.expansion.multipole import ( - VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion, + VolumeTaylorMultipoleExpansion, + H2DMultipoleExpansion, Y2DMultipoleExpansion, LaplaceConformingVolumeTaylorMultipoleExpansion, HelmholtzConformingVolumeTaylorMultipoleExpansion) from sumpy.expansion.local import ( - VolumeTaylorLocalExpansion, H2DLocalExpansion, + VolumeTaylorLocalExpansion, + H2DLocalExpansion, Y2DLocalExpansion, LaplaceConformingVolumeTaylorLocalExpansion, HelmholtzConformingVolumeTaylorLocalExpansion) @@ -101,6 +103,7 @@ def test_level_to_order_lookup(ctx_getter, lookup_func, extra_args): (HelmholtzKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), (HelmholtzKernel(3), HelmholtzConformingVolumeTaylorLocalExpansion, HelmholtzConformingVolumeTaylorMultipoleExpansion), + (YukawaKernel(2), Y2DLocalExpansion, Y2DMultipoleExpansion), ]) def test_sumpy_fmm(ctx_getter, knl, local_expn_class, mpole_expn_class): logging.basicConfig(level=logging.INFO) @@ -187,6 +190,15 @@ def test_sumpy_fmm(ctx_getter, knl, local_expn_class, mpole_expn_class): elif knl.dim == 2 and issubclass(local_expn_class, H2DLocalExpansion): order_values = [10, 12] + elif isinstance(knl, YukawaKernel): + extra_kwargs["lam"] = 2 + dtype = np.complex128 + + if knl.dim == 3: + order_values = [1, 2] + elif knl.dim == 2 and issubclass(local_expn_class, Y2DLocalExpansion): + order_values = [10, 12] + from functools import partial for order in order_values: out_kernels = [knl] diff --git a/test/test_misc.py b/test/test_misc.py index 301708043418ebd00da604d845ae8692418fa00e..9edd1966b279ef6655bcb92d548205cc94a48945 100644 --- a/test/test_misc.py +++ b/test/test_misc.py @@ -32,11 +32,44 @@ import pyopencl as cl # noqa: F401 from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) +from sumpy.kernel import BiharmonicKernel, YukawaKernel -@pytest.mark.parametrize("dim", [2, 3]) -def test_pde_check_biharmonic(ctx_factory, dim, order=5): - from sumpy.kernel import BiharmonicKernel - tctx = t.ToyContext(ctx_factory(), BiharmonicKernel(dim)) + +# {{{ pde check for kernels + +class BiharmonicKernelInfo: + def __init__(self, dim): + self.kernel = BiharmonicKernel(dim) + self.extra_kwargs = {} + + @staticmethod + def pde_func(cp, pot): + return cp.laplace(cp.laplace(pot)) + + nderivs = 4 + + +class YukawaKernelInfo: + def __init__(self, dim, lam): + self.kernel = YukawaKernel(dim) + self.lam = lam + self.extra_kwargs = {"lam": lam} + + def pde_func(self, cp, pot): + return cp.laplace(pot) - self.lam**2*pot + + nderivs = 2 + + +@pytest.mark.parametrize("knl_info", [ + BiharmonicKernelInfo(2), + BiharmonicKernelInfo(3), + YukawaKernelInfo(2, 5), + ]) +def test_pde_check_kernels(ctx_factory, knl_info, order=5): + dim = knl_info.kernel.dim + tctx = t.ToyContext(ctx_factory(), knl_info.kernel, + extra_source_kwargs=knl_info.extra_kwargs) pt_src = t.PointSources( tctx, @@ -51,13 +84,15 @@ def test_pde_check_biharmonic(ctx_factory, dim, order=5): cp = CalculusPatch(np.array([1, 0, 0])[:dim], h=h, order=order) pot = pt_src.eval(cp.points) - deriv = cp.laplace(cp.laplace(pot)) + pde = knl_info.pde_func(cp, pot) - err = la.norm(deriv) + err = la.norm(pde) eoc_rec.add_data_point(h, err) print(eoc_rec) - assert eoc_rec.order_estimate() > order-3-0.1 + assert eoc_rec.order_estimate() > order - knl_info.nderivs + 1 - 0.1 + +# }}} @pytest.mark.parametrize("dim", [1, 2, 3])