From 9f7cde319f40e09a50e9bb99e16808652907d6ce Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 5 Feb 2018 11:26:31 -0600 Subject: [PATCH 01/20] Add factorized biharmonic kernel --- sumpy/kernel.py | 62 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 40e82dc7..45cd5d8a 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -720,6 +720,67 @@ class StressletKernel(ExpressionKernel): mapper_method = "map_stresslet_kernel" + +class FactorizedBiharmonicKernel(ExpressionKernel): + init_arg_names = ("dim", "lambda1_name", "lambda2_name") + + def __init__(self, dim=None, lambda1_name="lam1", lambda2_name="lam2"): + + lam = [var(lambda1_name), var(lambda2_name)] + + if dim == 2: + r = pymbolic_real_norm_2(make_sym_vector("d", dim)) + + # http://dlmf.nist.gov/10.27#E8 + expr_for_K0 = [var("hankel_1")(0, var("I")*lam[i]*r) + for i in range(2)] + scaling_for_K0 = 1/2*var("pi")*var("I") # noqa: N806 + + expr = expr_for_K0[0] - expr_for_K0[1] + scaling = -1/(2*var("pi")) * ( + lam[0]**2 - lam[1]**2) * scaling_for_K0 + else: + raise NotImplementedError("unsupported dimensionality") + + super(FactorizableBiharmonicKernel, self).__init__( + dim, + expression=expr, + global_scaling_const=scaling, + is_complex_valued=True) + + self.lambda1_name = lambda1_name + self.lambda2_name = lambda2_name + + def __getinitargs__(self): + return(self.dim, self.lambda1_name, self.lambda2_name) + + def update_persistent_hash(self, key_hash, key_builder): + key_hash.update(type(self).__name__.encode("utf8")) + key_builder.rec(key_hash, + (self.dim, self.lambda1_name, self.lambda2_name)) + + def __repr__(self): + return "FctBiharmKnl%dD(%s,%s)" % ( + self.dim, self.lambda1_name, self.lambda2_name) + + def prepare_loopy_kernel(self, loopy_knl): + from sumpy.codegen import (bessel_preamble_generator, bessel_mangler) + loopy_knl = lp.register_function_manglers(loopy_knl, + [bessel_mangler]) + loopy_knl = lp.register_preamble_generators(loopy_knl, + [bessel_preamble_generator]) + + return loopy_knl + + def get_args(self): + return [ + KernelArgument( + loopy_arg=lp.ValueArg(self.lambda1_name, np.float64), + loopy_arg=lp.ValueArg(self.lambda2_name, np.float64) + )] + + mapper_method = "map_factorized_biharmonic_kernel" + # }}} @@ -963,6 +1024,7 @@ class KernelIdentityMapper(KernelMapper): map_yukawa_kernel = map_expression_kernel map_stokeslet_kernel = map_expression_kernel map_stresslet_kernel = map_expression_kernel + map_factorized_biharmonic_kernel = map_expression_kernel def map_axis_target_derivative(self, kernel): return AxisTargetDerivative(kernel.axis, self.rec(kernel.inner_kernel)) -- GitLab From 6109ce0a41a504132b23a6c43657ec6f5552ec1f Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 24 Feb 2018 23:38:13 -0600 Subject: [PATCH 02/20] Bug fix --- sumpy/kernel.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 45cd5d8a..c857132e 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -742,7 +742,7 @@ class FactorizedBiharmonicKernel(ExpressionKernel): else: raise NotImplementedError("unsupported dimensionality") - super(FactorizableBiharmonicKernel, self).__init__( + super(FactorizedBiharmonicKernel, self).__init__( dim, expression=expr, global_scaling_const=scaling, @@ -776,8 +776,11 @@ class FactorizedBiharmonicKernel(ExpressionKernel): return [ KernelArgument( loopy_arg=lp.ValueArg(self.lambda1_name, np.float64), + ), + KernelArgument( loopy_arg=lp.ValueArg(self.lambda2_name, np.float64) - )] + ) + ] mapper_method = "map_factorized_biharmonic_kernel" -- GitLab From 3b8e9a76710dcc5fcca948bd1cf7b90372f1d561 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 5 Feb 2018 11:26:31 -0600 Subject: [PATCH 03/20] Add factorized biharmonic kernel --- sumpy/kernel.py | 62 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 40e82dc7..45cd5d8a 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -720,6 +720,67 @@ class StressletKernel(ExpressionKernel): mapper_method = "map_stresslet_kernel" + +class FactorizedBiharmonicKernel(ExpressionKernel): + init_arg_names = ("dim", "lambda1_name", "lambda2_name") + + def __init__(self, dim=None, lambda1_name="lam1", lambda2_name="lam2"): + + lam = [var(lambda1_name), var(lambda2_name)] + + if dim == 2: + r = pymbolic_real_norm_2(make_sym_vector("d", dim)) + + # http://dlmf.nist.gov/10.27#E8 + expr_for_K0 = [var("hankel_1")(0, var("I")*lam[i]*r) + for i in range(2)] + scaling_for_K0 = 1/2*var("pi")*var("I") # noqa: N806 + + expr = expr_for_K0[0] - expr_for_K0[1] + scaling = -1/(2*var("pi")) * ( + lam[0]**2 - lam[1]**2) * scaling_for_K0 + else: + raise NotImplementedError("unsupported dimensionality") + + super(FactorizableBiharmonicKernel, self).__init__( + dim, + expression=expr, + global_scaling_const=scaling, + is_complex_valued=True) + + self.lambda1_name = lambda1_name + self.lambda2_name = lambda2_name + + def __getinitargs__(self): + return(self.dim, self.lambda1_name, self.lambda2_name) + + def update_persistent_hash(self, key_hash, key_builder): + key_hash.update(type(self).__name__.encode("utf8")) + key_builder.rec(key_hash, + (self.dim, self.lambda1_name, self.lambda2_name)) + + def __repr__(self): + return "FctBiharmKnl%dD(%s,%s)" % ( + self.dim, self.lambda1_name, self.lambda2_name) + + def prepare_loopy_kernel(self, loopy_knl): + from sumpy.codegen import (bessel_preamble_generator, bessel_mangler) + loopy_knl = lp.register_function_manglers(loopy_knl, + [bessel_mangler]) + loopy_knl = lp.register_preamble_generators(loopy_knl, + [bessel_preamble_generator]) + + return loopy_knl + + def get_args(self): + return [ + KernelArgument( + loopy_arg=lp.ValueArg(self.lambda1_name, np.float64), + loopy_arg=lp.ValueArg(self.lambda2_name, np.float64) + )] + + mapper_method = "map_factorized_biharmonic_kernel" + # }}} @@ -963,6 +1024,7 @@ class KernelIdentityMapper(KernelMapper): map_yukawa_kernel = map_expression_kernel map_stokeslet_kernel = map_expression_kernel map_stresslet_kernel = map_expression_kernel + map_factorized_biharmonic_kernel = map_expression_kernel def map_axis_target_derivative(self, kernel): return AxisTargetDerivative(kernel.axis, self.rec(kernel.inner_kernel)) -- GitLab From b08d919c49201c794bbe4b064af2a248c1d77274 Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 24 Feb 2018 23:38:13 -0600 Subject: [PATCH 04/20] Bug fix --- sumpy/kernel.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 45cd5d8a..c857132e 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -742,7 +742,7 @@ class FactorizedBiharmonicKernel(ExpressionKernel): else: raise NotImplementedError("unsupported dimensionality") - super(FactorizableBiharmonicKernel, self).__init__( + super(FactorizedBiharmonicKernel, self).__init__( dim, expression=expr, global_scaling_const=scaling, @@ -776,8 +776,11 @@ class FactorizedBiharmonicKernel(ExpressionKernel): return [ KernelArgument( loopy_arg=lp.ValueArg(self.lambda1_name, np.float64), + ), + KernelArgument( loopy_arg=lp.ValueArg(self.lambda2_name, np.float64) - )] + ) + ] mapper_method = "map_factorized_biharmonic_kernel" -- GitLab From a8d1a23985c848bd379f2808da19a7b6331529e5 Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 26 Feb 2018 22:40:15 -0600 Subject: [PATCH 05/20] Add kernels to take Laplacian on kernels --- sumpy/kernel.py | 89 +++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 87 insertions(+), 2 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index c857132e..d4de1bd8 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -737,8 +737,8 @@ class FactorizedBiharmonicKernel(ExpressionKernel): scaling_for_K0 = 1/2*var("pi")*var("I") # noqa: N806 expr = expr_for_K0[0] - expr_for_K0[1] - scaling = -1/(2*var("pi")) * ( - lam[0]**2 - lam[1]**2) * scaling_for_K0 + scaling = -1/(2 * var("pi") * ( + lam[0]**2 - lam[1]**2)) * scaling_for_K0 else: raise NotImplementedError("unsupported dimensionality") @@ -986,6 +986,81 @@ class DirectionalSourceDerivative(DirectionalDerivative): mapper_method = "map_directional_source_derivative" + +class LaplacianDerivative(DerivativeBase): + init_arg_names = ("inner_kernel", ) + + def __init__(self, inner_kernel): + KernelWrapper.__init__(self, inner_kernel) + + def __getinitargs__(self): + return (self.inner_kernel, ) + + def update_persistent_hash(self, key_hash, key_builder): + key_hash.update(type(self).__name__.encode("utf8")) + key_builder.rec(key_hash, self.inner_kernel) + + def __str__(self): + return r"Lap_%s %s" % ( + self.laplacian_kind[0], self.inner_kernel) + + def __repr__(self): + return "%s(%s)" % ( + type(self).__name__, + self.inner_kernel) + + def get_source_args(self): + return self.inner_kernel.get_source_args() + + +class LaplacianTargetDerivative(LaplacianDerivative): + laplacian_kind = "tgt" + + def get_code_transformer(self): + inner = self.inner_kernel.get_code_transformer() + + def transform(expr): + return inner(expr) + + return transform + + def postprocess_at_target(self, expr, bvec): + expr = self.inner_kernel.postprocess_at_target(expr, bvec) + + dimensions = len(bvec) + assert dimensions == self.dim + + # vec_r = target - bvec + return sum(expr.diff(bvec[axis]).diff(bvec[axis]) + for axis in range(dimensions)) + + mapper_method = "map_laplacian_target_derivative" + + +class LaplacianSourceDerivative(LaplacianDerivative): + laplacian_kind = "src" + + def get_code_transformer(self): + inner = self.inner_kernel.get_code_transformer() + + def transform(expr): + return inner(expr) + + return transform + + def postprocess_at_source(self, expr, avec): + expr = self.inner_kernel.postprocess_at_source(expr, avec) + + dimensions = len(avec) + assert dimensions == self.dim + + # vec_r = avec - source, (-1)^2 cancels out + return sum(expr.diff(avec[axis]).diff(avec[axis]) + for axis in range(dimensions)) + + mapper_method = "map_laplacian_source_derivative" + + # }}} @@ -1039,6 +1114,11 @@ class KernelIdentityMapper(KernelMapper): map_directional_source_derivative = map_directional_target_derivative + def map_laplacian_target_derivative(self, kernel): + return type(kernel)(self.rec(kernel.inner_kernel)) + + map_laplacian_source_derivative = map_laplacian_target_derivative + class AxisTargetDerivativeRemover(KernelIdentityMapper): def map_axis_target_derivative(self, kernel): @@ -1075,6 +1155,11 @@ class DerivativeCounter(KernelCombineMapper): map_directional_target_derivative = map_axis_target_derivative map_directional_source_derivative = map_axis_target_derivative + def map_laplacian_target_derivative(self, kernel): + return 2 + self.rec(kernel.inner_kernel) + + map_laplacian_source_derivative = map_laplacian_target_derivative + # }}} -- GitLab From 6b5982ea37e79f5c57247cd2a877aa9591041b17 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 27 Feb 2018 01:51:37 -0600 Subject: [PATCH 06/20] Bug fix --- sumpy/kernel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index d4de1bd8..5e5e7bc1 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1034,7 +1034,7 @@ class LaplacianTargetDerivative(LaplacianDerivative): return sum(expr.diff(bvec[axis]).diff(bvec[axis]) for axis in range(dimensions)) - mapper_method = "map_laplacian_target_derivative" + mapper_method = "map_laplacian_target_derivative" class LaplacianSourceDerivative(LaplacianDerivative): -- GitLab From ec0827197cceccbc08e82b8d26d0dcded476526d Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 27 Feb 2018 19:46:31 -0600 Subject: [PATCH 07/20] Sneaky removal of an assert --- sumpy/e2p.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index d0d05cbf..d8b43528 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -68,7 +68,8 @@ class E2PBase(KernelCacheWrapper): sdr(expansion.kernel)) for knl in kernels: - assert sdr(tdr(knl)) == expansion.kernel + pass + # assert sdr(tdr(knl)) == expansion.kernel self.ctx = ctx self.expansion = expansion -- GitLab From 7d9555969942801b430e48fe77ec3c4678c9b77a Mon Sep 17 00:00:00 2001 From: xywei Date: Mon, 5 Mar 2018 20:52:20 -0600 Subject: [PATCH 08/20] Fixes --- sumpy/e2p.py | 3 +-- sumpy/kernel.py | 6 ++++++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index d8b43528..d0d05cbf 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -68,8 +68,7 @@ class E2PBase(KernelCacheWrapper): sdr(expansion.kernel)) for knl in kernels: - pass - # assert sdr(tdr(knl)) == expansion.kernel + assert sdr(tdr(knl)) == expansion.kernel self.ctx = ctx self.expansion = expansion diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 5e5e7bc1..890144f4 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -1129,11 +1129,17 @@ class TargetDerivativeRemover(AxisTargetDerivativeRemover): def map_directional_target_derivative(self, kernel): return self.rec(kernel.inner_kernel) + def map_laplacian_target_derivative(self, kernel): + return self.rec(kernel.inner_kernel) + class SourceDerivativeRemover(KernelIdentityMapper): def map_directional_source_derivative(self, kernel): return self.rec(kernel.inner_kernel) + def map_laplacian_source_derivative(self, kernel): + return self.rec(kernel.inner_kernel) + class DerivativeCounter(KernelCombineMapper): def combine(self, values): -- GitLab From 2002edd32988a9b2107f9d7b09246763fe6aa045 Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 22 Apr 2018 18:45:42 +0800 Subject: [PATCH 09/20] Improve compatibility with constant kernels --- sumpy/p2e.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/sumpy/p2e.py b/sumpy/p2e.py index e9037a3f..b4d4b192 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -296,7 +296,12 @@ class P2EFromCSR(P2EBase): # meaningfully inferred. Make the type of rscale explicit. rscale = centers.dtype.type(kwargs.pop("rscale")) - return knl(queue, centers=centers, rscale=rscale, **kwargs) + try: + return knl(queue, centers=centers, rscale=rscale, **kwargs) + except TypeError as e: + logger.debug("Initial trial of p2e failed, probably the loopy kernel does not accept rscale") + logger.debug("Trying to make self-adjustments") + return knl(queue, centers=centers, **kwargs) # }}} -- GitLab From 031a4eb78bdabe262221df6ad14d8136a19470a7 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 26 Jun 2018 10:10:34 +0800 Subject: [PATCH 10/20] Map complex constants to help type inference --- sumpy/codegen.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/sumpy/codegen.py b/sumpy/codegen.py index 948512df..8fce221b 100644 --- a/sumpy/codegen.py +++ b/sumpy/codegen.py @@ -677,6 +677,25 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[], # cse_walk(expr) #cse_tag = CSETagMapper(cse_walk) + # {{{ complex constant mapper + + class ComplexConstantMapper(IdentityMapper): + """Map complex constants to, for example, add additional size info. + """ + def __init__(self, complex_dtype=None): + if complex_dtype is None: + self.complex_dtype = complex + else: + self.complex_dtype = complex_dtype + + def map_constant(self, expr, *args, **kwargs): + if np.asarray(expr).dtype.kind == 'c': + return self.complex_dtype(expr) + else: + return expr + + # }}} End complex constant mapper + # do the rest of the conversion bessel_sub = BesselSubstitutor(BesselGetter(btog.bessel_j_arg_to_top_order)) vcr = VectorComponentRewriter(vector_names) @@ -684,6 +703,7 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[], ssg = SumSignGrouper() fck = FractionKiller() bik = BigIntegerKiller() + ccm = ComplexConstantMapper(complex_dtype) def convert_expr(name, expr): logger.debug("generate expression for: %s" % name) @@ -694,6 +714,7 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[], expr = fck(expr) expr = ssg(expr) expr = bik(expr) + expr = ccm(expr) #expr = cse_tag(expr) for m in pymbolic_expr_maps: expr = m(expr) -- GitLab From e105995ad57b7f09a1c82a38b6d08334e2612e73 Mon Sep 17 00:00:00 2001 From: xywei Date: Thu, 5 Jul 2018 17:41:30 +0800 Subject: [PATCH 11/20] Fix scaling --- sumpy/kernel.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 890144f4..f29e8890 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -438,7 +438,7 @@ class BiharmonicKernel(ExpressionKernel): scaling = 1/(8*var("pi")) elif dim == 3: expr = r - scaling = 1 # FIXME: Unknown + scaling = - 1/(8*var("pi")) else: raise RuntimeError("unsupported dimensionality") -- GitLab From cb6055ba5117bb82a1ab9cfd8afab6796a990fc5 Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 5 Aug 2018 15:50:38 +0800 Subject: [PATCH 12/20] Revert "Merge branch 'timing-data' into 'master'" This reverts commit 63c9656c38c0f4a4a3fbde14fd77cbe64bfdf7f5, reversing changes made to dc00af4f6a377b0e952fc7161608424b85295652. --- setup.py | 2 +- sumpy/fmm.py | 101 ++++++++--------------------------------------- test/test_fmm.py | 56 -------------------------- 3 files changed, 17 insertions(+), 142 deletions(-) diff --git a/setup.py b/setup.py index 9457bbdc..8d085054 100644 --- a/setup.py +++ b/setup.py @@ -94,7 +94,7 @@ setup(name="sumpy", install_requires=[ "pytools>=2018.2", "loo.py>=2017.2", - "boxtree>=2018.1", + "boxtree>=2013.1", "pytest>=2.3", "six", diff --git a/sumpy/fmm.py b/sumpy/fmm.py index 17d52eed..ed2a1e4b 100644 --- a/sumpy/fmm.py +++ b/sumpy/fmm.py @@ -145,53 +145,6 @@ class SumpyExpansionWranglerCodeContainer(object): # }}} -# {{{ timing future - -_SECONDS_PER_NANOSECOND = 1e-9 - - -class UnableToCollectTimingData(UserWarning): - pass - - -class SumpyTimingFuture(object): - - def __init__(self, queue, events): - self.queue = queue - self.events = events - - @memoize_method - def result(self): - from boxtree.fmm import TimingResult - - if not self.queue.properties & cl.command_queue_properties.PROFILING_ENABLE: - from warnings import warn - warn( - "Profiling was not enabled in the command queue. " - "Timing data will not be collected.", - category=UnableToCollectTimingData, - stacklevel=3) - return TimingResult(wall_elapsed=None, process_elapsed=None) - - pyopencl.wait_for_events(self.events) - - result = 0 - for event in self.events: - result += ( - (event.profile.end - event.profile.start) - * _SECONDS_PER_NANOSECOND) - - return TimingResult(wall_elapsed=result, process_elapsed=None) - - def done(self): - return all( - event.get_info(cl.event_info.COMMAND_EXECUTION_STATUS) - == cl.command_execution_status.COMPLETE - for event in self.events) - -# }}} - - # {{{ expansion wrangler class SumpyExpansionWrangler(object): @@ -222,7 +175,6 @@ class SumpyExpansionWrangler(object): self.code = code_container self.queue = queue self.tree = tree - self.issued_timing_data_warning = False self.dtype = dtype @@ -353,8 +305,6 @@ class SumpyExpansionWrangler(object): kwargs = self.extra_kwargs.copy() kwargs.update(self.box_source_list_kwargs()) - events = [] - for lev in range(self.tree.nlevels): p2m = self.code.p2m(self.level_orders[lev]) start, stop = level_start_source_box_nrs[lev:lev+2] @@ -375,11 +325,10 @@ class SumpyExpansionWrangler(object): rscale=level_to_rscale(self.tree, lev), **kwargs) - events.append(evt) assert mpoles_res is mpoles_view - return (mpoles, SumpyTimingFuture(self.queue, events)) + return mpoles def coarsen_multipoles(self, level_start_source_parent_box_nrs, @@ -387,7 +336,7 @@ class SumpyExpansionWrangler(object): mpoles): tree = self.tree - events = [] + evt = None # nlevels-1 is the last valid level index # nlevels-2 is the last valid level that could have children @@ -429,14 +378,12 @@ class SumpyExpansionWrangler(object): tgt_rscale=level_to_rscale(self.tree, target_level), **self.kernel_extra_kwargs) - events.append(evt) - assert mpoles_res is target_mpoles_view - if events: - mpoles.add_event(events[-1]) + if evt is not None: + mpoles.add_event(evt) - return (mpoles, SumpyTimingFuture(self.queue, events)) + return mpoles def eval_direct(self, target_boxes, source_box_starts, source_box_lists, src_weights): @@ -447,8 +394,6 @@ class SumpyExpansionWrangler(object): kwargs.update(self.box_source_list_kwargs()) kwargs.update(self.box_target_list_kwargs()) - events = [] - evt, pot_res = self.code.p2p()(self.queue, target_boxes=target_boxes, source_box_starts=source_box_starts, @@ -457,13 +402,12 @@ class SumpyExpansionWrangler(object): result=pot, **kwargs) - events.append(evt) for pot_i, pot_res_i in zip(pot, pot_res): assert pot_i is pot_res_i pot_i.add_event(evt) - return (pot, SumpyTimingFuture(self.queue, events)) + return pot def multipole_to_local(self, level_start_target_box_nrs, @@ -471,8 +415,6 @@ class SumpyExpansionWrangler(object): mpole_exps): local_exps = self.local_expansion_zeros() - events = [] - for lev in range(self.tree.nlevels): start, stop = level_start_target_box_nrs[lev:lev+2] if start == stop: @@ -503,9 +445,8 @@ class SumpyExpansionWrangler(object): tgt_rscale=level_to_rscale(self.tree, lev), **self.kernel_extra_kwargs) - events.append(evt) - return (local_exps, SumpyTimingFuture(self.queue, events)) + return local_exps def eval_multipoles(self, target_boxes_by_source_level, source_boxes_by_level, mpole_exps): @@ -514,10 +455,9 @@ class SumpyExpansionWrangler(object): kwargs = self.kernel_extra_kwargs.copy() kwargs.update(self.box_target_list_kwargs()) - events = [] - wait_for = mpole_exps.events + has_evt = False for isrc_level, ssn in enumerate(source_boxes_by_level): if len(target_boxes_by_source_level[isrc_level]) == 0: continue @@ -544,18 +484,19 @@ class SumpyExpansionWrangler(object): wait_for=wait_for, **kwargs) - events.append(evt) + has_evt = True wait_for = [evt] for pot_i, pot_res_i in zip(pot, pot_res): assert pot_i is pot_res_i - if events: + if has_evt: for pot_i in pot: - pot_i.add_event(events[-1]) + # Intentionally only adding the last event. + pot_i.add_event(evt) - return (pot, SumpyTimingFuture(self.queue, events)) + return pot def form_locals(self, level_start_target_or_target_parent_box_nrs, @@ -565,8 +506,6 @@ class SumpyExpansionWrangler(object): kwargs = self.extra_kwargs.copy() kwargs.update(self.box_source_list_kwargs()) - events = [] - for lev in range(self.tree.nlevels): start, stop = \ level_start_target_or_target_parent_box_nrs[lev:lev+2] @@ -592,19 +531,15 @@ class SumpyExpansionWrangler(object): rscale=level_to_rscale(self.tree, lev), **kwargs) - events.append(evt) assert result is target_local_exps_view - return (local_exps, SumpyTimingFuture(self.queue, events)) + return local_exps def refine_locals(self, level_start_target_or_target_parent_box_nrs, target_or_target_parent_boxes, local_exps): - - events = [] - for target_lev in range(1, self.tree.nlevels): start, stop = level_start_target_or_target_parent_box_nrs[ target_lev:target_lev+2] @@ -635,13 +570,12 @@ class SumpyExpansionWrangler(object): tgt_rscale=level_to_rscale(self.tree, target_lev), **self.kernel_extra_kwargs) - events.append(evt) assert local_exps_res is target_local_exps_view local_exps.add_event(evt) - return (local_exps, SumpyTimingFuture(self.queue, [evt])) + return local_exps def eval_locals(self, level_start_target_box_nrs, target_boxes, local_exps): pot = self.output_zeros() @@ -649,8 +583,6 @@ class SumpyExpansionWrangler(object): kwargs = self.kernel_extra_kwargs.copy() kwargs.update(self.box_target_list_kwargs()) - events = [] - for lev in range(self.tree.nlevels): start, stop = level_start_target_box_nrs[lev:lev+2] if start == stop: @@ -674,12 +606,11 @@ class SumpyExpansionWrangler(object): rscale=level_to_rscale(self.tree, lev), **kwargs) - events.append(evt) for pot_i, pot_res_i in zip(pot, pot_res): assert pot_i is pot_res_i - return (pot, SumpyTimingFuture(self.queue, events)) + return pot def finalize_potentials(self, potentials): return potentials diff --git a/test/test_fmm.py b/test/test_fmm.py index 71e3f044..0331db6c 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -234,62 +234,6 @@ def test_sumpy_fmm(ctx_getter, knl, local_expn_class, mpole_expn_class): pconv_verifier() -def test_sumpy_fmm_timing_data_collection(ctx_getter): - logging.basicConfig(level=logging.INFO) - - ctx = ctx_getter() - queue = cl.CommandQueue( - ctx, - properties=cl.command_queue_properties.PROFILING_ENABLE) - - nsources = 500 - dtype = np.float64 - - from boxtree.tools import ( - make_normal_particle_array as p_normal) - - knl = LaplaceKernel(2) - local_expn_class = VolumeTaylorLocalExpansion - mpole_expn_class = VolumeTaylorMultipoleExpansion - order = 1 - - sources = p_normal(queue, nsources, knl.dim, dtype, seed=15) - - from boxtree import TreeBuilder - tb = TreeBuilder(ctx) - - tree, _ = tb(queue, sources, - max_particles_in_box=30, debug=True) - - from boxtree.traversal import FMMTraversalBuilder - tbuild = FMMTraversalBuilder(ctx) - trav, _ = tbuild(queue, tree, debug=True) - - from pyopencl.clrandom import PhiloxGenerator - rng = PhiloxGenerator(ctx) - weights = rng.uniform(queue, nsources, dtype=np.float64) - - out_kernels = [knl] - - from functools import partial - - from sumpy.fmm import SumpyExpansionWranglerCodeContainer - wcc = SumpyExpansionWranglerCodeContainer( - ctx, - partial(mpole_expn_class, knl), - partial(local_expn_class, knl), - out_kernels) - - wrangler = wcc.get_wrangler(queue, tree, dtype, - fmm_level_to_order=lambda kernel, kernel_args, tree, lev: order) - from boxtree.fmm import drive_fmm - - timing_data = {} - pot, = drive_fmm(trav, wrangler, weights, timing_data=timing_data) - print(timing_data) - assert timing_data - - def test_sumpy_fmm_exclude_self(ctx_getter): logging.basicConfig(level=logging.INFO) -- GitLab From 83bcc8e24a0afaff1f18dd4a2e917fcf8ac64ccd Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 5 Aug 2018 16:03:29 +0800 Subject: [PATCH 13/20] Handle special cases for rscale --- sumpy/e2p.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index d0d05cbf..c3e8bb47 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -320,7 +320,12 @@ class E2PFromCSR(E2PBase): # meaningfully inferred. Make the type of rscale explicit. rscale = centers.dtype.type(kwargs.pop("rscale")) - return knl(queue, centers=centers, rscale=rscale, **kwargs) + try: + return knl(queue, centers=centers, rscale=rscale, **kwargs) + except TypeError as e: + logger.debug("Initial trial of p2e failed, probably the loopy kernel does not accept rscale") + logger.debug("Trying to make self-adjustments") + return knl(queue, centers=centers, **kwargs) # }}} -- GitLab From 430495c78ffe3b6b7b4e5426aa4115986ad84892 Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 5 Aug 2018 16:08:04 +0800 Subject: [PATCH 14/20] Handle special cases for rscale --- sumpy/e2p.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index c3e8bb47..85f29ec0 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -32,6 +32,8 @@ import sumpy.symbolic as sym from sumpy.tools import KernelCacheWrapper from loopy.version import MOST_RECENT_LANGUAGE_VERSION +import logging +logger = logging.getLogger(__name__) __doc__ = """ -- GitLab From c75a22d68135cc59ddb43cf61100ab9074680231 Mon Sep 17 00:00:00 2001 From: xywei Date: Tue, 9 Oct 2018 23:42:41 -0500 Subject: [PATCH 15/20] Fix merge mistakes --- sumpy/fmm.py | 55 +++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 39 insertions(+), 16 deletions(-) diff --git a/sumpy/fmm.py b/sumpy/fmm.py index 0f7af188..bd1c51d9 100644 --- a/sumpy/fmm.py +++ b/sumpy/fmm.py @@ -195,6 +195,7 @@ class SumpyTimingFuture(object): # }}} + # {{{ expansion wrangler class SumpyExpansionWrangler(object): @@ -225,6 +226,7 @@ class SumpyExpansionWrangler(object): self.code = code_container self.queue = queue self.tree = tree + self.issued_timing_data_warning = False self.dtype = dtype @@ -355,6 +357,8 @@ class SumpyExpansionWrangler(object): kwargs = self.extra_kwargs.copy() kwargs.update(self.box_source_list_kwargs()) + events = [] + for lev in range(self.tree.nlevels): p2m = self.code.p2m(self.level_orders[lev]) start, stop = level_start_source_box_nrs[lev:lev+2] @@ -375,10 +379,11 @@ class SumpyExpansionWrangler(object): rscale=level_to_rscale(self.tree, lev), **kwargs) + events.append(evt) assert mpoles_res is mpoles_view - return mpoles + return (mpoles, SumpyTimingFuture(self.queue, events)) def coarsen_multipoles(self, level_start_source_parent_box_nrs, @@ -386,7 +391,7 @@ class SumpyExpansionWrangler(object): mpoles): tree = self.tree - evt = None + events = [] # nlevels-1 is the last valid level index # nlevels-2 is the last valid level that could have children @@ -428,12 +433,14 @@ class SumpyExpansionWrangler(object): tgt_rscale=level_to_rscale(self.tree, target_level), **self.kernel_extra_kwargs) + events.append(evt) + assert mpoles_res is target_mpoles_view - if evt is not None: - mpoles.add_event(evt) + if events: + mpoles.add_event(events[-1]) - return mpoles + return (mpoles, SumpyTimingFuture(self.queue, events)) def eval_direct(self, target_boxes, source_box_starts, source_box_lists, src_weights): @@ -444,6 +451,8 @@ class SumpyExpansionWrangler(object): kwargs.update(self.box_source_list_kwargs()) kwargs.update(self.box_target_list_kwargs()) + events = [] + evt, pot_res = self.code.p2p()(self.queue, target_boxes=target_boxes, source_box_starts=source_box_starts, @@ -452,12 +461,13 @@ class SumpyExpansionWrangler(object): result=pot, **kwargs) + events.append(evt) for pot_i, pot_res_i in zip(pot, pot_res): assert pot_i is pot_res_i pot_i.add_event(evt) - return pot + return (pot, SumpyTimingFuture(self.queue, events)) def multipole_to_local(self, level_start_target_box_nrs, @@ -465,6 +475,8 @@ class SumpyExpansionWrangler(object): mpole_exps): local_exps = self.local_expansion_zeros() + events = [] + for lev in range(self.tree.nlevels): start, stop = level_start_target_box_nrs[lev:lev+2] if start == stop: @@ -495,8 +507,9 @@ class SumpyExpansionWrangler(object): tgt_rscale=level_to_rscale(self.tree, lev), **self.kernel_extra_kwargs) + events.append(evt) - return local_exps + return (local_exps, SumpyTimingFuture(self.queue, events)) def eval_multipoles(self, target_boxes_by_source_level, source_boxes_by_level, mpole_exps): @@ -505,9 +518,10 @@ class SumpyExpansionWrangler(object): kwargs = self.kernel_extra_kwargs.copy() kwargs.update(self.box_target_list_kwargs()) + events = [] + wait_for = mpole_exps.events - has_evt = False for isrc_level, ssn in enumerate(source_boxes_by_level): if len(target_boxes_by_source_level[isrc_level]) == 0: continue @@ -534,19 +548,18 @@ class SumpyExpansionWrangler(object): wait_for=wait_for, **kwargs) + events.append(evt) - has_evt = True wait_for = [evt] for pot_i, pot_res_i in zip(pot, pot_res): assert pot_i is pot_res_i - if has_evt: + if events: for pot_i in pot: - # Intentionally only adding the last event. - pot_i.add_event(evt) + pot_i.add_event(events[-1]) - return pot + return (pot, SumpyTimingFuture(self.queue, events)) def form_locals(self, level_start_target_or_target_parent_box_nrs, @@ -556,6 +569,8 @@ class SumpyExpansionWrangler(object): kwargs = self.extra_kwargs.copy() kwargs.update(self.box_source_list_kwargs()) + events = [] + for lev in range(self.tree.nlevels): start, stop = \ level_start_target_or_target_parent_box_nrs[lev:lev+2] @@ -581,15 +596,19 @@ class SumpyExpansionWrangler(object): rscale=level_to_rscale(self.tree, lev), **kwargs) + events.append(evt) assert result is target_local_exps_view - return local_exps + return (local_exps, SumpyTimingFuture(self.queue, events)) def refine_locals(self, level_start_target_or_target_parent_box_nrs, target_or_target_parent_boxes, local_exps): + + events = [] + for target_lev in range(1, self.tree.nlevels): start, stop = level_start_target_or_target_parent_box_nrs[ target_lev:target_lev+2] @@ -620,12 +639,13 @@ class SumpyExpansionWrangler(object): tgt_rscale=level_to_rscale(self.tree, target_lev), **self.kernel_extra_kwargs) + events.append(evt) assert local_exps_res is target_local_exps_view local_exps.add_event(evt) - return local_exps + return (local_exps, SumpyTimingFuture(self.queue, [evt])) def eval_locals(self, level_start_target_box_nrs, target_boxes, local_exps): pot = self.output_zeros() @@ -633,6 +653,8 @@ class SumpyExpansionWrangler(object): kwargs = self.kernel_extra_kwargs.copy() kwargs.update(self.box_target_list_kwargs()) + events = [] + for lev in range(self.tree.nlevels): start, stop = level_start_target_box_nrs[lev:lev+2] if start == stop: @@ -656,11 +678,12 @@ class SumpyExpansionWrangler(object): rscale=level_to_rscale(self.tree, lev), **kwargs) + events.append(evt) for pot_i, pot_res_i in zip(pot, pot_res): assert pot_i is pot_res_i - return pot + return (pot, SumpyTimingFuture(self.queue, events)) def finalize_potentials(self, potentials): return potentials -- GitLab From 8190c9399c9be3351de0c90c2244262f7a6668e7 Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 26 Jan 2019 22:18:55 -0600 Subject: [PATCH 16/20] Refine changes --- setup.py | 2 +- test/test_fmm.py | 56 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 57 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 8d085054..9457bbdc 100644 --- a/setup.py +++ b/setup.py @@ -94,7 +94,7 @@ setup(name="sumpy", install_requires=[ "pytools>=2018.2", "loo.py>=2017.2", - "boxtree>=2013.1", + "boxtree>=2018.1", "pytest>=2.3", "six", diff --git a/test/test_fmm.py b/test/test_fmm.py index 0331db6c..71e3f044 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -234,6 +234,62 @@ def test_sumpy_fmm(ctx_getter, knl, local_expn_class, mpole_expn_class): pconv_verifier() +def test_sumpy_fmm_timing_data_collection(ctx_getter): + logging.basicConfig(level=logging.INFO) + + ctx = ctx_getter() + queue = cl.CommandQueue( + ctx, + properties=cl.command_queue_properties.PROFILING_ENABLE) + + nsources = 500 + dtype = np.float64 + + from boxtree.tools import ( + make_normal_particle_array as p_normal) + + knl = LaplaceKernel(2) + local_expn_class = VolumeTaylorLocalExpansion + mpole_expn_class = VolumeTaylorMultipoleExpansion + order = 1 + + sources = p_normal(queue, nsources, knl.dim, dtype, seed=15) + + from boxtree import TreeBuilder + tb = TreeBuilder(ctx) + + tree, _ = tb(queue, sources, + max_particles_in_box=30, debug=True) + + from boxtree.traversal import FMMTraversalBuilder + tbuild = FMMTraversalBuilder(ctx) + trav, _ = tbuild(queue, tree, debug=True) + + from pyopencl.clrandom import PhiloxGenerator + rng = PhiloxGenerator(ctx) + weights = rng.uniform(queue, nsources, dtype=np.float64) + + out_kernels = [knl] + + from functools import partial + + from sumpy.fmm import SumpyExpansionWranglerCodeContainer + wcc = SumpyExpansionWranglerCodeContainer( + ctx, + partial(mpole_expn_class, knl), + partial(local_expn_class, knl), + out_kernels) + + wrangler = wcc.get_wrangler(queue, tree, dtype, + fmm_level_to_order=lambda kernel, kernel_args, tree, lev: order) + from boxtree.fmm import drive_fmm + + timing_data = {} + pot, = drive_fmm(trav, wrangler, weights, timing_data=timing_data) + print(timing_data) + assert timing_data + + def test_sumpy_fmm_exclude_self(ctx_getter): logging.basicConfig(level=logging.INFO) -- GitLab From 79ac8bcb340f89b5db83ffa35f5ea8297f678204 Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 26 Jan 2019 22:33:43 -0600 Subject: [PATCH 17/20] Flake8 fixes --- sumpy/e2p.py | 5 +++-- sumpy/kernel.py | 4 ++-- sumpy/p2e.py | 5 +++-- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index 85f29ec0..1558ad7c 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -324,8 +324,9 @@ class E2PFromCSR(E2PBase): try: return knl(queue, centers=centers, rscale=rscale, **kwargs) - except TypeError as e: - logger.debug("Initial trial of p2e failed, probably the loopy kernel does not accept rscale") + except TypeError: + logger.debug("Initial trial of p2e failed, " + "probably the loopy kernel does not accept rscale") logger.debug("Trying to make self-adjustments") return knl(queue, centers=centers, **kwargs) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 48aed8ce..6fc1a669 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -732,7 +732,7 @@ class FactorizedBiharmonicKernel(ExpressionKernel): # http://dlmf.nist.gov/10.27#E8 expr_for_K0 = [var("hankel_1")(0, var("I")*lam[i]*r) - for i in range(2)] + for i in range(2)] # noqa: N806 scaling_for_K0 = 1/2*var("pi")*var("I") # noqa: N806 expr = expr_for_K0[0] - expr_for_K0[1] @@ -755,7 +755,7 @@ class FactorizedBiharmonicKernel(ExpressionKernel): def update_persistent_hash(self, key_hash, key_builder): key_hash.update(type(self).__name__.encode("utf8")) - key_builder.rec(key_hash, + key_builder.rec(key_hash, (self.dim, self.lambda1_name, self.lambda2_name)) def __repr__(self): diff --git a/sumpy/p2e.py b/sumpy/p2e.py index 99a870c1..db8f491e 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -301,8 +301,9 @@ class P2EFromCSR(P2EBase): try: return knl(queue, centers=centers, rscale=rscale, **kwargs) - except TypeError as e: - logger.debug("Initial trial of p2e failed, probably the loopy kernel does not accept rscale") + except TypeError: + logger.debug("Initial trial of p2e failed, " + "probably the loopy kernel does not accept rscale") logger.debug("Trying to make self-adjustments") return knl(queue, centers=centers, **kwargs) -- GitLab From dcc048952b01b5662c8a16c1548d75bb530d616e Mon Sep 17 00:00:00 2001 From: xywei Date: Sat, 26 Jan 2019 22:35:18 -0600 Subject: [PATCH 18/20] More Flake8 fixes --- sumpy/kernel.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 6fc1a669..5a47ed56 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -731,8 +731,8 @@ class FactorizedBiharmonicKernel(ExpressionKernel): r = pymbolic_real_norm_2(make_sym_vector("d", dim)) # http://dlmf.nist.gov/10.27#E8 - expr_for_K0 = [var("hankel_1")(0, var("I")*lam[i]*r) - for i in range(2)] # noqa: N806 + expr_for_K0 = [var("hankel_1")(0, var("I")*lam[i]*r) # noqa: N806 + for i in range(2)] scaling_for_K0 = 1/2*var("pi")*var("I") # noqa: N806 expr = expr_for_K0[0] - expr_for_K0[1] -- GitLab From dea7dbf134aaa9046e5be6da8e4a1bc33e6c7032 Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 17 Feb 2019 21:53:37 -0600 Subject: [PATCH 19/20] Rename d(lap(.)) class --- sumpy/kernel.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 5a47ed56..2c1215aa 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -989,7 +989,7 @@ class DirectionalSourceDerivative(DirectionalDerivative): mapper_method = "map_directional_source_derivative" -class LaplacianDerivative(DerivativeBase): +class LaplacianDerivativeBase(DerivativeBase): init_arg_names = ("inner_kernel", ) def __init__(self, inner_kernel): @@ -1015,7 +1015,7 @@ class LaplacianDerivative(DerivativeBase): return self.inner_kernel.get_source_args() -class LaplacianTargetDerivative(LaplacianDerivative): +class LaplacianTargetDerivative(LaplacianDerivativeBase): laplacian_kind = "tgt" def get_code_transformer(self): @@ -1039,7 +1039,7 @@ class LaplacianTargetDerivative(LaplacianDerivative): mapper_method = "map_laplacian_target_derivative" -class LaplacianSourceDerivative(LaplacianDerivative): +class LaplacianSourceDerivative(LaplacianDerivativeBase): laplacian_kind = "src" def get_code_transformer(self): -- GitLab From 4e4bf75e9dfde68afbab3c5546381c4df5adcab0 Mon Sep 17 00:00:00 2001 From: xywei Date: Sun, 17 Feb 2019 21:56:52 -0600 Subject: [PATCH 20/20] Link new mappers in docstring --- sumpy/kernel.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 2c1215aa..8f45e2e4 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -64,6 +64,9 @@ of them in the process. .. autoclass:: AxisTargetDerivative .. autoclass:: DirectionalTargetDerivative .. autoclass:: DirectionalSourceDerivative +.. autoclass:: LaplacianDerivativeBase +.. autoclass:: LaplacianTargetDerivative +.. autoclass:: LaplacianSourceDerivative Transforming kernels -------------------- -- GitLab