From 10809fe9fab43e9e4cd92e2842d1ae4537701032 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 23 Nov 2020 20:00:59 -0600 Subject: [PATCH 01/32] Additions enabling taqbx --- sumpy/expansion/local.py | 13 ++++++++----- sumpy/qbx.py | 30 +++++++++++++++++++++++++----- 2 files changed, 33 insertions(+), 10 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 436707e1..41398a19 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -114,12 +114,15 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): """ def coefficients_from_source(self, kernel, avec, bvec, rscale, sac=None): - from sumpy.tools import MiDerivativeTaker - ppkernel = kernel.postprocess_at_source(kernel.get_expression(avec), avec) + from sumpy.tools import MiDerivativeTaker, my_syntactic_subs - taker = MiDerivativeTaker(ppkernel, avec) - return [ - taker.diff(mi) * rscale ** sum(mi) + rad = avec + bvec # source to center + center to target = source to target + ppkernel = kernel.postprocess_at_source(kernel.get_expression(rad), rad) + + taker = MiDerivativeTaker(ppkernel, bvec) + + return [my_syntactic_subs(taker.diff(mi) * rscale ** sum(mi), + {bvc: 0 for bvc in bvec}) for mi in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec, rscale, sac=None): diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 9e556ccc..6be0633a 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -86,17 +86,22 @@ def expand(expansion_nr, sac, expansion, avec, bvec): class LayerPotentialBase(KernelComputation, KernelCacheWrapper): def __init__(self, ctx, expansions, strength_usage=None, value_dtypes=None, - name=None, device=None): + name=None, device=None, scaling_at_target=None): KernelComputation.__init__(self, ctx, expansions, strength_usage, value_dtypes, name, device) from pytools import single_valued self.dim = single_valued(knl.dim for knl in self.expansions) + self._vector_names = {"a", "b"} + + # if not None, a sympy expression for the scaling + self._scaling_at_target = scaling_at_target def get_cache_key(self): return (type(self).__name__, tuple(self.kernels), - tuple(self.strength_usage), tuple(self.value_dtypes)) + tuple(self.strength_usage), tuple(self.value_dtypes), + str(self._scaling_at_target)) @property def expansions(self): @@ -122,7 +127,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): from sumpy.codegen import to_loopy_insns loopy_insns = to_loopy_insns( sac.assignments.items(), - vector_names={"a", "b"}, + vector_names=self._vector_names, pymbolic_expr_maps=[ expn.kernel.get_code_transformer() for expn in self.expansions], retain_names=result_names, @@ -215,6 +220,20 @@ class LayerPotential(LayerPotentialBase): None, shape="ntargets", order="C") for i in range(len(self.kernels))]) + if self._scaling_at_target: + sat = "sat * " # inject the multiplicative scaling + from sumpy.codegen import to_loopy_insns + sat_insns = to_loopy_insns( + [("sat", self._scaling_at_target)], + vector_names=self._vector_names, + pymbolic_expr_maps=[ + expn.kernel.get_code_transformer() for expn in self.expansions], + retain_names=["sat"], + complex_dtype=np.complex128) # FIXME + else: + sat = "" + sat_insns = [] + loopy_knl = lp.make_kernel([""" {[itgt, isrc, idim]: \ 0 <= itgt < ntargets and \ @@ -227,11 +246,12 @@ class LayerPotential(LayerPotentialBase): + ["<> b[idim] = tgt[idim, itgt] - center[idim, itgt] {dup=idim}"] + ["<> rscale = expansion_radii[itgt]"] + loopy_insns + kernel_exprs + + sat_insns + [""" - result_{i}[itgt] = knl_{i}_scaling * \ + result_{i}[itgt] = knl_{i}_scaling * {sat} \ simul_reduce(sum, isrc, pair_result_{i}) \ {{id_prefix=write_lpot,inames=itgt}} - """.format(i=iknl) + """.format(i=iknl, sat=sat) for iknl in range(len(self.expansions))] + ["end"], arguments, -- GitLab From 8387357955e5c97b01aa3d5d405edf219f89e755 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 23 Nov 2020 20:14:25 -0600 Subject: [PATCH 02/32] Feed avec to postprocess_at_source --- sumpy/expansion/local.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 41398a19..c21eaa16 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -117,7 +117,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): from sumpy.tools import MiDerivativeTaker, my_syntactic_subs rad = avec + bvec # source to center + center to target = source to target - ppkernel = kernel.postprocess_at_source(kernel.get_expression(rad), rad) + ppkernel = kernel.postprocess_at_source(kernel.get_expression(rad), avec) taker = MiDerivativeTaker(ppkernel, bvec) -- GitLab From ae1c8146cd7ec0ac21b469628e5fad746020a7a0 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 23 Nov 2020 21:14:26 -0600 Subject: [PATCH 03/32] Make the changes backward compatible --- sumpy/expansion/local.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index c21eaa16..ea6baa26 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -116,14 +116,26 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): def coefficients_from_source(self, kernel, avec, bvec, rscale, sac=None): from sumpy.tools import MiDerivativeTaker, my_syntactic_subs - rad = avec + bvec # source to center + center to target = source to target - ppkernel = kernel.postprocess_at_source(kernel.get_expression(rad), avec) + if bvec is not None: + # source to center + center to target = source to target + rad = avec + bvec + ppkernel = kernel.postprocess_at_source( + kernel.get_expression(rad), avec) - taker = MiDerivativeTaker(ppkernel, bvec) + taker = MiDerivativeTaker(ppkernel, bvec) - return [my_syntactic_subs(taker.diff(mi) * rscale ** sum(mi), - {bvc: 0 for bvc in bvec}) - for mi in self.get_coefficient_identifiers()] + return [my_syntactic_subs(taker.diff(mi) * rscale ** sum(mi), + {bvc: 0 for bvc in bvec}) + for mi in self.get_coefficient_identifiers()] + + else: + # bvec may be unused at all + ppkernel = kernel.postprocess_at_source( + kernel.get_expression(avec), avec) + + taker = MiDerivativeTaker(ppkernel, avec) + return [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 -- GitLab From 7e1d322ff164102931855e51149a6143cc2c42d4 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 24 Nov 2020 23:03:16 -0600 Subject: [PATCH 04/32] Move rescaled kernel to sumpy, and support multi-kernel eval --- sumpy/kernel.py | 42 ++++++++++++++++++++++++++++++++++++++++++ sumpy/p2e.py | 2 -- sumpy/qbx.py | 25 +++++++++++++++---------- 3 files changed, 57 insertions(+), 12 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 16e740a1..13e6489a 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -974,6 +974,48 @@ class DirectionalSourceDerivative(DirectionalDerivative): # }}} +# {{{ rescaled kernels used by QBMAX + +class AsymptoticallyRescaledKernel(KernelWrapper): + init_arg_names = ("inner_kernel", "rescaled_expression", "scaling_expression") + + def __init__(self, inner_kernel, rescaled_expression, scaling_expression): + """ + :param rescaled_expression: A pymbolic expression. The kernel + expression after scaling. + :param scaling_expression: A pymbolic expression. The multiplier + used for scaling the inner_kernel. + """ + super().__init__(inner_kernel) + self.rescaled_expression = rescaled_expression + self.scaling_expression = scaling_expression + + def __repr__(self): + return "AsymKnl[%s]%dD" % (self.inner_kernel.__repr__(), self.dim) + + def __getinitargs__(self): + return (self.inner_kernel, self.rescaled_expression, self.scaling_expression) + + def get_expression(self, scaled_dist_vec): + from pymbolic.interop.sympy import PymbolicToSympyMapper + expr = PymbolicToSympyMapper()(self.rescaled_expression) + return expr + + def get_scaling_expression(self, scaled_dist_vec): + from pymbolic.interop.sympy import PymbolicToSympyMapper + expr = PymbolicToSympyMapper()(self.scaling_expression) + return expr + + def update_persistent_hash(self, key_hash, key_builder): + self.base_kernel.update_persistent_hash(key_hash, key_builder) + key_hash.update(str(self.rescaled_expression).encode("utf8")) + key_hash.update(str(self.scaling_expression).encode("utf8")) + + mapper_method = "map_asymptotically_rescaled_kernel" + +# }}} + + # {{{ kernel mappers class KernelMapper: diff --git a/sumpy/p2e.py b/sumpy/p2e.py index 93692e6e..d1dce966 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -313,5 +313,3 @@ class P2EFromCSR(P2EBase): return knl(queue, centers=centers, rscale=rscale, **kwargs) # }}} - -# vim: foldmethod=marker diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 6be0633a..43311e52 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -86,7 +86,7 @@ def expand(expansion_nr, sac, expansion, avec, bvec): class LayerPotentialBase(KernelComputation, KernelCacheWrapper): def __init__(self, ctx, expansions, strength_usage=None, value_dtypes=None, - name=None, device=None, scaling_at_target=None): + name=None, device=None, scales_at_target=None): KernelComputation.__init__(self, ctx, expansions, strength_usage, value_dtypes, name, device) @@ -95,13 +95,16 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): self.dim = single_valued(knl.dim for knl in self.expansions) self._vector_names = {"a", "b"} - # if not None, a sympy expression for the scaling - self._scaling_at_target = scaling_at_target + # an optional list of sympy expressions for the scales applied + # to each expansion at last + if scales_at_target: + assert len(scales_at_target) == len(expansions) + self._scales_at_target = scales_at_target def get_cache_key(self): return (type(self).__name__, tuple(self.kernels), tuple(self.strength_usage), tuple(self.value_dtypes), - str(self._scaling_at_target)) + str(self._scales_at_target)) @property def expansions(self): @@ -220,18 +223,20 @@ class LayerPotential(LayerPotentialBase): None, shape="ntargets", order="C") for i in range(len(self.kernels))]) - if self._scaling_at_target: - sat = "sat * " # inject the multiplicative scaling + if self._scales_at_target: + # inject the multiplicative scales of each output kernel + sats = [f"sat_{iknl} * " for iknl in range(len(self.expansions))] from sumpy.codegen import to_loopy_insns sat_insns = to_loopy_insns( - [("sat", self._scaling_at_target)], + [(f"sat_{iknl}", self._scale_at_target[iknl]) + for iknl in range(len(self.expansions))], vector_names=self._vector_names, pymbolic_expr_maps=[ expn.kernel.get_code_transformer() for expn in self.expansions], - retain_names=["sat"], + retain_names=[f"sat_{iknl}" for iknl in range(len(self.expansions))], complex_dtype=np.complex128) # FIXME else: - sat = "" + sats = ["", ] * len(self.expansions) sat_insns = [] loopy_knl = lp.make_kernel([""" @@ -251,7 +256,7 @@ class LayerPotential(LayerPotentialBase): result_{i}[itgt] = knl_{i}_scaling * {sat} \ simul_reduce(sum, isrc, pair_result_{i}) \ {{id_prefix=write_lpot,inames=itgt}} - """.format(i=iknl, sat=sat) + """.format(i=iknl, sat=sats[iknl]) for iknl in range(len(self.expansions))] + ["end"], arguments, -- GitLab From 6bb4f448e0e786e7a68dfea038c149588307ebd7 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 30 Nov 2020 12:00:38 -0600 Subject: [PATCH 05/32] Handle expansions and kernels --- sumpy/kernel.py | 14 ++++++++----- sumpy/qbx.py | 54 +++++++++++++++++++++++++++++++------------------ 2 files changed, 43 insertions(+), 25 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 13e6489a..e6e015d1 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -996,20 +996,24 @@ class AsymptoticallyRescaledKernel(KernelWrapper): def __getinitargs__(self): return (self.inner_kernel, self.rescaled_expression, self.scaling_expression) - def get_expression(self, scaled_dist_vec): + def get_expression(self, scaled_dist_vec=None): from pymbolic.interop.sympy import PymbolicToSympyMapper expr = PymbolicToSympyMapper()(self.rescaled_expression) return expr - def get_scaling_expression(self, scaled_dist_vec): + def get_scaling_expression(self, scaled_dist_vec=None): from pymbolic.interop.sympy import PymbolicToSympyMapper expr = PymbolicToSympyMapper()(self.scaling_expression) return expr def update_persistent_hash(self, key_hash, key_builder): - self.base_kernel.update_persistent_hash(key_hash, key_builder) - key_hash.update(str(self.rescaled_expression).encode("utf8")) - key_hash.update(str(self.scaling_expression).encode("utf8")) + for name, value in zip(self.init_arg_names, self.__getinitargs__()): + if name in ["rescaled_expression", "scaling_expression"]: + from pymbolic.mapper.persistent_hash import ( + PersistentHashWalkMapper as PersistentHashWalkMapper) + PersistentHashWalkMapper(key_hash)(value) + else: + key_builder.rec(key_hash, value) mapper_method = "map_asymptotically_rescaled_kernel" diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 43311e52..66b85610 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -86,7 +86,7 @@ def expand(expansion_nr, sac, expansion, avec, bvec): class LayerPotentialBase(KernelComputation, KernelCacheWrapper): def __init__(self, ctx, expansions, strength_usage=None, value_dtypes=None, - name=None, device=None, scales_at_target=None): + name=None, device=None): KernelComputation.__init__(self, ctx, expansions, strength_usage, value_dtypes, name, device) @@ -95,16 +95,9 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): self.dim = single_valued(knl.dim for knl in self.expansions) self._vector_names = {"a", "b"} - # an optional list of sympy expressions for the scales applied - # to each expansion at last - if scales_at_target: - assert len(scales_at_target) == len(expansions) - self._scales_at_target = scales_at_target - def get_cache_key(self): return (type(self).__name__, tuple(self.kernels), - tuple(self.strength_usage), tuple(self.value_dtypes), - str(self._scales_at_target)) + tuple(self.strength_usage), tuple(self.value_dtypes)) @property def expansions(self): @@ -223,20 +216,41 @@ class LayerPotential(LayerPotentialBase): None, shape="ntargets", order="C") for i in range(len(self.kernels))]) - if self._scales_at_target: - # inject the multiplicative scales of each output kernel - sats = [f"sat_{iknl} * " for iknl in range(len(self.expansions))] + # "expansions" = kernels + sats = [] + sat_insns_data = [] + for iknl in range(len(self.kernels)): + + from sumpy.expansion import ExpansionBase + if isinstance(self.expansions[iknl], ExpansionBase): + knl = self.expansions[iknl].kernel + else: + from sumpy.kernel import Kernel + assert isinstance(self.expansions[iknl], Kernel) + knl = self.expansions[iknl] + + try: + sat_expr = knl.get_scaling_expression() + except AttributeError: + sat_expr = None + + if sat_expr: + # inject the multiplicative scales of each output kernel + sats.append(f"sat_{iknl} * ") + sat_insns_data.append((f"sat_{iknl}", sat_expr)) + else: + sats.append("") + + if sat_insns_data: from sumpy.codegen import to_loopy_insns sat_insns = to_loopy_insns( - [(f"sat_{iknl}", self._scale_at_target[iknl]) - for iknl in range(len(self.expansions))], - vector_names=self._vector_names, - pymbolic_expr_maps=[ - expn.kernel.get_code_transformer() for expn in self.expansions], - retain_names=[f"sat_{iknl}" for iknl in range(len(self.expansions))], - complex_dtype=np.complex128) # FIXME + sat_insns_data, + vector_names=self._vector_names, + pymbolic_expr_maps=[ + expn.kernel.get_code_transformer() for expn in self.expansions], + retain_names=[datpair[0] for datpair in sat_insns_data], + complex_dtype=np.complex128) # FIXME else: - sats = ["", ] * len(self.expansions) sat_insns = [] loopy_knl = lp.make_kernel([""" -- GitLab From 06589bcafaafee942dd98ccae769caa846aa0b0d Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 30 Nov 2020 12:01:54 -0600 Subject: [PATCH 06/32] Flake8 fix --- sumpy/qbx.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 66b85610..5edcd403 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -247,7 +247,8 @@ class LayerPotential(LayerPotentialBase): sat_insns_data, vector_names=self._vector_names, pymbolic_expr_maps=[ - expn.kernel.get_code_transformer() for expn in self.expansions], + expn.kernel.get_code_transformer() + for expn in self.expansions], retain_names=[datpair[0] for datpair in sat_insns_data], complex_dtype=np.complex128) # FIXME else: -- GitLab From ccd1a8e6d023086017145ea3857a6b8b32f1219b Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 30 Nov 2020 21:17:49 -0600 Subject: [PATCH 07/32] Codegen support for fmm --- sumpy/e2p.py | 9 ++++++--- sumpy/kernel.py | 6 ++++++ sumpy/p2e.py | 22 +++++++++++++++++++--- 3 files changed, 31 insertions(+), 6 deletions(-) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index 2512aca8..cd7c44fe 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -56,14 +56,16 @@ class E2PBase(KernelCacheWrapper): if device is None: device = ctx.devices[0] - from sumpy.kernel import SourceDerivativeRemover, TargetDerivativeRemover + from sumpy.kernel import (SourceDerivativeRemover, TargetDerivativeRemover, + AsymptoticScalingRemover) sdr = SourceDerivativeRemover() tdr = TargetDerivativeRemover() + asr = AsymptoticScalingRemover() expansion = expansion.with_kernel( sdr(expansion.kernel)) for knl in kernels: - assert sdr(tdr(knl)) == expansion.kernel + assert sdr(tdr(knl)) == asr(expansion.kernel) self.ctx = ctx self.expansion = expansion @@ -72,6 +74,7 @@ class E2PBase(KernelCacheWrapper): self.device = device self.dim = expansion.dim + self._vector_names = {"b"} def get_loopy_insns_and_result_names(self): from sumpy.symbolic import make_sym_vector @@ -98,7 +101,7 @@ class E2PBase(KernelCacheWrapper): from sumpy.codegen import to_loopy_insns loopy_insns = to_loopy_insns( sac.assignments.items(), - vector_names={"b"}, + vector_names=self._vector_names, pymbolic_expr_maps=[self.expansion.get_code_transformer()], retain_names=result_names, complex_dtype=np.complex128 # FIXME diff --git a/sumpy/kernel.py b/sumpy/kernel.py index e6e015d1..8203ff2e 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1058,6 +1058,7 @@ class KernelIdentityMapper(KernelMapper): map_yukawa_kernel = map_expression_kernel map_stokeslet_kernel = map_expression_kernel map_stresslet_kernel = map_expression_kernel + map_asymptotically_rescaled_kernel = map_expression_kernel def map_axis_target_derivative(self, kernel): return AxisTargetDerivative(kernel.axis, self.rec(kernel.inner_kernel)) @@ -1105,6 +1106,11 @@ class DerivativeCounter(KernelCombineMapper): map_directional_target_derivative = map_axis_target_derivative map_directional_source_derivative = map_axis_target_derivative + +class AsymptoticScalingRemover(KernelIdentityMapper): + def map_asymptotically_rescaled_kernel(self, kernel): + return self.rec(kernel.inner_kernel) + # }}} diff --git a/sumpy/p2e.py b/sumpy/p2e.py index d1dce966..1fcd4855 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -68,13 +68,29 @@ class P2EBase(KernelComputation, KernelCacheWrapper): expansion = expansion.with_kernel( TargetDerivativeRemover()(expansion.kernel)) + from sumpy.kernel import AsymptoticallyRescaledKernel + base_knl = TargetDerivativeRemover()(expansion.kernel) + if base_knl != expansion.kernel: + if isinstance(base_knl, AsymptoticallyRescaledKernel): + # TODO: target derivatives of scaled kernels using Leibniz's rule + # (fg)' = f'g + fg' ==> f' = [(fg)' - fg'] / g + raise ValueError( + "derivatives of scaled kernel is not currently supported.") + self.expansion = expansion self.dim = expansion.dim + self._vector_names = {"a"} - def get_loopy_instructions(self): + def get_loopy_instructions(self, isotropic_expansion=True): from sumpy.symbolic import make_sym_vector avec = make_sym_vector("a", self.dim) + if isotropic_expansion: + bvec = None + else: + # non-isotropic expansions like QBMAX depend on bvec + bvec = make_sym_vector("b", self.dim) + import sumpy.symbolic as sp rscale = sp.Symbol("rscale") @@ -84,7 +100,7 @@ class P2EBase(KernelComputation, KernelCacheWrapper): coeff_names = [] for knl_idx, kernel in enumerate(self.kernels): for i, coeff_i in enumerate( - self.expansion.coefficients_from_source(kernel, avec, None, rscale, + self.expansion.coefficients_from_source(kernel, avec, bvec, rscale, sac) ): sac.add_assignment(f"coeff{i}_{knl_idx}", coeff_i) @@ -98,7 +114,7 @@ class P2EBase(KernelComputation, KernelCacheWrapper): from sumpy.codegen import to_loopy_insns return to_loopy_insns( sac.assignments.items(), - vector_names={"a"}, + vector_names=self._vector_names, pymbolic_expr_maps=code_transformers, retain_names=coeff_names, complex_dtype=np.complex128 # FIXME -- GitLab From 46e30f9fcbe0f98b3a18446f576ae52b6391cf44 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 7 Dec 2020 17:33:59 -0600 Subject: [PATCH 08/32] Refactor to make room for kernel derivatives --- sumpy/kernel.py | 259 ++++++++++++++++++++++++++++++++++++++++++++---- sumpy/p2p.py | 3 +- sumpy/qbx.py | 2 +- 3 files changed, 241 insertions(+), 23 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 8203ff2e..6228b967 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -23,11 +23,12 @@ THE SOFTWARE. import loopy as lp import numpy as np +import sumpy.symbolic as sym from pymbolic.mapper import IdentityMapper, CSECachingMapperMixin from sumpy.symbolic import pymbolic_real_norm_2 from pymbolic.primitives import make_sym_vector from pymbolic import var - +from pytools import memoize_method __doc__ = """ Kernel interface @@ -976,39 +977,255 @@ class DirectionalSourceDerivative(DirectionalDerivative): # {{{ rescaled kernels used by QBMAX -class AsymptoticallyRescaledKernel(KernelWrapper): - init_arg_names = ("inner_kernel", "rescaled_expression", "scaling_expression") +class _AsymptoticallyRescaledKernelExpressionFactory: + """The rescaled kernels depend on expansion methods. This class isolates + mathematical asymptotic expansions that depend solely on the kernel and + abstract geometry. + + .. attribute :: dim + + The dimension of the ambient space. + + .. attribute :: kernel + + A :class:`~sumpy.kernel.Kernel` object. + + .. attribute :: expr + + The asymptotic expansions. A scalar expression that depends on distance + from the boundary (denoted by the symbol "dist"), satisfying the + normalizing condition: + + expr(dist=0) = 1 + + .. attribute :: geometric_order + + Order of geometric approximation when computing the distance to + boundary. For example, when order=1, the distance is calculated using a + first-order approximation to the boundary using QBX disk's radius and + the point of tangency. + """ + + def __init__(self, kernel, expr, geometric_order=1): + self.dim = kernel.dim + self.kernel = kernel + self.expr = expr + self.geometric_order = geometric_order + + @memoize_method + def _get_scaling_at_target(self, expn_class): + """Target-dependent scaling expression, used for scaling back the + potential after summing up the expansion. + """ + full_scaling = self._get_scaling_for_expansion(expn_class) + + from sumpy.expansion.local import ( + LineTaylorLocalExpansion, VolumeTaylorLocalExpansionBase) + if issubclass(expn_class, LineTaylorLocalExpansion): + sym_map = {sym.sym.Symbol("tau"): 1} # target is at tau=1 + elif issubclass(expn_class, VolumeTaylorLocalExpansionBase): + sym_map = {} + else: + raise ValueError("unsupported expansion class") + + return full_scaling.subs(sym_map) + + @memoize_method + def _get_scaling_for_expansion(self, expn_class): + """Target-dependent scaling expression, including the virtual variables + like "tau" for taking expansions. + """ + sym_map = self._get_symbol_map(expn_class) + expr_h = self._approximate_expr() + scaling = expr_h.subs(sym_map) + return scaling + + @memoize_method + def _get_rescaled_kernel_expr(self, expn_class): + """Apply scaling to the kernel expression. + """ + kernel_expr = self._get_kernel_expression_at_source(expn_class) + scaling = self._get_scaling_for_expansion(expn_class) + return kernel_expr / scaling + + @memoize_method + def _approximate_expr(self): + """Transform the expression by approximating the distance with + geometric information known to the QBX routine. + """ + normal, target, dist, rad = ( + self._nvec, self._tvec, self._dist, self._radius) + + if self.geometric_order == 1: + dist_approx = rad - normal.dot(target) + else: + # need curvature information for higher order approximations + raise NotImplementedError + + return self.expr.subs(dist, dist_approx) - def __init__(self, inner_kernel, rescaled_expression, scaling_expression): + @memoize_method + def _get_symbol_map(self, expn_class): + """Get interpretation of ``nvec`` and ``tvec`` that depends on the + expansion class. """ - :param rescaled_expression: A pymbolic expression. The kernel - expression after scaling. + from sumpy.expansion.local import ( + LineTaylorLocalExpansion, VolumeTaylorLocalExpansionBase) + + bvec = sym.make_sym_vector("b", self.dim) + + if issubclass(expn_class, LineTaylorLocalExpansion): + # in a line expansion, bvec is the normal, and the actual target is + # identified with tau. + + bvec_norm = sym.sym.sqrt(bvec.dot(bvec)) + sym_map = { + self._radius: bvec_norm + } + + # self._nvec: bvec / bvec_norm + for lhs, rhs in zip(self._nvec, bvec): + sym_map[lhs] = rhs / bvec_norm + + # self._tvec: sym.Symbol("tau") * bvec + for lhs, rhs in zip(self._tvec, bvec): + sym_map[lhs] = sym.sym.Symbol("tau") * rhs + + elif issubclass(expn_class, VolumeTaylorLocalExpansionBase): + sym_map = { + # disk radius is used as rscale + self._radius: sym.sym.Symbol("rscale"), + } + + # self._nvec: normal / norm(normal) + normal = sym.make_sym_vector("normal", self.dim) + normal_norm = sym.sym.sqrt(normal.dot(normal)) + for lhs, rhs in zip(self._nvec, normal): + sym_map[lhs] = rhs / normal_norm + + # self._tvec: bvec + for lhs, rhs in zip(self._tvec, bvec): + sym_map[lhs] = rhs + + else: + raise ValueError("unsupported expansion class") + + return sym_map + + @memoize_method + def _get_kernel_expression_at_source(self, expn_class): + """Kernel expression used by ``coefficients_from_source``. + """ + from sumpy.expansion.local import ( + LineTaylorLocalExpansion, VolumeTaylorLocalExpansionBase) + + avec = sym.make_sym_vector("a", self.dim) + bvec = sym.make_sym_vector("b", self.dim) + + if issubclass(expn_class, LineTaylorLocalExpansion): + tau = sym.sym.Symbol("tau") + avec_line = avec + tau*bvec + return self.kernel.get_expression(avec_line) + elif issubclass(expn_class, VolumeTaylorLocalExpansionBase): + rad = avec + bvec + return self.kernel.postprocess_at_source( + self.kernel.get_expression(rad), avec) + else: + raise ValueError( + f"unsupported expansion class for QBMAX: {expn_class.__name__}") + + @property + def _dist(self): + """The distance of the target point from the boundary. + """ + return sym.sym.Symbol("dist") + + @property + def _nvec(self): + """The unit normal vector obtained by normalizing the vector from the + QBX expansion center to the point of tangency. + """ + return sym.make_sym_vector("nvec", self.dim) + + @property + def _tvec(self): + """The vector from the QBX expansion center to the evaluation target + point. + """ + return sym.make_sym_vector("tvec", self.dim) + + @property + def _radius(self): + """Radius of the QBX disk, aka distance of the expansion center from + the boundary. + """ + return sym.sym.Symbol("radius") + + +class AsymptoticallyInformedKernel(KernelWrapper): + """An asymptotically informed kernel should behave the same as its inner kernel + if the handler does not know how to take advantage of the added information. + """ + init_arg_names = ("inner_kernel", "scaling_expression", "geometric_order") + + def __init__(self, inner_kernel, scaling_expression, geometric_order=1): + """ + :param inner_kernel: A raw PDE kernel being rescaled. :param scaling_expression: A pymbolic expression. The multiplier - used for scaling the inner_kernel. + used for scaling the inner_kernel. The scaling is a function + of "dist". + :param geometric_order: An integer for the geometric order. """ super().__init__(inner_kernel) - self.rescaled_expression = rescaled_expression self.scaling_expression = scaling_expression + self.geometric_order = geometric_order + self.asymptotics = _AsymptoticallyRescaledKernelExpressionFactory( + inner_kernel, scaling_expression, geometric_order) + + self._expn_class = None def __repr__(self): - return "AsymKnl[%s]%dD" % (self.inner_kernel.__repr__(), self.dim) + return "AsymKnl[%s]" % (self.inner_kernel.__repr__(),) def __getinitargs__(self): - return (self.inner_kernel, self.rescaled_expression, self.scaling_expression) - - def get_expression(self, scaled_dist_vec=None): - from pymbolic.interop.sympy import PymbolicToSympyMapper - expr = PymbolicToSympyMapper()(self.rescaled_expression) - return expr - - def get_scaling_expression(self, scaled_dist_vec=None): - from pymbolic.interop.sympy import PymbolicToSympyMapper - expr = PymbolicToSympyMapper()(self.scaling_expression) - return expr + return (self.inner_kernel, self.scaling_expression, + self.geometric_order) + + def set_expansion_class(self, expn_class): + self._expn_class = expn_class + + def get_expression(self, scaled_dist_vec, expn_class=None): + + # In this special case, we allow falling back to the inner kernel + # to be compatible with p2p. In the meanwhile, be careful not to + # include the scaling into the p2p result. + from sumpy.symbolic import make_sym_vector + dvec = make_sym_vector("d", self.dim) + if scaled_dist_vec == dvec: + return self.inner_kernel.get_expression(scaled_dist_vec) + + if not expn_class: + expn_class = self._expn_class + if not expn_class: + # fallback to the inner kernel + return self.inner_kernel.get_expression(scaled_dist_vec) + return self.asymptotics._get_rescaled_kernel_expr(expn_class) + + def get_scaling_expression(self, scaled_dist_vec, expn_class=None): + if scaled_dist_vec: + # the asymptotic expression has its own opinion on the dist vec + raise ValueError + if not expn_class: + expn_class = self._expn_class + if not expn_class: + raise ValueError( + "In order to generate expression, the asymptotic expansion " + "needs to know the expansion class.") + return self.asymptotics._get_scaling_at_target(expn_class) def update_persistent_hash(self, key_hash, key_builder): for name, value in zip(self.init_arg_names, self.__getinitargs__()): - if name in ["rescaled_expression", "scaling_expression"]: + if name in ["scaling_expression"]: from pymbolic.mapper.persistent_hash import ( PersistentHashWalkMapper as PersistentHashWalkMapper) PersistentHashWalkMapper(key_hash)(value) diff --git a/sumpy/p2p.py b/sumpy/p2p.py index 746a9ef7..5560b305 100644 --- a/sumpy/p2p.py +++ b/sumpy/p2p.py @@ -69,6 +69,7 @@ class P2PBase(KernelComputation, KernelCacheWrapper): from pytools import single_valued self.dim = single_valued(knl.dim for knl in self.kernels) + self._vector_names = {"d"} def get_cache_key(self): return (type(self).__name__, tuple(self.kernels), self.exclude_self, @@ -95,7 +96,7 @@ class P2PBase(KernelComputation, KernelCacheWrapper): from sumpy.codegen import to_loopy_insns loopy_insns = to_loopy_insns(sac.assignments.items(), - vector_names={"d"}, + vector_names=self._vector_names, pymbolic_expr_maps=[ knl.get_code_transformer() for knl in self.kernels], retain_names=result_names, diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 5edcd403..73fa0640 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -230,7 +230,7 @@ class LayerPotential(LayerPotentialBase): knl = self.expansions[iknl] try: - sat_expr = knl.get_scaling_expression() + sat_expr = knl.get_scaling_expression(None) except AttributeError: sat_expr = None -- GitLab From 701f4e531400d2f6cac93a52f3bf647988ce4c82 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 7 Dec 2020 17:59:43 -0600 Subject: [PATCH 09/32] Improve caching --- sumpy/kernel.py | 32 +++++++++++++++++++++++++++++--- 1 file changed, 29 insertions(+), 3 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 6228b967..3c487bd5 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1009,9 +1009,22 @@ class _AsymptoticallyRescaledKernelExpressionFactory: def __init__(self, kernel, expr, geometric_order=1): self.dim = kernel.dim self.kernel = kernel - self.expr = expr self.geometric_order = geometric_order + from sympy import Expr + if isinstance(expr, Expr): + from sumpy.symbolic import SympyToPymbolicMapper + self._expr = SympyToPymbolicMapper()(expr) + else: + from pymbolic.primitives import Expression + assert isinstance(expr, Expression) + self._expr = expr + + @property + def expr(self): + from sumpy.symbolic import PymbolicToSympyMapperWithSymbols + return PymbolicToSympyMapperWithSymbols()(self._expr) + @memoize_method def _get_scaling_at_target(self, expn_class): """Target-dependent scaling expression, used for scaling back the @@ -1171,13 +1184,22 @@ class AsymptoticallyInformedKernel(KernelWrapper): def __init__(self, inner_kernel, scaling_expression, geometric_order=1): """ :param inner_kernel: A raw PDE kernel being rescaled. - :param scaling_expression: A pymbolic expression. The multiplier + :param scaling_expression: A pymbolic/sympy expression. The multiplier used for scaling the inner_kernel. The scaling is a function of "dist". :param geometric_order: An integer for the geometric order. """ super().__init__(inner_kernel) - self.scaling_expression = scaling_expression + + from sympy import Expr + if isinstance(scaling_expression, Expr): + from sumpy.symbolic import SympyToPymbolicMapper + self.scaling_expression = SympyToPymbolicMapper()(scaling_expression) + else: + from pymbolic.primitives import Expression + assert isinstance(scaling_expression, Expression) + self.scaling_expression = scaling_expression + self.geometric_order = geometric_order self.asymptotics = _AsymptoticallyRescaledKernelExpressionFactory( inner_kernel, scaling_expression, geometric_order) @@ -1191,6 +1213,10 @@ class AsymptoticallyInformedKernel(KernelWrapper): return (self.inner_kernel, self.scaling_expression, self.geometric_order) + def replace_inner_kernel(self, new_inner_kernel): + return type(self)(new_inner_kernel, + self.scaling_expression, self.geometric_order) + def set_expansion_class(self, expn_class): self._expn_class = expn_class -- GitLab From 5402e6ad4424d84359419f497f17929897aaeed3 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 7 Dec 2020 18:50:50 -0600 Subject: [PATCH 10/32] Remove premature postprocess_at_source --- sumpy/kernel.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 3c487bd5..95636b7a 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1141,8 +1141,7 @@ class _AsymptoticallyRescaledKernelExpressionFactory: return self.kernel.get_expression(avec_line) elif issubclass(expn_class, VolumeTaylorLocalExpansionBase): rad = avec + bvec - return self.kernel.postprocess_at_source( - self.kernel.get_expression(rad), avec) + return self.kernel.get_expression(rad) else: raise ValueError( f"unsupported expansion class for QBMAX: {expn_class.__name__}") -- GitLab From 0cec98de1d77268726ccd567a6a24f996574c0a2 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 7 Dec 2020 21:06:29 -0600 Subject: [PATCH 11/32] More fixes --- sumpy/e2p.py | 11 ++++++----- sumpy/kernel.py | 43 +++++++++++++++++++++++++++++++++---------- sumpy/p2e.py | 4 ++-- 3 files changed, 41 insertions(+), 17 deletions(-) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index cd7c44fe..20e01afe 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -56,16 +56,17 @@ class E2PBase(KernelCacheWrapper): if device is None: device = ctx.devices[0] - from sumpy.kernel import (SourceDerivativeRemover, TargetDerivativeRemover, - AsymptoticScalingRemover) + from sumpy.kernel import ( + SourceDerivativeRemover, TargetDerivativeRemover, + AsymptoticScalingRemover) sdr = SourceDerivativeRemover() tdr = TargetDerivativeRemover() - asr = AsymptoticScalingRemover() + adr = AsymptoticScalingRemover() expansion = expansion.with_kernel( - sdr(expansion.kernel)) + adr(sdr(expansion.kernel))) for knl in kernels: - assert sdr(tdr(knl)) == asr(expansion.kernel) + assert adr(sdr(tdr(knl))) == expansion.kernel self.ctx = ctx self.expansion = expansion diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 95636b7a..7fbda53c 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1220,20 +1220,28 @@ class AsymptoticallyInformedKernel(KernelWrapper): self._expn_class = expn_class def get_expression(self, scaled_dist_vec, expn_class=None): - - # In this special case, we allow falling back to the inner kernel - # to be compatible with p2p. In the meanwhile, be careful not to - # include the scaling into the p2p result. + """To get rescaled expression, pass in scaled_dist_vec = *None*. + """ + # In these special cases, we allow falling back to the inner kernel + # to be compatible with p2p, p2e and e2p. + # In the meanwhile, be careful not to include the scaling into the p2p + # result. from sumpy.symbolic import make_sym_vector + avec = make_sym_vector("a", self.dim) + bvec = make_sym_vector("b", self.dim) dvec = make_sym_vector("d", self.dim) - if scaled_dist_vec == dvec: + if scaled_dist_vec in [dvec, bvec, avec]: return self.inner_kernel.get_expression(scaled_dist_vec) + # Return properly rescaled expression based on expansion class + assert scaled_dist_vec == avec + bvec if not expn_class: expn_class = self._expn_class if not expn_class: - # fallback to the inner kernel - return self.inner_kernel.get_expression(scaled_dist_vec) + raise ValueError( + "In order to generate expression, the asymptotic expansion " + "needs to know the expansion class.") + return self.asymptotics._get_rescaled_kernel_expr(expn_class) def get_scaling_expression(self, scaled_dist_vec, expn_class=None): @@ -1257,7 +1265,7 @@ class AsymptoticallyInformedKernel(KernelWrapper): else: key_builder.rec(key_hash, value) - mapper_method = "map_asymptotically_rescaled_kernel" + mapper_method = "map_asymptotically_informed_kernel" # }}} @@ -1300,7 +1308,7 @@ class KernelIdentityMapper(KernelMapper): map_yukawa_kernel = map_expression_kernel map_stokeslet_kernel = map_expression_kernel map_stresslet_kernel = map_expression_kernel - map_asymptotically_rescaled_kernel = map_expression_kernel + map_asymptotically_informed_kernel = map_expression_kernel def map_axis_target_derivative(self, kernel): return AxisTargetDerivative(kernel.axis, self.rec(kernel.inner_kernel)) @@ -1317,6 +1325,13 @@ class AxisTargetDerivativeRemover(KernelIdentityMapper): def map_axis_target_derivative(self, kernel): return self.rec(kernel.inner_kernel) + def map_asymptotically_informed_kernel(self, kernel): + mapped = type(kernel)( + self.rec(kernel.inner_kernel), + kernel.scaling_expression, kernel.geometric_order) + mapped.set_expansion_class(kernel._expn_class) + return mapped + class TargetDerivativeRemover(AxisTargetDerivativeRemover): def map_directional_target_derivative(self, kernel): @@ -1327,6 +1342,13 @@ class SourceDerivativeRemover(KernelIdentityMapper): def map_directional_source_derivative(self, kernel): return self.rec(kernel.inner_kernel) + def map_asymptotically_informed_kernel(self, kernel): + mapped = type(kernel)( + self.rec(kernel.inner_kernel), + kernel.scaling_expression, kernel.geometric_order) + mapped.set_expansion_class(kernel._expn_class) + return mapped + class DerivativeCounter(KernelCombineMapper): def combine(self, values): @@ -1341,6 +1363,7 @@ class DerivativeCounter(KernelCombineMapper): map_yukawa_kernel = map_expression_kernel map_stokeslet_kernel = map_expression_kernel map_stresslet_kernel = map_expression_kernel + map_asymptotically_informed_kernel = map_expression_kernel def map_axis_target_derivative(self, kernel): return 1 + self.rec(kernel.inner_kernel) @@ -1350,7 +1373,7 @@ class DerivativeCounter(KernelCombineMapper): class AsymptoticScalingRemover(KernelIdentityMapper): - def map_asymptotically_rescaled_kernel(self, kernel): + def map_asymptotically_informed_kernel(self, kernel): return self.rec(kernel.inner_kernel) # }}} diff --git a/sumpy/p2e.py b/sumpy/p2e.py index 1fcd4855..27621759 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -68,10 +68,10 @@ class P2EBase(KernelComputation, KernelCacheWrapper): expansion = expansion.with_kernel( TargetDerivativeRemover()(expansion.kernel)) - from sumpy.kernel import AsymptoticallyRescaledKernel + from sumpy.kernel import AsymptoticallyInformedKernel base_knl = TargetDerivativeRemover()(expansion.kernel) if base_knl != expansion.kernel: - if isinstance(base_knl, AsymptoticallyRescaledKernel): + if isinstance(base_knl, AsymptoticallyInformedKernel): # TODO: target derivatives of scaled kernels using Leibniz's rule # (fg)' = f'g + fg' ==> f' = [(fg)' - fg'] / g raise ValueError( -- GitLab From c74cb73c00731fc825c4158e9c47888f910c8703 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 7 Dec 2020 21:41:59 -0600 Subject: [PATCH 12/32] Fix multipole.coefficients_from_source --- sumpy/expansion/multipole.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index e316bdfc..6247315c 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -54,7 +54,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): """ def coefficients_from_source(self, kernel, avec, bvec, rscale, sac=None): - from sumpy.kernel import DirectionalSourceDerivative + from sumpy.kernel import ( + DirectionalSourceDerivative, AsymptoticallyInformedKernel) if kernel is None: kernel = self.kernel @@ -63,11 +64,15 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): if not self.use_rscale: rscale = 1 - if isinstance(kernel, DirectionalSourceDerivative): + if isinstance(kernel, AsymptoticallyInformedKernel): + tmp_kernel = kernel.inner_kernel + else: + tmp_kernel = kernel + + if isinstance(tmp_kernel, DirectionalSourceDerivative): from sumpy.symbolic import make_sym_vector dir_vecs = [] - tmp_kernel = kernel while isinstance(tmp_kernel, DirectionalSourceDerivative): dir_vecs.append(make_sym_vector(tmp_kernel.dir_vec_name, kernel.dim)) tmp_kernel = tmp_kernel.inner_kernel -- GitLab From f4605c9736ac4023e4eda2401fb874d2774bb5ee Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 7 Dec 2020 22:27:52 -0600 Subject: [PATCH 13/32] Remove unreasonable assertion --- sumpy/kernel.py | 1 - sumpy/p2e.py | 7 ++++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 7fbda53c..b21a1236 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1234,7 +1234,6 @@ class AsymptoticallyInformedKernel(KernelWrapper): return self.inner_kernel.get_expression(scaled_dist_vec) # Return properly rescaled expression based on expansion class - assert scaled_dist_vec == avec + bvec if not expn_class: expn_class = self._expn_class if not expn_class: diff --git a/sumpy/p2e.py b/sumpy/p2e.py index 27621759..54086cda 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -73,7 +73,12 @@ class P2EBase(KernelComputation, KernelCacheWrapper): if base_knl != expansion.kernel: if isinstance(base_knl, AsymptoticallyInformedKernel): # TODO: target derivatives of scaled kernels using Leibniz's rule - # (fg)' = f'g + fg' ==> f' = [(fg)' - fg'] / g + # (fg)' = f'g + fg' ==> f' = [(fg)' - fg'] / g = (fg)'/g - fg'/g, + # i.e., + # f' ~ QBX[(fg)']/g - QBX[fg]*g'/g^2, + # where + # QBX[(fg)'] ~ QBX[fg]', + # can be computed by differentiating QBMAX local expansions. raise ValueError( "derivatives of scaled kernel is not currently supported.") -- GitLab From a2ea5e81ad105a1de21668fe350dadbdf7820a83 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 14 Dec 2020 14:30:30 -0600 Subject: [PATCH 14/32] Refactor --- sumpy/kernel.py | 45 +++++++++++++++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 14 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index b21a1236..723394d3 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1178,15 +1178,18 @@ class AsymptoticallyInformedKernel(KernelWrapper): """An asymptotically informed kernel should behave the same as its inner kernel if the handler does not know how to take advantage of the added information. """ - init_arg_names = ("inner_kernel", "scaling_expression", "geometric_order") + init_arg_names = ("inner_kernel", "scaling_expression", "geometric_order", + "expansion_class") - def __init__(self, inner_kernel, scaling_expression, geometric_order=1): + def __init__(self, inner_kernel, scaling_expression, geometric_order=1, + expansion_class=None): """ :param inner_kernel: A raw PDE kernel being rescaled. :param scaling_expression: A pymbolic/sympy expression. The multiplier used for scaling the inner_kernel. The scaling is a function of "dist". :param geometric_order: An integer for the geometric order. + :param expansion_class: Expansion class to be used with the kernel. """ super().__init__(inner_kernel) @@ -1203,21 +1206,25 @@ class AsymptoticallyInformedKernel(KernelWrapper): self.asymptotics = _AsymptoticallyRescaledKernelExpressionFactory( inner_kernel, scaling_expression, geometric_order) - self._expn_class = None + self.expansion_class = expansion_class def __repr__(self): - return "AsymKnl[%s]" % (self.inner_kernel.__repr__(),) + return "AsymKnl[%s,%s]" % (self.inner_kernel.__repr__(), + repr(self.expansion_class)) def __getinitargs__(self): return (self.inner_kernel, self.scaling_expression, - self.geometric_order) + self.geometric_order, self.expansion_class) def replace_inner_kernel(self, new_inner_kernel): return type(self)(new_inner_kernel, - self.scaling_expression, self.geometric_order) + self.scaling_expression, self.geometric_order, + self.expansion_class) - def set_expansion_class(self, expn_class): - self._expn_class = expn_class + def replace_expansion_class(self, expn_class): + return type(self)(self.inner_kernel, + self.scaling_expression, self.geometric_order, + expn_class) def get_expression(self, scaled_dist_vec, expn_class=None): """To get rescaled expression, pass in scaled_dist_vec = *None*. @@ -1235,7 +1242,7 @@ class AsymptoticallyInformedKernel(KernelWrapper): # Return properly rescaled expression based on expansion class if not expn_class: - expn_class = self._expn_class + expn_class = self.expansion_class if not expn_class: raise ValueError( "In order to generate expression, the asymptotic expansion " @@ -1248,7 +1255,7 @@ class AsymptoticallyInformedKernel(KernelWrapper): # the asymptotic expression has its own opinion on the dist vec raise ValueError if not expn_class: - expn_class = self._expn_class + expn_class = self.expansion_class if not expn_class: raise ValueError( "In order to generate expression, the asymptotic expansion " @@ -1327,8 +1334,8 @@ class AxisTargetDerivativeRemover(KernelIdentityMapper): def map_asymptotically_informed_kernel(self, kernel): mapped = type(kernel)( self.rec(kernel.inner_kernel), - kernel.scaling_expression, kernel.geometric_order) - mapped.set_expansion_class(kernel._expn_class) + kernel.scaling_expression, kernel.geometric_order, + kernel.expansion_class) return mapped @@ -1344,8 +1351,8 @@ class SourceDerivativeRemover(KernelIdentityMapper): def map_asymptotically_informed_kernel(self, kernel): mapped = type(kernel)( self.rec(kernel.inner_kernel), - kernel.scaling_expression, kernel.geometric_order) - mapped.set_expansion_class(kernel._expn_class) + kernel.scaling_expression, kernel.geometric_order, + kernel.expansion_class) return mapped @@ -1375,6 +1382,16 @@ class AsymptoticScalingRemover(KernelIdentityMapper): def map_asymptotically_informed_kernel(self, kernel): return self.rec(kernel.inner_kernel) + +class ExpaniosnClassReplacer(KernelIdentityMapper): + def __init__(self, expansion_class): + self.expansion_class = expansion_class + + def map_asymptotically_informed_kernel(self, kernel): + mapped = kernel.replace_expansion_class(self.expansion_class) + mapped_inner = self.rec(kernel.inner_kernel) + return mapped.replace_inner_kernel(mapped_inner) + # }}} -- GitLab From 0ea360e6e01199fe87ea64b0937b4b3cfa760846 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 14 Dec 2020 20:31:15 -0600 Subject: [PATCH 15/32] Fix kernel caching --- sumpy/kernel.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 723394d3..311562da 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1268,6 +1268,8 @@ class AsymptoticallyInformedKernel(KernelWrapper): from pymbolic.mapper.persistent_hash import ( PersistentHashWalkMapper as PersistentHashWalkMapper) PersistentHashWalkMapper(key_hash)(value) + elif name in ["expansion_class"]: + key_builder.rec(key_hash, repr(value)) else: key_builder.rec(key_hash, value) -- GitLab From 95133c13bde750833980a4630952943e732604e5 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 5 Jul 2021 12:12:27 -0500 Subject: [PATCH 16/32] Preliminary fixes --- sumpy/qbx.py | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/sumpy/qbx.py b/sumpy/qbx.py index f85347df..df61a764 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -240,24 +240,13 @@ class LayerPotential(LayerPotentialBase): None, shape="ntargets", order="C") for i in range(len(self.target_kernels))]) - # "expansions" = kernels sats = [] sat_insns_data = [] - for iknl in range(len(self.kernels)): - - from sumpy.expansion import ExpansionBase - if isinstance(self.expansions[iknl], ExpansionBase): - knl = self.expansions[iknl].kernel - else: - from sumpy.kernel import Kernel - assert isinstance(self.expansions[iknl], Kernel) - knl = self.expansions[iknl] - + for iknl, knl in enumerate(self.target_kernels): try: sat_expr = knl.get_scaling_expression(None) except AttributeError: sat_expr = None - if sat_expr: # inject the multiplicative scales of each output kernel sats.append(f"sat_{iknl} * ") @@ -272,7 +261,7 @@ class LayerPotential(LayerPotentialBase): vector_names=self._vector_names, pymbolic_expr_maps=[ expn.kernel.get_code_transformer() - for expn in self.expansions], + for expn in self.target_kernels], retain_names=[datpair[0] for datpair in sat_insns_data], complex_dtype=np.complex128) # FIXME else: -- GitLab From da39e83d3052126281970a47b7557cba571960b4 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 5 Jul 2021 23:30:14 -0500 Subject: [PATCH 17/32] QBMAX support --- sumpy/kernel.py | 17 +++++++++++++++-- sumpy/p2e.py | 36 ++++++++++++++++++++---------------- sumpy/qbx.py | 4 +--- sumpy/tools.py | 2 +- 4 files changed, 37 insertions(+), 22 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 71ba2be3..2954d891 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1467,6 +1467,9 @@ class _AsymptoticallyRescaledKernelExpressionFactory: class AsymptoticallyInformedKernel(KernelWrapper): """An asymptotically informed kernel should behave the same as its inner kernel if the handler does not know how to take advantage of the added information. + Specifically, the kernel keeps a record of what type of expansion will be used + on it. This information needs to be kept updated by the FMM driver for correct + scaling. """ init_arg_names = ("inner_kernel", "scaling_expression", "geometric_order", "expansion_class") @@ -1516,11 +1519,21 @@ class AsymptoticallyInformedKernel(KernelWrapper): self.scaling_expression, self.geometric_order, expn_class) + def get_base_kernel(self): + """The base kernel is wrapped with the same asymptotic information. + """ + return type(self)( + self.inner_kernel.get_base_kernel(), self.scaling_expression, + self.geometric_order, self.expansion_class) + def get_expression(self, scaled_dist_vec, expn_class=None): - """To get rescaled expression, pass in scaled_dist_vec = *None*. + """To get the asymptotically rescaled expression, pass in + scaled_dist_vec = *None*. Otherwise, it delegates the request to its + inner kernel. """ # In these special cases, we allow falling back to the inner kernel - # to be compatible with p2p, p2e and e2p. + # to be compatible with p2p, p2e and e2p. It is important to let + # other forms like a+tau*b be rescaled. # In the meanwhile, be careful not to include the scaling into the p2p # result. from sumpy.symbolic import make_sym_vector diff --git a/sumpy/p2e.py b/sumpy/p2e.py index ff3f144a..18f35d4a 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -56,8 +56,9 @@ class P2EBase(KernelComputation, KernelCacheWrapper): number of strength arrays that need to be passed. Default: all kernels use the same strength. """ - from sumpy.kernel import (TargetTransformationRemover, - SourceTransformationRemover) + from sumpy.kernel import ( + TargetTransformationRemover, + SourceTransformationRemover, AsymptoticallyInformedKernel) txr = TargetTransformationRemover() sxr = SourceTransformationRemover() @@ -70,26 +71,28 @@ class P2EBase(KernelComputation, KernelCacheWrapper): for knl in kernels: assert txr(knl) == knl - assert sxr(knl) == expansion.kernel + if sxr(knl) != expansion.kernel: + assert isinstance(sxr(knl), AsymptoticallyInformedKernel) + assert isinstance(expansion.kernel, AsymptoticallyInformedKernel) + assert sxr(knl).inner_kernel == expansion.kernel.inner_kernel + + kernels = [knl.replace_expansion_class(type(expansion)) + if isinstance(knl, AsymptoticallyInformedKernel) + else knl + for knl in kernels] KernelComputation.__init__(self, ctx=ctx, target_kernels=[], source_kernels=kernels, strength_usage=strength_usage, value_dtypes=None, name=name, device=device) - from sumpy.kernel import AsymptoticallyInformedKernel, TargetDerivativeRemover - base_knl = TargetDerivativeRemover()(expansion.kernel) - if base_knl != expansion.kernel: - if isinstance(base_knl, AsymptoticallyInformedKernel): - # TODO: target derivatives of scaled kernels using Leibniz's rule - # (fg)' = f'g + fg' ==> f' = [(fg)' - fg'] / g = (fg)'/g - fg'/g, - # i.e., - # f' ~ QBX[(fg)']/g - QBX[fg]*g'/g^2, - # where - # QBX[(fg)'] ~ QBX[fg]', - # can be computed by differentiating QBMAX local expansions. - raise ValueError( - "derivatives of scaled kernel is not currently supported.") + # TODO: target derivatives of scaled kernels using Leibniz's rule + # (fg)' = f'g + fg' ==> f' = [(fg)' - fg'] / g = (fg)'/g - fg'/g, + # i.e., + # f' ~ QBX[(fg)']/g - QBX[fg]*g'/g^2, + # where + # QBX[(fg)'] ~ QBX[fg]', + # can be computed by differentiating QBMAX local expansions. self.expansion = expansion self.dim = expansion.dim @@ -112,6 +115,7 @@ class P2EBase(KernelComputation, KernelCacheWrapper): sac = SymbolicAssignmentCollection() strengths = [sp.Symbol(f"strength_{i}") for i in self.strength_usage] + coeffs = self.expansion.coefficients_from_source_vec(self.source_kernels, avec, bvec, rscale, strengths, sac=sac) diff --git a/sumpy/qbx.py b/sumpy/qbx.py index df61a764..a805e086 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -244,6 +244,7 @@ class LayerPotential(LayerPotentialBase): sat_insns_data = [] for iknl, knl in enumerate(self.target_kernels): try: + knl = knl.replace_expansion_class(type(self.expansion)) sat_expr = knl.get_scaling_expression(None) except AttributeError: sat_expr = None @@ -259,9 +260,6 @@ class LayerPotential(LayerPotentialBase): sat_insns = to_loopy_insns( sat_insns_data, vector_names=self._vector_names, - pymbolic_expr_maps=[ - expn.kernel.get_code_transformer() - for expn in self.target_kernels], retain_names=[datpair[0] for datpair in sat_insns_data], complex_dtype=np.complex128) # FIXME else: diff --git a/sumpy/tools.py b/sumpy/tools.py index c4bb40d5..58564d2c 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -426,7 +426,7 @@ class DifferentiatedExprDerivativeTaker: for extra_mi, coeff in self.derivative_transformation.items()) if subs: - result = result.sub(subs) + result = result.subs(subs) return result * save_intermediate(1 / self.taker.rscale ** max_order) -- GitLab From 36ccedc4af1ab0addeb74c71bcc18a591e266ff5 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 6 Jul 2021 09:16:39 -0500 Subject: [PATCH 18/32] Flake8 fixes --- sumpy/expansion/multipole.py | 3 +-- sumpy/kernel.py | 1 - 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index fff3524d..408b90a9 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -57,8 +57,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): compresses them. This is more efficient that calculating full coefficients, compressing and then summing. """ - from sumpy.kernel import (KernelWrapper, - DirectionalSourceDerivative, AsymptoticallyInformedKernel) + from sumpy.kernel import KernelWrapper from sumpy.tools import mi_power, mi_factorial if not self.use_rscale: diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 2954d891..44a3ddda 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -26,7 +26,6 @@ import numpy as np import sumpy.symbolic as sym 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, parse from pytools import memoize_method -- GitLab From c29c1b7986e993ca820ac1be7f55b14948a1ca82 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 9 Jul 2021 09:41:03 -0500 Subject: [PATCH 19/32] Refactor mapper --- sumpy/kernel.py | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 44a3ddda..0afd801c 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1636,6 +1636,13 @@ class KernelIdentityMapper(KernelMapper): map_directional_source_derivative = map_directional_target_derivative + def map_asymptotically_informed_kernel(self, kernel): + mapped = type(kernel)( + self.rec(kernel.inner_kernel), + kernel.scaling_expression, kernel.geometric_order, + kernel.expansion_class) + return mapped + class AxisSourceDerivativeRemover(KernelIdentityMapper): def map_axis_source_derivative(self, kernel): @@ -1646,13 +1653,6 @@ class AxisTargetDerivativeRemover(KernelIdentityMapper): def map_axis_target_derivative(self, kernel): return self.rec(kernel.inner_kernel) - def map_asymptotically_informed_kernel(self, kernel): - mapped = type(kernel)( - self.rec(kernel.inner_kernel), - kernel.scaling_expression, kernel.geometric_order, - kernel.expansion_class) - return mapped - class TargetDerivativeRemover(AxisTargetDerivativeRemover): def map_directional_target_derivative(self, kernel): @@ -1663,13 +1663,6 @@ class SourceDerivativeRemover(AxisSourceDerivativeRemover): def map_directional_source_derivative(self, kernel): return self.rec(kernel.inner_kernel) - def map_asymptotically_informed_kernel(self, kernel): - mapped = type(kernel)( - self.rec(kernel.inner_kernel), - kernel.scaling_expression, kernel.geometric_order, - kernel.expansion_class) - return mapped - class TargetTransformationRemover(TargetDerivativeRemover): def map_target_point_multiplier(self, kernel): -- GitLab From 3c9757a7c35a8a0ba2a6b837cfb12e18c742aaa0 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 9 Jul 2021 14:48:52 -0500 Subject: [PATCH 20/32] Fix D with volume-direct --- sumpy/expansion/local.py | 79 ++++++++++++++++++++++++++-------------- sumpy/kernel.py | 4 +- 2 files changed, 54 insertions(+), 29 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index b6d74e39..095726f7 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -178,45 +178,70 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): if not self.use_rscale: rscale = 1 base_kernel = single_valued(knl.get_base_kernel() for knl in kernels) + # TODO: use a kernel mapper to find asymptotic info wrapper that is not the outmost layer + from sumpy.kernel import AsymptoticallyInformedKernel + using_qbmax = isinstance(base_kernel, AsymptoticallyInformedKernel) + + if using_qbmax and bvec is not None: + # bvec is None when, e.g., + # 1) using FMM without QBMAX + # 2) using FMM with QBMAX, during P2L + # + # With a directional expansion (like line taylor), we write: + # source to target = ( + # source to center (avec) + center to target (bvec)). + # Kernel evaluated wrt rad, postprocessed wrt avec, + # then take derivatives along the direction of bvec, and + # finally let bvec=0 + rad = avec + bvec + post_diff_subs = {bvc: 0 for bvc in bvec} + else: + base_taker = base_kernel.get_derivative_taker(avec, rscale, sac) + post_diff_subs = dict() + result = [0]*len(self) for knl, weight in zip(kernels, weights): - if bvec is not None: - # for qbmax - # source to center + center to target = source to target - # Kernel evaluated wrt rad, postprocessed wrt avec, - # then take derivatives along the direction of bvec, and - # finally let bvec=0 + # apply source derivatives, if any + if using_qbmax and bvec is not None: + # Since DifferentiatedExprDerivativeTaker does not work with + # anisotropic expansions, using vanilla expression here is + # necessary. That is, the base_taker is directly constructed + # from \partial_{n(y)}G(x-y) from sumpy.tools import ExprDerivativeTaker - rad = avec + bvec - base_taker = ExprDerivativeTaker( - base_kernel.get_expression(rad), bvec, rscale, sac) + ppkernel = knl.postprocess_at_source(knl.get_expression(rad), avec) + base_taker = ExprDerivativeTaker(ppkernel, bvec, rscale, sac) + + from pytools.tag import tag_dataclass + + @tag_dataclass + class AnisotropicDifferentiatedExprDerivativeTaker: + taker: ExprDerivativeTaker + + def diff(self, mi, save_intermediate=lambda x: x, subs=None): + result = self.taker.diff(mi) + if subs: + result = result.subs(subs) + return result * save_intermediate(1) + + taker = AnisotropicDifferentiatedExprDerivativeTaker(base_taker) + + else: taker = knl.postprocess_at_source(base_taker, avec) + # Following is a hack to make sure cse works. + if 1: def save_temp(x): return add_to_sac(sac, weight * x) for i, mi in enumerate(self.get_coefficient_identifiers()): - result[i] += taker.diff( - mi, save_temp, subs={bvc: 0 for bvc in bvec}) + result[i] += taker.diff(mi, save_temp, subs=post_diff_subs) else: - # bvec may be unused at all (None), default to usual qbx - # Kernel evaluated and post-processed wrt avec - base_taker = base_kernel.get_derivative_taker(avec, rscale, sac) - taker = knl.postprocess_at_source(base_taker, avec) - # Following is a hack to make sure cse works. - if 1: - def save_temp(x): - return add_to_sac(sac, weight * x) - - for i, mi in enumerate(self.get_coefficient_identifiers()): - result[i] += taker.diff(mi, save_temp) - else: - def save_temp(x): - return add_to_sac(sac, x) + def save_temp(x): + return add_to_sac(sac, x) - for i, mi in enumerate(self.get_coefficient_identifiers()): - result[i] += weight * taker.diff(mi, save_temp) + for i, mi in enumerate(self.get_coefficient_identifiers()): + result[i] += weight * taker.diff(mi, save_temp, subs=post_diff_subs) return result def coefficients_from_source(self, kernel, avec, bvec, rscale, sac=None): diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 0afd801c..e6343d91 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1306,7 +1306,7 @@ class _AsymptoticallyRescaledKernelExpressionFactory: self._expr = SympyToPymbolicMapper()(expr) else: from pymbolic.primitives import Expression - assert isinstance(expr, Expression) + assert isinstance(expr, (Expression, int)) self._expr = expr @property @@ -1491,7 +1491,7 @@ class AsymptoticallyInformedKernel(KernelWrapper): self.scaling_expression = SympyToPymbolicMapper()(scaling_expression) else: from pymbolic.primitives import Expression - assert isinstance(scaling_expression, Expression) + assert isinstance(scaling_expression, (Expression, int)) self.scaling_expression = scaling_expression self.geometric_order = geometric_order -- GitLab From 47d37374847ae850546b379c60ded6b279f8a553 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 9 Jul 2021 18:46:29 -0500 Subject: [PATCH 21/32] More fixes --- sumpy/expansion/local.py | 13 +++++++++++-- sumpy/kernel.py | 20 ++++++++++++++++++++ 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 095726f7..464afecb 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -178,7 +178,8 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): if not self.use_rscale: rscale = 1 base_kernel = single_valued(knl.get_base_kernel() for knl in kernels) - # TODO: use a kernel mapper to find asymptotic info wrapper that is not the outmost layer + # TODO: use a kernel mapper to find asymptotic info wrapper that is not + # the outmost layer from sumpy.kernel import AsymptoticallyInformedKernel using_qbmax = isinstance(base_kernel, AsymptoticallyInformedKernel) @@ -241,7 +242,8 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): return add_to_sac(sac, x) for i, mi in enumerate(self.get_coefficient_identifiers()): - result[i] += weight * taker.diff(mi, save_temp, subs=post_diff_subs) + result[i] += weight * taker.diff( + mi, save_temp, subs=post_diff_subs) return result def coefficients_from_source(self, kernel, avec, bvec, rscale, sac=None): @@ -259,6 +261,13 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): bvec_scaled = [b*rscale**-1 for b in bvec] from sumpy.tools import mi_power, mi_factorial + from sumpy.kernel import AsymptoticScalingDetector + asd = AsymptoticScalingDetector() + using_qbmax = (asd(kernel).get_base_kernel().dim == 0) + if using_qbmax: + # FIXME: do target derivative for QBMAX here + pass + result = sum( coeff * mi_power(bvec_scaled, mi, evaluate=False) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index e6343d91..d12ad8c8 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1564,6 +1564,13 @@ class AsymptoticallyInformedKernel(KernelWrapper): "needs to know the expansion class.") return self.asymptotics._get_scaling_at_target(expn_class) + def get_pde_as_diff_op(self): + r"""Nonlinear scaling breaks the PDE. + TODO: compression is still possible for the bulk FMM, by using full + Taylor expansion for only qbmax locals. + """ + raise NotImplementedError + def update_persistent_hash(self, key_hash, key_builder): for name, value in zip(self.init_arg_names, self.__getinitargs__()): if name in ["scaling_expression"]: @@ -1700,6 +1707,19 @@ class AsymptoticScalingRemover(KernelIdentityMapper): return self.rec(kernel.inner_kernel) +class AsymptoticScalingDetector(KernelIdentityMapper): + """Swaps base kernel with a placeholder to signal existence of asymptotic + information. This is needed because wrappers are typically ordered from + inside out as: + - base kernel + - source derivatives + - asymptotic information + - target derivatives + """ + def map_asymptotically_informed_kernel(self, kernel): + return ExpressionKernel(0, None, None, None) + + class ExpaniosnClassReplacer(KernelIdentityMapper): def __init__(self, expansion_class): self.expansion_class = expansion_class -- GitLab From 2a451e282fe1e714171dee3a3992bb7b97343c10 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 9 Jul 2021 22:04:46 -0500 Subject: [PATCH 22/32] Add a qbmax test: direct eval, line taylor --- sumpy/qbx.py | 5 +++++ test/test_qbx.py | 21 +++++++++++++-------- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/sumpy/qbx.py b/sumpy/qbx.py index a805e086..3cb8fde6 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -71,6 +71,11 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): value_dtypes=None, name=None, device=None, source_kernels=None, target_kernels=None): + from sumpy.kernel import ExpaniosnClassReplacer + ecr = ExpaniosnClassReplacer(type(expansion)) + source_kernels = tuple(ecr(sk) for sk in source_kernels) + target_kernels = tuple(ecr(tk) for tk in target_kernels) + from pytools import single_valued KernelComputation.__init__(self, ctx=ctx, target_kernels=target_kernels, strength_usage=strength_usage, source_kernels=source_kernels, diff --git a/test/test_qbx.py b/test/test_qbx.py index 379ef0e2..0eedda27 100644 --- a/test/test_qbx.py +++ b/test/test_qbx.py @@ -33,12 +33,20 @@ from sumpy.expansion.local import ( LineTaylorLocalExpansion, VolumeTaylorLocalExpansion) import pytest +from sumpy.kernel import LaplaceKernel, AsymptoticallyInformedKernel -@pytest.mark.parametrize("expn_class", [ - LineTaylorLocalExpansion, - VolumeTaylorLocalExpansion, +import sympy as sp +dist = sp.Symbol("dist") +asexpr = sp.exp(-dist) # no practical use, just to test the machinary + + +@pytest.mark.parametrize("expn_class,lknl", [ + (LineTaylorLocalExpansion, LaplaceKernel(2)), + (VolumeTaylorLocalExpansion, LaplaceKernel(2)), + (LineTaylorLocalExpansion, AsymptoticallyInformedKernel( + LaplaceKernel(2), asexpr)), ]) -def test_direct_qbx_vs_eigval(ctx_factory, expn_class): +def test_direct_qbx_vs_eigval(ctx_factory, expn_class, lknl): """This evaluates a single layer potential on a circle using a known eigenvalue/eigenvector combination. """ @@ -48,15 +56,12 @@ def test_direct_qbx_vs_eigval(ctx_factory, expn_class): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - from sumpy.kernel import LaplaceKernel - lknl = LaplaceKernel(2) - order = 12 from sumpy.qbx import LayerPotential lpot = LayerPotential(ctx, expansion=expn_class(lknl, order), - target_kernels=(lknl,), source_kernels=(lknl,)) + target_kernels=(lknl,), source_kernels=(lknl,)) mode_nr = 25 -- GitLab From ba3f60f1bdd2975165ad5b7c392c8a09c7f90c11 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 9 Jul 2021 22:43:51 -0500 Subject: [PATCH 23/32] Leave space for target derivatives --- sumpy/expansion/local.py | 38 ++++++++++++++++++++++++++++++-------- 1 file changed, 30 insertions(+), 8 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 464afecb..b9e75451 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -261,21 +261,43 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): bvec_scaled = [b*rscale**-1 for b in bvec] from sumpy.tools import mi_power, mi_factorial - from sumpy.kernel import AsymptoticScalingDetector + from sumpy.kernel import AsymptoticScalingDetector, DerivativeCounter asd = AsymptoticScalingDetector() + dcr = DerivativeCounter() using_qbmax = (asd(kernel).get_base_kernel().dim == 0) + nder = dcr(kernel) if using_qbmax: # FIXME: do target derivative for QBMAX here + # NOTE: only first order target derivative is implemented + if nder > 1: + raise NotImplementedError + # + # Potential = (int[f(x)*G(x,y)dx] * g(y) * g^{-1}(y))' + # = Expand{int[f(x)*G(x,y)dx] * g(y)}' / g(y) + # - Expand{int[f(x)*G(x,y)dx] * g(y)} * g'(y)/g^2(y) + #:= sum A_n(y) / g(y) - sum B_n(y) * g'(y)/g^2(y) + # = sum[A_n(y)/g(y) - B_n(y)*g'(y)/g^2(y)] + # = sum[A_n(y) - B_n(y)*g'(y)/g(y)] / g(y) + # ^: A_n are coeffs of usual target derivative + # ^: B_n are original coeffs + + # get the map from B_n to A_n + + # compute g'(y)/g(y) + + # assemble result pass + else: + result = sum( + coeff + * mi_power(bvec_scaled, mi, evaluate=False) + / mi_factorial(mi) + for coeff, mi in zip( + evaluated_coeffs, self.get_full_coefficient_identifiers())) + result = kernel.postprocess_at_target(result, bvec) - result = sum( - coeff - * mi_power(bvec_scaled, mi, evaluate=False) - / mi_factorial(mi) - for coeff, mi in zip( - evaluated_coeffs, self.get_full_coefficient_identifiers())) + return result - return kernel.postprocess_at_target(result, bvec) def m2l_global_precompute_nexpr(self, src_expansion): """Returns number of expressions in M2L global precomputation step. -- GitLab From f8198b130f677dd7fa363ac11fe81d8364b158aa Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 9 Jul 2021 22:44:36 -0500 Subject: [PATCH 24/32] Flake8 fix --- sumpy/expansion/local.py | 1 - 1 file changed, 1 deletion(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index b9e75451..c3ababe3 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -298,7 +298,6 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): return result - def m2l_global_precompute_nexpr(self, src_expansion): """Returns number of expressions in M2L global precomputation step. """ -- GitLab From 1b7d3b994ef6a2782bf5cdae65ffca01272394f3 Mon Sep 17 00:00:00 2001 From: xywei Date: Fri, 9 Jul 2021 23:21:22 -0500 Subject: [PATCH 25/32] Add elasticity kernel support to derivative counter --- sumpy/kernel.py | 1 + 1 file changed, 1 insertion(+) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index d12ad8c8..7d525d67 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1691,6 +1691,7 @@ class DerivativeCounter(KernelCombineMapper): map_helmholtz_kernel = map_expression_kernel map_yukawa_kernel = map_expression_kernel map_line_of_compression_kernel = map_expression_kernel + map_elasticity_kernel = map_expression_kernel map_stresslet_kernel = map_expression_kernel map_asymptotically_informed_kernel = map_expression_kernel -- GitLab From 143221df67c6a2d16aabe4c3da87ddbc73b4c1d0 Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 10 Jul 2021 12:07:42 -0500 Subject: [PATCH 26/32] Add QBMAX support for S' --- sumpy/expansion/local.py | 57 +++++++++++++++++++++++++++++++--------- sumpy/kernel.py | 21 ++++++++++++--- 2 files changed, 61 insertions(+), 17 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index c3ababe3..d0f9881e 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -180,8 +180,9 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): base_kernel = single_valued(knl.get_base_kernel() for knl in kernels) # TODO: use a kernel mapper to find asymptotic info wrapper that is not # the outmost layer - from sumpy.kernel import AsymptoticallyInformedKernel - using_qbmax = isinstance(base_kernel, AsymptoticallyInformedKernel) + from sumpy.kernel import AsymptoticScalingGetter + asg = AsymptoticScalingGetter() + using_qbmax = (asg(base_kernel) is not None) if using_qbmax and bvec is not None: # bvec is None when, e.g., @@ -261,17 +262,21 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): bvec_scaled = [b*rscale**-1 for b in bvec] from sumpy.tools import mi_power, mi_factorial - from sumpy.kernel import AsymptoticScalingDetector, DerivativeCounter - asd = AsymptoticScalingDetector() + from sumpy.kernel import AsymptoticScalingGetter, DerivativeCounter + asg = AsymptoticScalingGetter() dcr = DerivativeCounter() - using_qbmax = (asd(kernel).get_base_kernel().dim == 0) + scaling_expr = asg(kernel) nder = dcr(kernel) - if using_qbmax: - # FIXME: do target derivative for QBMAX here - # NOTE: only first order target derivative is implemented + if (scaling_expr is not None) and (nder > 0): + # Target derivatives for qbmax require postprocessing with + # the Leibniz formula. + # + # NOTE: only first order target derivative is implemented. + # The kernel passed to this method during QBX should only have + # target derivatives. if nder > 1: raise NotImplementedError - # + # Potential = (int[f(x)*G(x,y)dx] * g(y) * g^{-1}(y))' # = Expand{int[f(x)*G(x,y)dx] * g(y)}' / g(y) # - Expand{int[f(x)*G(x,y)dx] * g(y)} * g'(y)/g^2(y) @@ -281,12 +286,38 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # ^: A_n are coeffs of usual target derivative # ^: B_n are original coeffs - # get the map from B_n to A_n + # get the map coefficients at the target + expr_dict = {(0,)*self.dim: 1} + expr_dict = kernel.get_derivative_transformation_at_target(expr_dict) + + def orig_term(coeff, mi): # sympy expr for B_mi + # B_mi = (taylor_coeff * rscale^mi) * (((x-center)/rscale)^mi/mi!) + # ^: coeff + return (coeff * mi_power(bvec_scaled, mi, evaluate=False) + / mi_factorial(mi)) + + def tder_term(coeff, mi): # sympy expr for A_mi + expr = orig_term(coeff, mi) + return sum( + dcoeff * kernel._diff(expr, bvec, mi) + for mi, dcoeff in expr_dict.items()) + + # get sympy expr for g'/g + g = 1 / scaling_expr + gpog = sum( + coeff * kernel._diff(g, bvec, mi) / g + for mi, coeff in expr_dict.items()) + + def leibniz_term(coeff, mi): + ami = tder_term(coeff, mi) + bmi = orig_term(coeff, mi) + return ami - bmi * gpog - # compute g'(y)/g(y) + result = sum( + leibniz_term(coeff, mi) + for coeff, mi in zip( + evaluated_coeffs, self.get_full_coefficient_identifiers())) - # assemble result - pass else: result = sum( coeff diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 7d525d67..aa473251 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1708,17 +1708,30 @@ class AsymptoticScalingRemover(KernelIdentityMapper): return self.rec(kernel.inner_kernel) -class AsymptoticScalingDetector(KernelIdentityMapper): - """Swaps base kernel with a placeholder to signal existence of asymptotic - information. This is needed because wrappers are typically ordered from +class AsymptoticScalingGetter(KernelCombineMapper): + """This is needed because wrappers are typically ordered from inside out as: - base kernel - source derivatives - asymptotic information - target derivatives """ + def combine(self, values): + return values[0] + + def map_expression_kernel(self, kernel): + return None + + 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_line_of_compression_kernel = map_expression_kernel + map_elasticity_kernel = map_expression_kernel + map_stresslet_kernel = map_expression_kernel + def map_asymptotically_informed_kernel(self, kernel): - return ExpressionKernel(0, None, None, None) + return kernel.get_scaling_expression(None) class ExpaniosnClassReplacer(KernelIdentityMapper): -- GitLab From 665124eee2aaefb78e8db7a297136f36ff6b3a95 Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 10 Jul 2021 12:54:35 -0500 Subject: [PATCH 27/32] Add switch for qbmax in local epansions --- sumpy/e2p.py | 5 +++-- sumpy/expansion/local.py | 6 ++++-- sumpy/expansion/multipole.py | 2 +- sumpy/qbx.py | 11 +++++++++-- 4 files changed, 17 insertions(+), 7 deletions(-) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index d68de9f3..3033f0c4 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -77,7 +77,7 @@ class E2PBase(KernelCacheWrapper): self.dim = expansion.dim self._vector_names = {"b"} - def get_loopy_insns_and_result_names(self): + def get_loopy_insns_and_result_names(self, use_qbmax=False): from sumpy.symbolic import make_sym_vector bvec = make_sym_vector("b", self.dim) @@ -92,7 +92,8 @@ class E2PBase(KernelCacheWrapper): result_names = [ sac.assign_unique("result_%d_p" % i, - self.expansion.evaluate(knl, coeff_exprs, bvec, rscale, sac=sac)) + self.expansion.evaluate(knl, coeff_exprs, bvec, rscale, sac=sac, + use_qbmax=use_qbmax)) for i, knl in enumerate(self.kernels) ] diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index d0f9881e..9e1b526a 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -251,7 +251,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): return self.coefficients_from_source_vec((kernel,), avec, bvec, rscale=rscale, weights=(1,), sac=sac) - def evaluate(self, kernel, coeffs, bvec, rscale, sac=None): + def evaluate(self, kernel, coeffs, bvec, rscale, sac=None, use_qbmax=False): if not self.use_rscale: rscale = 1 @@ -267,7 +267,9 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): dcr = DerivativeCounter() scaling_expr = asg(kernel) nder = dcr(kernel) - if (scaling_expr is not None) and (nder > 0): + # When using FMM, we want to turn off qbmax for the far field, + # thus `use_qbmax` must be explicitly set to `True` for qbmax. + if use_qbmax and (scaling_expr is not None) and (nder > 0): # Target derivatives for qbmax require postprocessing with # the Leibniz formula. # diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 408b90a9..387792cf 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -85,7 +85,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): return self.coefficients_from_source_vec((kernel,), avec, bvec, rscale, (1,), sac=sac) - def evaluate(self, kernel, coeffs, bvec, rscale, sac=None): + def evaluate(self, kernel, coeffs, bvec, rscale, sac=None, use_qbmax=False): if not self.use_rscale: rscale = 1 diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 3cb8fde6..36d3b10d 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -69,7 +69,7 @@ def stringify_expn_index(i): class LayerPotentialBase(KernelComputation, KernelCacheWrapper): def __init__(self, ctx, expansion, strength_usage=None, value_dtypes=None, name=None, device=None, - source_kernels=None, target_kernels=None): + source_kernels=None, target_kernels=None, _use_qbmax=False): from sumpy.kernel import ExpaniosnClassReplacer ecr = ExpaniosnClassReplacer(type(expansion)) @@ -84,6 +84,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): self.dim = single_valued(knl.dim for knl in self.target_kernels) self.expansion = expansion self._vector_names = {"a", "b"} + self._use_qbmax = _use_qbmax def get_cache_key(self): return (type(self).__name__, self.expansion, tuple(self.target_kernels), @@ -120,8 +121,14 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): coefficients[self.expansion.get_storage_index(i)])) for i in self.expansion.get_coefficient_identifiers()] + if self._use_qbmax is None or self._use_qbmax: + use_qbmax = True + else: + use_qbmax = False + return sac.assign_unique("expn%d_result" % expansion_nr, - self.expansion.evaluate(tgt_knl, assigned_coeffs, bvec, rscale)) + self.expansion.evaluate(tgt_knl, assigned_coeffs, bvec, rscale, + use_qbmax=use_qbmax)) def get_loopy_insns_and_result_names(self): from sumpy.symbolic import make_sym_vector -- GitLab From 88697713f851452d611386ece64bdc8f95223176 Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 10 Jul 2021 12:59:39 -0500 Subject: [PATCH 28/32] Fix expansion.evaluate() interface --- sumpy/expansion/__init__.py | 2 +- sumpy/expansion/local.py | 4 ++-- sumpy/expansion/multipole.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 8bfb85e9..41071aaf 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -141,7 +141,7 @@ class ExpansionBase: result[i] += weight * coeffs[i] return result - def evaluate(self, kernel, coeffs, bvec, rscale, sac=None): + def evaluate(self, kernel, coeffs, bvec, rscale, sac=None, use_qbmax=False): """ :return: a :mod:`sympy` expression corresponding to the evaluated expansion with the coefficients diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 9e1b526a..3097a1c4 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -142,7 +142,7 @@ class LineTaylorLocalExpansion(LocalExpansionBase): .subs("tau", 0) for i in self.get_coefficient_identifiers()] - def evaluate(self, tgt_kernel, coeffs, bvec, rscale, sac=None): + def evaluate(self, tgt_kernel, coeffs, bvec, rscale, sac=None, use_qbmax=False): # no point in heeding rscale here--just ignore it from pytools import factorial return sym.Add(*( @@ -720,7 +720,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): * sym.exp(sym.I * c * source_angle_rel_center), avec) for c in self.get_coefficient_identifiers()] - def evaluate(self, kernel, coeffs, bvec, rscale, sac=None): + def evaluate(self, kernel, coeffs, bvec, rscale, sac=None, use_qbmax=False): if not self.use_rscale: rscale = 1 diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 387792cf..4cdc0e71 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -420,7 +420,7 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): avec) for c in self.get_coefficient_identifiers()] - def evaluate(self, kernel, coeffs, bvec, rscale, sac=None): + def evaluate(self, kernel, coeffs, bvec, rscale, sac=None, use_qbmax=False): if not self.use_rscale: rscale = 1 -- GitLab From e0161aebea41b2ed7ef4838a47b321b94ae7276c Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 11 Jul 2021 12:23:47 -0500 Subject: [PATCH 29/32] Add sumpy qbmax support for volume expansion (self eval), add tests --- sumpy/qbx.py | 42 +++++++++++++++++++++++++++-------- test/test_qbx.py | 58 +++++++++++++++++++++++++++++++++++------------- 2 files changed, 75 insertions(+), 25 deletions(-) diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 36d3b10d..cb5e70a0 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -241,8 +241,6 @@ class LayerPotential(LayerPotentialBase): @memoize_method def get_kernel(self): - loopy_insns, result_names = self.get_loopy_insns_and_result_names() - kernel_exprs = self.get_kernel_exprs(result_names) arguments = ( self.get_default_src_tgt_arguments() + [lp.GlobalArg("strength_%d" % i, @@ -254,12 +252,15 @@ class LayerPotential(LayerPotentialBase): sats = [] sat_insns_data = [] + for iknl, knl in enumerate(self.target_kernels): - try: - knl = knl.replace_expansion_class(type(self.expansion)) - sat_expr = knl.get_scaling_expression(None) - except AttributeError: - sat_expr = None + from sumpy.kernel import AsymptoticScalingGetter, ExpaniosnClassReplacer + ecr = ExpaniosnClassReplacer(type(self.expansion)) + knl = ecr(knl) + + asg = AsymptoticScalingGetter() + sat_expr = asg(knl) + if sat_expr: # inject the multiplicative scales of each output kernel sats.append(f"sat_{iknl} * ") @@ -267,7 +268,26 @@ class LayerPotential(LayerPotentialBase): else: sats.append("") + qbmax_local_temp_insns = [] if sat_insns_data: + # Volume expansion needs normal vector information. + # Though using volume expansion, it only evaluates self interaction. + # Centers and targets are 1-to-1. + from sumpy.expansion import VolumeTaylorExpansionBase + if isinstance(self.expansion, VolumeTaylorExpansionBase): + qbmax_local_temp_insns = [""" + <> normal[idim] = \ + nodes_of_qbx_balls_tangency[idim, itgt//2] - \ + center[idim, itgt] {dup=idim} + """] + + arguments += ( + lp.GlobalArg("nodes_of_qbx_balls_tangency", None, + shape=(self.dim, "ntargets"), dim_tags="sep,C"), + ) + + self._vector_names.add("normal") + from sumpy.codegen import to_loopy_insns sat_insns = to_loopy_insns( sat_insns_data, @@ -277,6 +297,9 @@ class LayerPotential(LayerPotentialBase): else: sat_insns = [] + loopy_insns, result_names = self.get_loopy_insns_and_result_names() + kernel_exprs = self.get_kernel_exprs(result_names) + loopy_knl = lp.make_kernel([""" {[itgt, isrc, idim]: \ 0 <= itgt < ntargets and \ @@ -290,8 +313,8 @@ class LayerPotential(LayerPotentialBase): + ["<> rscale = expansion_radii[itgt]"] + [f"<> strength_{i}_isrc = strength_{i}[isrc]" for i in range(self.strength_count)] + + sat_insns + qbmax_local_temp_insns + loopy_insns + kernel_exprs - + sat_insns + [""" result_{i}[itgt] = knl_{i}_scaling * {sat} \ simul_reduce(sum, isrc, pair_result_{i}) \ @@ -311,6 +334,8 @@ class LayerPotential(LayerPotentialBase): for knl in self.target_kernels + self.source_kernels: loopy_knl = knl.prepare_loopy_kernel(loopy_knl) + print(loopy_knl) + return loopy_knl def __call__(self, queue, targets, sources, centers, strengths, expansion_radii, @@ -319,7 +344,6 @@ class LayerPotential(LayerPotentialBase): :arg strengths: are required to have area elements and quadrature weights already multiplied in. """ - knl = self.get_cached_optimized_kernel( targets_is_obj_array=is_obj_array_like(targets), sources_is_obj_array=is_obj_array_like(sources), diff --git a/test/test_qbx.py b/test/test_qbx.py index 0eedda27..d23f74d7 100644 --- a/test/test_qbx.py +++ b/test/test_qbx.py @@ -45,6 +45,8 @@ asexpr = sp.exp(-dist) # no practical use, just to test the machinary (VolumeTaylorLocalExpansion, LaplaceKernel(2)), (LineTaylorLocalExpansion, AsymptoticallyInformedKernel( LaplaceKernel(2), asexpr)), + (VolumeTaylorLocalExpansion, AsymptoticallyInformedKernel( + LaplaceKernel(2), asexpr)), ]) def test_direct_qbx_vs_eigval(ctx_factory, expn_class, lknl): """This evaluates a single layer potential on a circle using a known @@ -60,8 +62,11 @@ def test_direct_qbx_vs_eigval(ctx_factory, expn_class, lknl): from sumpy.qbx import LayerPotential - lpot = LayerPotential(ctx, expansion=expn_class(lknl, order), - target_kernels=(lknl,), source_kernels=(lknl,)) + use_qbmax = isinstance(lknl, AsymptoticallyInformedKernel) + lpot = LayerPotential( + ctx, expansion=expn_class(lknl, order), + target_kernels=(lknl,), source_kernels=(lknl,), + _use_qbmax=use_qbmax) mode_nr = 25 @@ -89,9 +94,16 @@ def test_direct_qbx_vs_eigval(ctx_factory, expn_class, lknl): expansion_radii = np.ones(n) * radius + if use_qbmax and (expn_class is VolumeTaylorLocalExpansion): + qbmax_extra_kwargs = {"nodes_of_qbx_balls_tangency": unit_circle} + else: + qbmax_extra_kwargs = dict() + strengths = (sigma * h,) - evt, (result_qbx,) = lpot(queue, targets, sources, centers, strengths, - expansion_radii=expansion_radii) + evt, (result_qbx,) = lpot( + queue, targets, sources, centers, strengths, + expansion_radii=expansion_radii, + **qbmax_extra_kwargs) eocrec.add_data_point(h, np.max(np.abs(result_ref - result_qbx))) @@ -101,11 +113,15 @@ def test_direct_qbx_vs_eigval(ctx_factory, expn_class, lknl): assert eocrec.order_estimate() > order - slack -@pytest.mark.parametrize("expn_class", [ - LineTaylorLocalExpansion, - VolumeTaylorLocalExpansion, +@pytest.mark.parametrize("expn_class, lknl", [ + (LineTaylorLocalExpansion, LaplaceKernel(2)), + (VolumeTaylorLocalExpansion, LaplaceKernel(2)), + (LineTaylorLocalExpansion, AsymptoticallyInformedKernel( + LaplaceKernel(2), asexpr)), + (VolumeTaylorLocalExpansion, AsymptoticallyInformedKernel( + LaplaceKernel(2), asexpr)), ]) -def test_direct_qbx_vs_eigval_with_tgt_deriv(ctx_factory, expn_class): +def test_direct_qbx_vs_eigval_with_tgt_deriv(ctx_factory, expn_class, lknl): """This evaluates a single layer potential on a circle using a known eigenvalue/eigenvector combination. """ @@ -115,17 +131,22 @@ def test_direct_qbx_vs_eigval_with_tgt_deriv(ctx_factory, expn_class): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - from sumpy.kernel import LaplaceKernel, AxisTargetDerivative - lknl = LaplaceKernel(2) + from sumpy.kernel import AxisTargetDerivative 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,)) + # explicitly enable _use_qbmax when directly using LayerPotential + use_qbmax = isinstance(lknl, AsymptoticallyInformedKernel) + lpot_dx = LayerPotential( + ctx, expansion=expn_class(lknl, order), + target_kernels=(AxisTargetDerivative(0, lknl),), source_kernels=(lknl,), + _use_qbmax=use_qbmax) + lpot_dy = LayerPotential( + ctx, expansion=expn_class(lknl, order), + target_kernels=(AxisTargetDerivative(1, lknl),), source_kernels=(lknl,), + _use_qbmax=use_qbmax) mode_nr = 15 @@ -156,10 +177,15 @@ def test_direct_qbx_vs_eigval_with_tgt_deriv(ctx_factory, expn_class): strengths = (sigma * h,) + if use_qbmax and (expn_class is VolumeTaylorLocalExpansion): + qbmax_extra_kwargs = {"nodes_of_qbx_balls_tangency": unit_circle} + else: + qbmax_extra_kwargs = dict() + evt, (result_qbx_dx,) = lpot_dx(queue, targets, sources, centers, strengths, - expansion_radii=expansion_radii) + expansion_radii=expansion_radii, **qbmax_extra_kwargs) evt, (result_qbx_dy,) = lpot_dy(queue, targets, sources, centers, strengths, - expansion_radii=expansion_radii) + expansion_radii=expansion_radii, **qbmax_extra_kwargs) normals = unit_circle result_qbx = normals[0] * result_qbx_dx + normals[1] * result_qbx_dy -- GitLab From ebfcb080f7d86bce5e826f3978f4b81a4a3135b3 Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 11 Jul 2021 17:35:21 -0500 Subject: [PATCH 30/32] Remove print statement --- sumpy/qbx.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/sumpy/qbx.py b/sumpy/qbx.py index cb5e70a0..be9d55a6 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -334,8 +334,6 @@ class LayerPotential(LayerPotentialBase): for knl in self.target_kernels + self.source_kernels: loopy_knl = knl.prepare_loopy_kernel(loopy_knl) - print(loopy_knl) - return loopy_knl def __call__(self, queue, targets, sources, centers, strengths, expansion_radii, -- GitLab From 35eeeb075a17e743455077d30694a5a3b0163fe6 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 12 Jul 2021 14:07:29 -0500 Subject: [PATCH 31/32] Refactor source derivative --- sumpy/expansion/local.py | 38 +++++++++++--------------------------- sumpy/kernel.py | 16 ++++++++++++---- sumpy/tools.py | 23 +++++++++++++++++++++++ 3 files changed, 46 insertions(+), 31 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 3097a1c4..aece0bb8 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -182,54 +182,38 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # the outmost layer from sumpy.kernel import AsymptoticScalingGetter asg = AsymptoticScalingGetter() - using_qbmax = (asg(base_kernel) is not None) - if using_qbmax and bvec is not None: - # bvec is None when, e.g., - # 1) using FMM without QBMAX - # 2) using FMM with QBMAX, during P2L - # + # bvec is None when + # 1) when using FMM without QBMAX + # 2) when using FMM with QBMAX, but during GIGAQBX's far-field P2L + using_qbmax = (asg(base_kernel) is not None) and (bvec is not None) + + if using_qbmax: # With a directional expansion (like line taylor), we write: # source to target = ( # source to center (avec) + center to target (bvec)). - # Kernel evaluated wrt rad, postprocessed wrt avec, + # Kernel evaluated wrt rad = avec + bvec, postprocessed wrt avec, # then take derivatives along the direction of bvec, and # finally let bvec=0 - rad = avec + bvec + base_taker = base_kernel.get_expression(avec + bvec) post_diff_subs = {bvc: 0 for bvc in bvec} else: base_taker = base_kernel.get_derivative_taker(avec, rscale, sac) post_diff_subs = dict() result = [0]*len(self) - for knl, weight in zip(kernels, weights): # apply source derivatives, if any - if using_qbmax and bvec is not None: + if using_qbmax: # Since DifferentiatedExprDerivativeTaker does not work with # anisotropic expansions, using vanilla expression here is # necessary. That is, the base_taker is directly constructed # from \partial_{n(y)}G(x-y) from sumpy.tools import ExprDerivativeTaker - ppkernel = knl.postprocess_at_source(knl.get_expression(rad), avec) + ppkernel = knl.postprocess_at_source(base_taker, avec, bvec) base_taker = ExprDerivativeTaker(ppkernel, bvec, rscale, sac) - from pytools.tag import tag_dataclass - - @tag_dataclass - class AnisotropicDifferentiatedExprDerivativeTaker: - taker: ExprDerivativeTaker - - def diff(self, mi, save_intermediate=lambda x: x, subs=None): - result = self.taker.diff(mi) - if subs: - result = result.subs(subs) - return result * save_intermediate(1) - - taker = AnisotropicDifferentiatedExprDerivativeTaker(base_taker) - - else: - taker = knl.postprocess_at_source(base_taker, avec) + taker = knl.postprocess_at_source(base_taker, avec, bvec) # Following is a hack to make sure cse works. if 1: diff --git a/sumpy/kernel.py b/sumpy/kernel.py index aa473251..fb58346b 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -211,7 +211,7 @@ class Kernel: expr = expr.diff(vec[i], mi[i]) return expr - def postprocess_at_source(self, expr, avec): + def postprocess_at_source(self, expr, avec, bvec=None): """Transform a kernel evaluation or expansion expression in a place where the vector a (something - source) is known. ("something" may be an expansion center or a target.) @@ -219,12 +219,20 @@ class Kernel: The typical use of this function is to apply source-variable derivatives to the kernel. """ - from sumpy.tools import (ExprDerivativeTaker, - DifferentiatedExprDerivativeTaker) + from sumpy.tools import ( + ExprDerivativeTaker, + DifferentiatedExprDerivativeTaker, + AnisotropicDifferentiatedExprDerivativeTaker) + expr_dict = {(0,)*self.dim: 1} expr_dict = self.get_derivative_transformation_at_source(expr_dict) + if isinstance(expr, ExprDerivativeTaker): - return DifferentiatedExprDerivativeTaker(expr, expr_dict) + if bvec is not None: + # qbmax + return AnisotropicDifferentiatedExprDerivativeTaker(expr) + else: + return DifferentiatedExprDerivativeTaker(expr, expr_dict) result = 0 for mi, coeff in expr_dict.items(): diff --git a/sumpy/tools.py b/sumpy/tools.py index 58564d2c..d6a11c9e 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -430,6 +430,29 @@ class DifferentiatedExprDerivativeTaker: return result * save_intermediate(1 / self.taker.rscale ** max_order) + +@tag_dataclass +class AnisotropicDifferentiatedExprDerivativeTaker: + """Implements the :class:`ExprDerivativeTaker` interface for an expression + that is itself a linear combination of derivatives of a base expression. To + take the actual derivatives, it makes use of an underlying derivative taker + *taker*. + + Unlike :class:`DifferentiatedExprDerivativeTaker`, the derivatives inside + the taker are already evaluated. Because for anisotropic (target-dependent) + qbmax expansions, derivative transforms are nonlinear. + + .. attribute:: taker + A :class:`ExprDerivativeTaker` for the base expression. + """ + taker: ExprDerivativeTaker + + def diff(self, mi, save_intermediate=lambda x: x, subs=None): + result = self.taker.diff(mi) + if subs: + result = result.subs(subs) + return result * save_intermediate(1) + # }}} -- GitLab From 93464217f3a0e1d363cc22c9356f1285907a88b5 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 20 Jul 2021 15:06:11 -0500 Subject: [PATCH 32/32] Fix pytential tests with (use_fmm=False, use_qbmax=False) --- sumpy/expansion/local.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index aece0bb8..f66d283a 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -200,6 +200,8 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): else: base_taker = base_kernel.get_derivative_taker(avec, rscale, sac) post_diff_subs = dict() + # to signal postprocess_at_source that we are not using qbmax + bvec = None result = [0]*len(self) for knl, weight in zip(kernels, weights): -- GitLab