From b67ab80e45e8e20626c7899edaed0084b1c963c4 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 20 May 2017 18:04:52 -0500 Subject: [PATCH 1/8] Mini progress towards expansion toys [ci skip] --- sumpy/toys.py | 70 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 sumpy/toys.py diff --git a/sumpy/toys.py b/sumpy/toys.py new file mode 100644 index 00000000..ac0cf39b --- /dev/null +++ b/sumpy/toys.py @@ -0,0 +1,70 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2017 Andreas Kloeckner" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import six +from six.moves import range + +import numpy as np +import loopy as lp + +import logging +logger = logging.getLogger(__name__) + + +class ToyContext(object): + def __init__(self, cl_context, kernel): + self.code_cache = {} + + +class PotentialSource(object): + def __init__(self, toy_ctx, center): + self.toy_ctx = toy_ctx + self.center + + def eval(self, toy_ctx, targets): + raise NotImplementedError() + + def plot(self, extent): + pass + + +class PointSources(PotentialSource): + """ + .. attribute:: points + + ``[ndim, npoints]`` + """ + + def __init__(self, toy_ctx, points, weights, center=None): + super(PointSources, self).__init__(toy_ctx, center) + + +class Multipole(PotentialSource): + def __init__(self, toy_ctx, center, coeffs): + super(Multipole, self).__init__(toy_ctx) + + +class Local(PotentialSource): + def __init__(self, toy_ctx, center, coeffs): + super(Local, self).__init__(toy_ctx) -- GitLab From aeeeb20490d4d29a4e37ecdb8c50eea7efd3047d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 21 May 2017 18:07:34 -0700 Subject: [PATCH 2/8] More work on expansion toy --- sumpy/toys.py | 209 +++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 190 insertions(+), 19 deletions(-) diff --git a/sumpy/toys.py b/sumpy/toys.py index ac0cf39b..7333f53a 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -22,31 +22,138 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -import six -from six.moves import range +import six # noqa: F401 +from six.moves import range # noqa: F401 -import numpy as np -import loopy as lp +from pytools import memoize_method + +import numpy as np # noqa: F401 +import loopy as lp # noqa: F401 +import pyopencl as cl import logging logger = logging.getLogger(__name__) +# {{{ context + class ToyContext(object): - def __init__(self, cl_context, kernel): - self.code_cache = {} + def __init__(self, cl_context, kernel, + mpole_expn_class=None, + local_expn_class=None, + extra_source_kwargs=None): + self.cl_context = cl_context + self.queue = cl.CommandQueue(self.cl_context) + self.kernel = kernel + + if mpole_expn_class is None: + from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion + mpole_expn_class = VolumeTaylorMultipoleExpansion + if local_expn_class is None: + from sumpy.expansion.local import VolumeTaylorLocalExpansion + local_expn_class = VolumeTaylorLocalExpansion + + if extra_source_kwargs is None: + extra_source_kwargs = {} + + self.mpole_expn_class = mpole_expn_class + self.local_expn_class = local_expn_class + self.extra_source_kwargs = extra_source_kwargs + + @memoize_method + def get_p2p(self): + from sumpy.p2p import P2P + return P2P(self.cl_context, [self.kernel], exclude_self=False) + + @memoize_method + def get_p2m(self, order): + from sumpy import P2EFromSingleBox + return P2EFromSingleBox(self.cl_context, + self.mpole_expn_class(self.kernel, order), + [self.kernel]) + + @memoize_method + def get_p2l(self, order): + from sumpy import P2EFromSingleBox + return P2EFromSingleBox(self.cl_context, + self.local_expn_class(self.kernel, order), + [self.kernel]) + +# }}} + + +# {{{ helpers + +def _p2e(psource, center, order, p2e, expn_class): + source_boxes = np.array([0], dtype=np.int32) + box_source_starts = np.array([0], dtype=np.int32) + box_source_counts_nonchild = np.array( + [psource.points.shape[-1]], dtype=np.int32) + + toy_ctx = psource.toy_ctx + center = np.asarray(center) + centers = np.array(center, dtype=np.float64).reshape( + toy_ctx.kernel.dim, 1) + + evt, (coeffs,) = p2e(toy_ctx.queue, + source_boxes=source_boxes, + box_source_starts=box_source_starts, + box_source_counts_nonchild=box_source_counts_nonchild, + centers=centers, + sources=psource.points, + strengths=psource.weights, + nboxes=1, + tgt_base_ibox=0, + + #flags="print_hl_cl", + out_host=True, **toy_ctx.extra_source_kwargs) + + return expn_class(toy_ctx, center, order, coeffs[0]) + + +def _e2p(psource, targets, e2p): + ntargets = targets.shape[-1] + + _boxes = np.array([0], dtype=np.int32) + + box_target_starts = np.array([0], dtype=np.int32) + box_target_counts_nonchild = np.array([ntargets], dtype=np.int32) + + coeffs = np.array([psource.coeffs]) + toy_ctx = psource.toy_ctx + evt, (pot,) = e2p( + toy_ctx.queue, + src_expansions=coeffs, + src_base_ibox=0, + target_boxes=source_boxes, + box_target_starts=box_target_starts, + box_target_counts_nonchild=box_target_counts_nonchild, + centers=centers, + targets=targets, + #flags="print_hl_cl", + out_host=True, **extra_kwargs) + +# }}} class PotentialSource(object): - def __init__(self, toy_ctx, center): + def __init__(self, toy_ctx): self.toy_ctx = toy_ctx - self.center - def eval(self, toy_ctx, targets): + def eval(self, targets): raise NotImplementedError() - def plot(self, extent): - pass + def __add__(self, other): + if not isinstance(other, PointSources): + raise TypeError() + + return LinearCombination((1, 1), (self, other)) + + def __sub__(self, other): + if not isinstance(other, PointSources): + raise TypeError() + + return LinearCombination((1, -1), (self, other)) class PointSources(PotentialSource): @@ -56,15 +163,79 @@ class PointSources(PotentialSource): ``[ndim, npoints]`` """ - def __init__(self, toy_ctx, points, weights, center=None): - super(PointSources, self).__init__(toy_ctx, center) + def __init__(self, toy_ctx, points, weights): + super(PointSources, self).__init__(toy_ctx) + + self.points = points + self.weights = weights + + def eval(self, targets): + evt, (potential,) = self.toy_ctx.get_p2p()( + self.toy_ctx.queue, targets, self.points, [self.weights], + out_host=True) + + return potential + + +class ExpansionPotentialSource(PotentialSource): + def __init__(self, toy_ctx, center, order, coeffs): + super(ExpansionPotentialSource, self).__init__(toy_ctx) + + +class MultipoleExpansion(ExpansionPotentialSource): + pass + + +class LocalExpansion(ExpansionPotentialSource): + def eval(self, targets): + raise NotImplementedError() + + +def multipole_expand(psource, center, order): + if isinstance(psource, PointSources): + return _p2e(psource, center, order, psource.toy_ctx.get_p2m(order), + MultipoleExpansion) + + elif isinstance(psource, MultipoleExpansion): + raise NotImplementedError() + else: + raise TypeError("do not know how to expand '%s'" + % type(psource).__name__) + + +def local_expand(psource, center, order): + if isinstance(psource, PointSources): + return _p2e(psource, center, order, psource.toy_ctx.get_p2l(order), + LocalExpansion) + + elif isinstance(psource, MultipoleExpansion): + raise NotImplementedError() + elif isinstance(psource, LocalExpansion): + raise NotImplementedError() + else: + raise TypeError("do not know how to expand '%s'" + % type(psource).__name__) + + +class LinearCombination(PotentialSource): + def __init__(self, coeffs, psources): + from pytools import single_valued + super(LinearCombination, self).__init__( + single_valued(psource.toy_ctx for psource in psources)) + + self.coeffs = coeffs + self.psources = psources + + def eval(self, targets): + result = 0 + for coeff, psource in zip(self.coeffs, self.psources): + result += coeff * psource.eval(targets) + return result -class Multipole(PotentialSource): - def __init__(self, toy_ctx, center, coeffs): - super(Multipole, self).__init__(toy_ctx) +def logplot(fp, psource, **kwargs): + fp.show_scalar_in_matplotlib( + np.log10(np.abs(psource.eval(fp.points) + 1e-15)), **kwargs) -class Local(PotentialSource): - def __init__(self, toy_ctx, center, coeffs): - super(Local, self).__init__(toy_ctx) +# vim: foldmethod=marker -- GitLab From 36aa84cd366b5618eba3b0886d5d0b27e3d84aed Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 21 May 2017 18:08:07 -0700 Subject: [PATCH 3/8] Add some sample code for expansion toy --- examples/expansion-toys.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 examples/expansion-toys.py diff --git a/examples/expansion-toys.py b/examples/expansion-toys.py new file mode 100644 index 00000000..0e8cbd49 --- /dev/null +++ b/examples/expansion-toys.py @@ -0,0 +1,28 @@ +import pyopencl as cl +import sumpy.toys as t +import numpy as np +from sumpy.visualization import FieldPlotter +import matplotlib.pyplot as plt + + +def main(): + from sumpy.kernel import LaplaceKernel + tctx = t.ToyContext(cl.create_some_context(), LaplaceKernel(2)) + + pt_src = t.PointSources( + tctx, + np.random.rand(2, 50) - 0.5, + np.ones(50)) + + fp = FieldPlotter([0, 0], extent=4) + + if 0: + t.logplot(fp, pt_src, cmap="jet") + plt.colorbar() + plt.show() + + t.local_expand(pt_src, [3, 0], 5) + + +if __name__ == "__main__": + main() -- GitLab From 2bdd64386b8dcdfcce2d2f5ff6b7e3a6dcf533bf Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 25 May 2017 22:41:16 -0700 Subject: [PATCH 4/8] Expansion toybox: Add expressions, M/L eval --- examples/expansion-toys.py | 14 ++++- sumpy/toys.py | 109 +++++++++++++++++++++++++++++++------ 2 files changed, 103 insertions(+), 20 deletions(-) diff --git a/examples/expansion-toys.py b/examples/expansion-toys.py index 0e8cbd49..5c0714b9 100644 --- a/examples/expansion-toys.py +++ b/examples/expansion-toys.py @@ -14,14 +14,24 @@ def main(): np.random.rand(2, 50) - 0.5, np.ones(50)) - fp = FieldPlotter([0, 0], extent=4) + fp = FieldPlotter([3, 0], extent=8) if 0: t.logplot(fp, pt_src, cmap="jet") plt.colorbar() plt.show() - t.local_expand(pt_src, [3, 0], 5) + lexp = t.multipole_expand(pt_src, [0, 0], 5) + + diff = lexp - pt_src + + if 1: + t.logplot(fp, diff, cmap="jet", vmin=-3) + plt.colorbar() + plt.show() + + + if __name__ == "__main__": diff --git a/sumpy/toys.py b/sumpy/toys.py index 7333f53a..6b049973 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -26,6 +26,7 @@ import six # noqa: F401 from six.moves import range # noqa: F401 from pytools import memoize_method +from numbers import Number import numpy as np # noqa: F401 import loopy as lp # noqa: F401 @@ -79,6 +80,20 @@ class ToyContext(object): self.local_expn_class(self.kernel, order), [self.kernel]) + @memoize_method + def get_m2p(self, order): + from sumpy import E2PFromSingleBox + return E2PFromSingleBox(self.cl_context, + self.mpole_expn_class(self.kernel, order), + [self.kernel]) + + @memoize_method + def get_l2p(self, order): + from sumpy import E2PFromSingleBox + return E2PFromSingleBox(self.cl_context, + self.local_expn_class(self.kernel, order), + [self.kernel]) + # }}} @@ -114,24 +129,29 @@ def _p2e(psource, center, order, p2e, expn_class): def _e2p(psource, targets, e2p): ntargets = targets.shape[-1] - _boxes = np.array([0], dtype=np.int32) + boxes = np.array([0], dtype=np.int32) box_target_starts = np.array([0], dtype=np.int32) box_target_counts_nonchild = np.array([ntargets], dtype=np.int32) - coeffs = np.array([psource.coeffs]) toy_ctx = psource.toy_ctx + centers = np.array(psource.center, dtype=np.float64).reshape( + toy_ctx.kernel.dim, 1) + + coeffs = np.array([psource.coeffs]) evt, (pot,) = e2p( toy_ctx.queue, src_expansions=coeffs, src_base_ibox=0, - target_boxes=source_boxes, + target_boxes=boxes, box_target_starts=box_target_starts, box_target_counts_nonchild=box_target_counts_nonchild, centers=centers, targets=targets, #flags="print_hl_cl", - out_host=True, **extra_kwargs) + out_host=True, **toy_ctx.extra_source_kwargs) + + return pot # }}} @@ -143,17 +163,56 @@ class PotentialSource(object): def eval(self, targets): raise NotImplementedError() + def __neg__(self): + return -1*self + def __add__(self, other): - if not isinstance(other, PointSources): - raise TypeError() + if isinstance(other, (Number, np.number)): + other = ConstantPotential(self.toy_ctx, other) + elif not isinstance(other, PotentialSource): + return NotImplemented - return LinearCombination((1, 1), (self, other)) + return Sum((self, other)) + + __radd__ = __add__ def __sub__(self, other): - if not isinstance(other, PointSources): - raise TypeError() + return self.__add__(-other) + + def __rsub__(self, other): + return (-self).__add__(other) + + def __mul__(self, other): + if isinstance(other, (Number, np.number)): + other = ConstantPotential(self.toy_ctx, other) + elif not isinstance(other, PotentialSource): + return NotImplemented + + return Product((self, other)) - return LinearCombination((1, -1), (self, other)) + __rmul__ = __mul__ + + +class ConstantPotential(PotentialSource): + def __init__(self, toy_ctx, value): + super(ConstantPotential, self).__init__(toy_ctx) + self.value = np.array(value) + + def eval(self, targets): + pot = np.empty(targets.shape[-1], dtype=self.value.dtype) + pot.fill(self.value) + return pot + + +class OneOnBallPotential(PotentialSource): + def __init__(self, toy_ctx, center, radius): + super(PointSources, self).__init__(toy_ctx) + self.center = center + self.radius = radius + + def eval(self, targets): + dist_vec = targets - self.center[:, np.newaxis] + return (np.sum(dist_vec**2, axis=0) < self.radius).astype(np.float64) class PointSources(PotentialSource): @@ -180,15 +239,19 @@ class PointSources(PotentialSource): class ExpansionPotentialSource(PotentialSource): def __init__(self, toy_ctx, center, order, coeffs): super(ExpansionPotentialSource, self).__init__(toy_ctx) + self.center = center + self.order = order + self.coeffs = coeffs class MultipoleExpansion(ExpansionPotentialSource): - pass + def eval(self, targets): + return _e2p(self, targets, self.toy_ctx.get_m2p(self.order)) class LocalExpansion(ExpansionPotentialSource): def eval(self, targets): - raise NotImplementedError() + return _e2p(self, targets, self.toy_ctx.get_l2p(self.order)) def multipole_expand(psource, center, order): @@ -217,19 +280,29 @@ def local_expand(psource, center, order): % type(psource).__name__) -class LinearCombination(PotentialSource): - def __init__(self, coeffs, psources): +class PotentialExpressionNode(PotentialSource): + def __init__(self, psources): from pytools import single_valued - super(LinearCombination, self).__init__( + super(PotentialExpressionNode, self).__init__( single_valued(psource.toy_ctx for psource in psources)) - self.coeffs = coeffs self.psources = psources + +class Sum(PotentialExpressionNode): def eval(self, targets): result = 0 - for coeff, psource in zip(self.coeffs, self.psources): - result += coeff * psource.eval(targets) + for psource in self.psources: + result = result + psource.eval(targets) + + return result + + +class Product(PotentialExpressionNode): + def eval(self, targets): + result = 1 + for psource in self.psources: + result = result * psource.eval(targets) return result -- GitLab From 53b6d98022bef3a0e2608d15507e97567f53432f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 08:23:18 -0700 Subject: [PATCH 5/8] Complete the set of available translations in the toy --- examples/expansion-toys.py | 7 +- sumpy/toys.py | 138 ++++++++++++++++++++++++++++++------- 2 files changed, 118 insertions(+), 27 deletions(-) diff --git a/examples/expansion-toys.py b/examples/expansion-toys.py index 5c0714b9..72be2c7e 100644 --- a/examples/expansion-toys.py +++ b/examples/expansion-toys.py @@ -21,8 +21,13 @@ def main(): plt.colorbar() plt.show() - lexp = t.multipole_expand(pt_src, [0, 0], 5) + mexp = t.multipole_expand(pt_src, [0, 0], 9) + mexp2 = t.multipole_expand(mexp, [0, 0.25]) + lexp = t.local_expand(mexp, [3, 0]) + lexp2 = t.local_expand(lexp, [3, 1]) + diff = mexp - pt_src + diff = mexp2 - pt_src diff = lexp - pt_src if 1: diff --git a/sumpy/toys.py b/sumpy/toys.py index 6b049973..f61a9090 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -94,6 +94,27 @@ class ToyContext(object): self.local_expn_class(self.kernel, order), [self.kernel]) + @memoize_method + def get_m2m(self, from_order, to_order): + from sumpy import E2EFromCSR + return E2EFromCSR(self.cl_context, + self.mpole_expn_class(self.kernel, from_order), + self.mpole_expn_class(self.kernel, to_order)) + + @memoize_method + def get_m2l(self, from_order, to_order): + from sumpy import E2EFromCSR + return E2EFromCSR(self.cl_context, + self.mpole_expn_class(self.kernel, from_order), + self.local_expn_class(self.kernel, to_order)) + + @memoize_method + def get_l2l(self, from_order, to_order): + from sumpy import E2EFromCSR + return E2EFromCSR(self.cl_context, + self.local_expn_class(self.kernel, from_order), + self.local_expn_class(self.kernel, to_order)) + # }}} @@ -153,9 +174,48 @@ def _e2p(psource, targets, e2p): return pot + +def _e2e(psource, to_center, to_order, e2e, expn_class): + toy_ctx = psource.toy_ctx + + target_boxes = np.array([1], dtype=np.int32) + src_box_starts = np.array([0, 1], dtype=np.int32) + src_box_lists = np.array([0], dtype=np.int32) + + centers = (np.array( + [ + # box 0: source + psource.center, + + # box 1: target + to_center, + ], + dtype=np.float64)).T.copy() + + coeffs = np.array([psource.coeffs]) + + evt, (to_coeffs,) = e2e( + toy_ctx.queue, + src_expansions=coeffs, + src_base_ibox=0, + tgt_base_ibox=0, + ntgt_level_boxes=2, + + target_boxes=target_boxes, + + src_box_starts=src_box_starts, + src_box_lists=src_box_lists, + centers=centers, + #flags="print_hl_cl", + out_host=True, **toy_ctx.extra_source_kwargs) + + return expn_class(toy_ctx, to_center, to_order, to_coeffs[1]) + # }}} +# {{{ potential source classes + class PotentialSource(object): def __init__(self, toy_ctx): self.toy_ctx = toy_ctx @@ -254,32 +314,6 @@ class LocalExpansion(ExpansionPotentialSource): return _e2p(self, targets, self.toy_ctx.get_l2p(self.order)) -def multipole_expand(psource, center, order): - if isinstance(psource, PointSources): - return _p2e(psource, center, order, psource.toy_ctx.get_p2m(order), - MultipoleExpansion) - - elif isinstance(psource, MultipoleExpansion): - raise NotImplementedError() - else: - raise TypeError("do not know how to expand '%s'" - % type(psource).__name__) - - -def local_expand(psource, center, order): - if isinstance(psource, PointSources): - return _p2e(psource, center, order, psource.toy_ctx.get_p2l(order), - LocalExpansion) - - elif isinstance(psource, MultipoleExpansion): - raise NotImplementedError() - elif isinstance(psource, LocalExpansion): - raise NotImplementedError() - else: - raise TypeError("do not know how to expand '%s'" - % type(psource).__name__) - - class PotentialExpressionNode(PotentialSource): def __init__(self, psources): from pytools import single_valued @@ -306,6 +340,58 @@ class Product(PotentialExpressionNode): return result +# }}} + + +def multipole_expand(psource, center, order=None): + if isinstance(psource, PointSources): + if order is None: + raise ValueError("order may not be None") + + return _p2e(psource, center, order, psource.toy_ctx.get_p2m(order), + MultipoleExpansion) + + elif isinstance(psource, MultipoleExpansion): + if order is None: + order = psource.order + + return _e2e(psource, center, order, + psource.toy_ctx.get_m2m(psource.order, order), + MultipoleExpansion) + + else: + raise TypeError("do not know how to expand '%s'" + % type(psource).__name__) + + +def local_expand(psource, center, order=None): + if isinstance(psource, PointSources): + if order is None: + raise ValueError("order may not be None") + + return _p2e(psource, center, order, psource.toy_ctx.get_p2l(order), + LocalExpansion) + + elif isinstance(psource, MultipoleExpansion): + if order is None: + order = psource.order + + return _e2e(psource, center, order, + psource.toy_ctx.get_m2l(psource.order, order), + LocalExpansion) + + elif isinstance(psource, LocalExpansion): + if order is None: + order = psource.order + + return _e2e(psource, center, order, + psource.toy_ctx.get_l2l(psource.order, order), + LocalExpansion) + + else: + raise TypeError("do not know how to expand '%s'" + % type(psource).__name__) + def logplot(fp, psource, **kwargs): fp.show_scalar_in_matplotlib( -- GitLab From c4442e2c39edc180b6a04a564561a82d36e4a997 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 10:29:32 -0700 Subject: [PATCH 6/8] Expansion toys: restrict_inner, restrict_outer, l_inf --- examples/expansion-toys.py | 8 +++----- sumpy/toys.py | 32 +++++++++++++++++++++++++++++--- 2 files changed, 32 insertions(+), 8 deletions(-) diff --git a/examples/expansion-toys.py b/examples/expansion-toys.py index 72be2c7e..782bbcf5 100644 --- a/examples/expansion-toys.py +++ b/examples/expansion-toys.py @@ -28,16 +28,14 @@ def main(): diff = mexp - pt_src diff = mexp2 - pt_src - diff = lexp - pt_src + diff = lexp2 - pt_src + print(t.l_inf(diff, 1.2, center=lexp2.center)) if 1: - t.logplot(fp, diff, cmap="jet", vmin=-3) + t.logplot(fp, diff, cmap="jet", vmin=-3, vmax=0) plt.colorbar() plt.show() - - - if __name__ == "__main__": main() diff --git a/sumpy/toys.py b/sumpy/toys.py index f61a9090..64431538 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -266,8 +266,8 @@ class ConstantPotential(PotentialSource): class OneOnBallPotential(PotentialSource): def __init__(self, toy_ctx, center, radius): - super(PointSources, self).__init__(toy_ctx) - self.center = center + super(OneOnBallPotential, self).__init__(toy_ctx) + self.center = np.asarray(center) self.radius = radius def eval(self, targets): @@ -299,7 +299,7 @@ class PointSources(PotentialSource): class ExpansionPotentialSource(PotentialSource): def __init__(self, toy_ctx, center, order, coeffs): super(ExpansionPotentialSource, self).__init__(toy_ctx) - self.center = center + self.center = np.asarray(center) self.order = order self.coeffs = coeffs @@ -397,4 +397,30 @@ def logplot(fp, psource, **kwargs): fp.show_scalar_in_matplotlib( np.log10(np.abs(psource.eval(fp.points) + 1e-15)), **kwargs) + +def restrict_inner(psource, radius, center=None): + if center is None: + center = psource.center + + return psource * OneOnBallPotential(psource.toy_ctx, center, radius) + + +def restrict_outer(psource, radius, center=None): + if center is None: + center = psource.center + + return psource * (1-OneOnBallPotential(psource.toy_ctx, center, radius)) + + +def l_inf(psource, radius, center=None, npoints=100): + if center is None: + center = psource.center + + restr = psource * OneOnBallPotential(psource.toy_ctx, center, radius) + + from sumpy.visualization import FieldPlotter + fp = FieldPlotter(center, extent=2*radius, npoints=npoints) + return np.max(np.abs(restr.eval(fp.points))) + + # vim: foldmethod=marker -- GitLab From f26c9834a647dd2d82768954824fecbe0ec74531 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 12:06:25 -0700 Subject: [PATCH 7/8] Toys: Fix OneOnBallPotential --- sumpy/toys.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/toys.py b/sumpy/toys.py index 64431538..512dc3a1 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -272,7 +272,7 @@ class OneOnBallPotential(PotentialSource): def eval(self, targets): dist_vec = targets - self.center[:, np.newaxis] - return (np.sum(dist_vec**2, axis=0) < self.radius).astype(np.float64) + return (np.sum(dist_vec**2, axis=0) < self.radius**2).astype(np.float64) class PointSources(PotentialSource): -- GitLab From bf879de9a788ed2bf7818d66e7309dac7758c4a6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 26 May 2017 12:06:50 -0700 Subject: [PATCH 8/8] Toys: l_inf debug, draw_box, draw_circle --- sumpy/toys.py | 43 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 41 insertions(+), 2 deletions(-) diff --git a/sumpy/toys.py b/sumpy/toys.py index 512dc3a1..4641d697 100644 --- a/sumpy/toys.py +++ b/sumpy/toys.py @@ -412,7 +412,7 @@ def restrict_outer(psource, radius, center=None): return psource * (1-OneOnBallPotential(psource.toy_ctx, center, radius)) -def l_inf(psource, radius, center=None, npoints=100): +def l_inf(psource, radius, center=None, npoints=100, debug=False): if center is None: center = psource.center @@ -420,7 +420,46 @@ def l_inf(psource, radius, center=None, npoints=100): from sumpy.visualization import FieldPlotter fp = FieldPlotter(center, extent=2*radius, npoints=npoints) - return np.max(np.abs(restr.eval(fp.points))) + z = restr.eval(fp.points) + if debug: + fp.show_scalar_in_matplotlib( + np.log10(np.abs(z + 1e-15))) + import matplotlib.pyplot as pt + pt.colorbar() + pt.show() + + return np.max(np.abs(z)) + + +def draw_box(center, radius, **kwargs): + center = np.asarray(center) + + el = center - radius + eh = center + radius + + import matplotlib.pyplot as pt + import matplotlib.patches as mpatches + from matplotlib.path import Path + + pathdata = [ + (Path.MOVETO, (el[0], el[1])), + (Path.LINETO, (eh[0], el[1])), + (Path.LINETO, (eh[0], eh[1])), + (Path.LINETO, (el[0], eh[1])), + (Path.CLOSEPOLY, (el[0], el[1])), + ] + + codes, verts = zip(*pathdata) + path = Path(verts, codes) + patch = mpatches.PathPatch(path, **kwargs) + pt.gca().add_patch(patch) + + +def draw_circle(center, radius, **kwargs): + center = np.asarray(center) + + import matplotlib.pyplot as plt + plt.gca().add_patch(plt.Circle((center[0], center[1]), radius, **kwargs)) # vim: foldmethod=marker -- GitLab