From 62f9126d34e2c6c025734cb343f239a5c437a589 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 23 Dec 2016 15:38:00 -0600 Subject: [PATCH 1/3] Codegen: Add a pass that converts big integers into floats. This allows us to use Taylor expansions for the Laplace FMM with order >= 9. Without this, loopy will complain that it can't infer the type of large ints. --- sumpy/codegen.py | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/sumpy/codegen.py b/sumpy/codegen.py index 9c8a4958..0cb14028 100644 --- a/sumpy/codegen.py +++ b/sumpy/codegen.py @@ -430,6 +430,41 @@ class FractionKiller(CSECachingMapperMixin, IdentityMapper): # }}} +# {{{ convert big integers into floats + +from loopy.tools import is_integer + + +class BigIntegerKiller(CSECachingMapperMixin, IdentityMapper): + + def __init__(self, warn_on_first=True, int_type=np.int64, float_type=np.float64): + IdentityMapper.__init__(self) + self.bits = 64 + self.warn = warn_on_first + self.float_type = float_type + self.iinfo = np.iinfo(int_type) + + def map_constant(self, expr): + if not is_integer(expr): + return IdentityMapper.map_constant(self, expr) + + if self.iinfo.min <= expr <= self.iinfo.max: + return expr + + if self.warn: + from warnings import warn + warn("Converting '%d' to float: this may result in " + "loss of digits." % expr) + # Suppress further warnings. + self.warn = False + + return self.float_type(expr) + + map_common_subexpression_uncached = IdentityMapper.map_common_subexpression + +# }}} + + # {{{ vector component rewriter INDEXED_VAR_RE = re.compile("^([a-zA-Z_]+)([0-9]+)$") @@ -534,6 +569,7 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[], pwr = PowerRewriter() ssg = SumSignGrouper() fck = FractionKiller() + bik = BigIntegerKiller() def convert_expr(name, expr): logger.debug("generate expression for: %s" % name) @@ -543,6 +579,7 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[], expr = pwr(expr) expr = fck(expr) expr = ssg(expr) + expr = bik(expr) #expr = cse_tag(expr) for m in pymbolic_expr_maps: expr = m(expr) -- GitLab From 85653a69b1ccac1c266bc6a2a0e09dc3e6163349 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 26 Dec 2016 19:41:17 -0600 Subject: [PATCH 2/3] Specify types of nsources in expansions. --- sumpy/p2p.py | 4 ++-- sumpy/qbx.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/sumpy/p2p.py b/sumpy/p2p.py index 7234180f..0345f904 100644 --- a/sumpy/p2p.py +++ b/sumpy/p2p.py @@ -154,8 +154,8 @@ class P2P(P2PBase): shape=(self.dim, "nsources")), lp.GlobalArg("targets", None, shape=(self.dim, "ntargets")), - lp.ValueArg("nsources", None), - lp.ValueArg("ntargets", None), + lp.ValueArg("nsources", np.int32), + lp.ValueArg("ntargets", np.int32), lp.GlobalArg("strength", None, shape="nstrengths,nsources"), lp.GlobalArg("result", None, shape="nresults,ntargets", dim_tags="sep,C") diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 0b3dc035..629f9092 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -111,8 +111,8 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): shape=(self.dim, "ntargets"), order="C"), lp.GlobalArg("center", None, shape=(self.dim, "ntargets"), order="C"), - lp.ValueArg("nsources", None), - lp.ValueArg("ntargets", None), + lp.ValueArg("nsources", np.int32), + lp.ValueArg("ntargets", np.int32), ] @memoize_method -- GitLab From e250e5bba96e2ba251449980ee4e59f678dd1bf6 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 26 Dec 2016 19:41:32 -0600 Subject: [PATCH 3/3] Fix stresslet kernel to be FMMable. --- sumpy/kernel.py | 41 +++++++++++++++++++++++++---------------- 1 file changed, 25 insertions(+), 16 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 39f490a3..4a3987b5 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -506,26 +506,22 @@ class StressletKernel(ExpressionKernel): if dim == 2: d = make_sym_vector("d", dim) - n = make_sym_vector(stresslet_vector_name, dim) r = pymbolic_real_norm_2(d) expr = ( - sum(n[axis]*d[axis] for axis in range(dim)) - * - d[icomp]*d[jcomp]/r**4 + -var("log")(r)*(1 if icomp == jcomp else 0) + + + d[icomp]*d[jcomp]/r**2 ) - scaling = 1/(var("pi")) - + scaling = 1/(4*var("pi")) elif dim == 3: d = make_sym_vector("d", dim) - n = make_sym_vector(stresslet_vector_name, dim) r = pymbolic_real_norm_2(d) expr = ( - sum(n[axis]*d[axis] for axis in range(dim)) - * - d[icomp]*d[jcomp]/r**5 + (1/r)*(1 if icomp == jcomp else 0) + + + d[icomp]*d[jcomp]/r**3 ) - scaling = -3/(4*var("pi")) - + scaling = -1/(8*var("pi")) elif dim is None: expr = None scaling = None @@ -562,15 +558,28 @@ class StressletKernel(ExpressionKernel): return [ KernelArgument( loopy_arg=lp.ValueArg(self.viscosity_mu_name, np.float64), - ), + )] + + def get_source_args(self): + return [ KernelArgument( loopy_arg=lp.GlobalArg(self.stresslet_vector_name, None, shape=(self.dim, "nsources"), - dim_tags="sep,C")) - ] + dim_tags="sep,C"))] + + def postprocess_at_source(self, expr, avec): + dimensions = len(avec) + assert dimensions == self.dim - def get_code_tranformer(self): + from sumpy.symbolic import make_sympy_vector + n = make_sympy_vector(self.stresslet_vector_name, dimensions) + + # avec = center-src -> minus sign from chain rule + return sum(-n[axis]*expr.diff(avec[axis]) + for axis in range(dimensions)) + + def get_code_transformer(self): from sumpy.codegen import VectorComponentRewriter vcr = VectorComponentRewriter([self.stresslet_vector_name]) from pymbolic.primitives import Variable -- GitLab