diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py index 55c370563a1ea71a4da73334c2f5d0d82b826560..8cb016222d3bae0b35b9ae0a5dab58f2a2f2766e 100644 --- a/examples/coeff-scaling-experiment.py +++ b/examples/coeff-scaling-experiment.py @@ -13,13 +13,17 @@ def main(): from sumpy.kernel import ( # noqa: F401 YukawaKernel, HelmholtzKernel, LaplaceKernel, + BiharmonicKernel, StokesletKernel, StressletKernel, AxisTargetDerivative) tctx = t.ToyContext( cl.create_some_context(), - #LaplaceKernel(dim), + LaplaceKernel(dim), #AxisTargetDerivative(0, LaplaceKernel(dim)), #YukawaKernel(dim), extra_source_kwargs={"lam": 5}, - HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3}, + #HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3}, + #BiharmonicKernel(dim), + #StokesletKernel(dim, 1, 1), extra_source_kwargs={"mu": 0.3}, + #StressletKernel(dim, 1, 1, 0), extra_source_kwargs={"mu": 0.3}, ) np.random.seed(12) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 3507509f3cdfa1430200ef8957d7828ee3e9f97a..f2aec7060db5938f46eb2245a07e591a55b49b3e 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -217,6 +217,9 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): return list(range(-self.order, self.order+1)) def coefficients_from_source(self, avec, bvec, rscale): + if not self.kernel.supports_rscale: + rscale = 1 + from sumpy.symbolic import sym_real_norm_2 hankel_1 = sym.Function("hankel_1") @@ -232,6 +235,9 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): for l in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec, rscale): + if not self.kernel.supports_rscale: + rscale = 1 + from sumpy.symbolic import sym_real_norm_2 bessel_j = sym.Function("bessel_j") bvec_len = sym_real_norm_2(bvec) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 20e1cb7d6b1980ad488237de0d0bebe7af73cc2b..ca09139733a927a373ed43e75da347cbdb37d5cd 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -60,6 +60,9 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): from sumpy.tools import mi_power, mi_factorial + if not self.kernel.supports_rscale: + rscale = 1 + avec = avec/rscale if isinstance(kernel, DirectionalSourceDerivative): @@ -94,6 +97,9 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): self.derivative_wrangler.get_stored_mpole_coefficients_from_full(result)) def evaluate(self, coeffs, bvec, rscale): + if not self.kernel.supports_rscale: + rscale = 1 + taker = self.get_kernel_derivative_taker(bvec, rscale) result = sym.Add(*tuple( coeff diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 32123dae9d8c0b6d3551d4538fb0fcae65049400..4ae60dbba7cdb65ba20b9b330d5e51f1b61a4d77 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -32,8 +32,6 @@ from pymbolic.primitives import make_sym_vector from pymbolic import var -# FIXME: Rename ProxyExpressionKernel back to ExpressionKernel once we're done. - __doc__ = """ Kernel interface ---------------- @@ -102,6 +100,7 @@ class Kernel(object): """Basic kernel interface. .. attribute:: is_complex_valued + .. attribute:: supports_rscale .. attribute:: dim .. automethod:: get_base_kernel @@ -116,6 +115,8 @@ class Kernel(object): .. automethod:: get_source_args """ + supports_rscale = False + def __init__(self, dim=None): self.dim = dim @@ -276,11 +277,6 @@ class Kernel(object): class ExpressionKernel(Kernel): - def get_expression(self, dist_vec): - raise NotImplementedError - - -class ProxyExpressionKernel(Kernel): is_complex_valued = False init_arg_names = ("dim", "expression", "kernel_scaling", "is_complex_valued") @@ -352,12 +348,12 @@ class ProxyExpressionKernel(Kernel): mapper_method = "map_expression_kernel" -one_kernel_2d = ProxyExpressionKernel( +one_kernel_2d = ExpressionKernel( dim=2, proxy_expression=1, global_scaling_const=1, is_complex_valued=False) -one_kernel_3d = ProxyExpressionKernel( +one_kernel_3d = ExpressionKernel( dim=3, proxy_expression=1, global_scaling_const=1, @@ -366,8 +362,9 @@ one_kernel_3d = ProxyExpressionKernel( # {{{ PDE kernels -class LaplaceKernel(ProxyExpressionKernel): +class LaplaceKernel(ExpressionKernel): init_arg_names = ("dim",) + supports_rscale = True def __init__(self, dim=None): # See (Kress LIE, Thm 6.2) for scaling @@ -413,6 +410,7 @@ class LaplaceKernel(ProxyExpressionKernel): class BiharmonicKernel(ExpressionKernel): init_arg_names = ("dim",) + supports_rscale = False def __init__(self, dim=None): r = pymbolic_real_norm_2(make_sym_vector("d", dim)) @@ -426,13 +424,16 @@ class BiharmonicKernel(ExpressionKernel): else: raise RuntimeError("unsupported dimensionality") - ExpressionKernel.__init__( - self, + super(BiharmonicKernel, self).__init__( dim, - expression=expr, - scaling=scaling, + proxy_expression=expr, + global_scaling_const=scaling, is_complex_valued=False) + def adjust_proxy_expression(self, expr, rscale, nderivatives): + assert rscale == 1 + return expr + def __getinitargs__(self): return (self.dim,) @@ -442,8 +443,9 @@ class BiharmonicKernel(ExpressionKernel): mapper_method = "map_biharmonic_kernel" -class HelmholtzKernel(ProxyExpressionKernel): +class HelmholtzKernel(ExpressionKernel): init_arg_names = ("dim", "helmholtz_k_name", "allow_evanescent") + supports_rscale = True def __init__(self, dim=None, helmholtz_k_name="k", allow_evanescent=False): @@ -484,8 +486,10 @@ class HelmholtzKernel(ProxyExpressionKernel): if self.dim == 2: return result - else: + elif self.dim == 3: return result / rscale**(nderivatives+1) + else: + raise RuntimeError("unsupported dimensionality") def __getinitargs__(self): return (self.dim, self.helmholtz_k_name, @@ -525,6 +529,7 @@ class HelmholtzKernel(ProxyExpressionKernel): class YukawaKernel(ExpressionKernel): init_arg_names = ("dim", "yukawa_lambda_name") + supports_rscale = True def __init__(self, dim=None, yukawa_lambda_name="lam"): """ @@ -544,15 +549,27 @@ class YukawaKernel(ExpressionKernel): else: raise RuntimeError("unsupported dimensionality") - ExpressionKernel.__init__( - self, + super(YukawaKernel, self).__init__( dim, - expression=expr, - scaling=scaling, + proxy_expression=expr, + global_scaling_const=scaling, is_complex_valued=True) self.yukawa_lambda_name = yukawa_lambda_name + def adjust_proxy_expression(self, expr, rscale, nderivatives): + import sympy as sp + result = expr.subs( + sp.Symbol(self.yukawa_lambda_name), + sp.Symbol(self.yukawa_lambda_name) * rscale) + + if self.dim == 2: + return result + elif self.dim == 3: + return result / rscale**(nderivatives+1) + else: + raise RuntimeError("unsupported dimensionality") + def __getinitargs__(self): return (self.dim, self.yukawa_lambda_name) @@ -584,6 +601,7 @@ class YukawaKernel(ExpressionKernel): class StokesletKernel(ExpressionKernel): init_arg_names = ("dim", "icomp", "jcomp", "viscosity_mu_name") + supports_rscale = False def __init__(self, dim, icomp, jcomp, viscosity_mu_name="mu"): r""" @@ -623,13 +641,16 @@ class StokesletKernel(ExpressionKernel): self.icomp = icomp self.jcomp = jcomp - ExpressionKernel.__init__( - self, + super(StokesletKernel, self).__init__( dim, - expression=expr, - scaling=scaling, + proxy_expression=expr, + global_scaling_const=scaling, is_complex_valued=False) + def adjust_proxy_expression(self, expr, rscale, nderivatives): + assert rscale == 1 + return expr + def __getinitargs__(self): return (self.dim, self.icomp, self.jcomp, self.viscosity_mu_name) @@ -652,6 +673,7 @@ class StokesletKernel(ExpressionKernel): class StressletKernel(ExpressionKernel): init_arg_names = ("dim", "icomp", "jcomp", "kcomp", "viscosity_mu_name") + supports_rscale = False def __init__(self, dim=None, icomp=None, jcomp=None, kcomp=None, viscosity_mu_name="mu"): @@ -689,13 +711,16 @@ class StressletKernel(ExpressionKernel): self.jcomp = jcomp self.kcomp = kcomp - ExpressionKernel.__init__( - self, + super(StressletKernel, self).__init__( dim, - expression=expr, - scaling=scaling, + proxy_expression=expr, + global_scaling_const=scaling, is_complex_valued=False) + def adjust_proxy_expression(self, expr, rscale, nderivatives): + assert rscale == 1 + return expr + def __getinitargs__(self): return (self.dim, self.icomp, self.jcomp, self.kcomp, self.viscosity_mu_name)