From 8b376d86eb917f85bf51fda6fd750543aff62d0f Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Fri, 15 Apr 2016 15:30:35 -0500 Subject: [PATCH 1/2] Add 2D Stresslet kernel --- sumpy/kernel.py | 107 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 106 insertions(+), 1 deletion(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 3f745b70..a06e9702 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -222,7 +222,7 @@ class ExpressionKernel(Kernel): raise ValueError("dist_vec length does not match expected dimension") expr = expr.subs([ - ("d_%d" % i, dist_vec_i) + ("d%d" % i, dist_vec_i) for i, dist_vec_i in enumerate(dist_vec) ]) @@ -427,6 +427,84 @@ class StokesletKernel(ExpressionKernel): mapper_method = "map_stokeslet_kernel" + +class StressletKernel(ExpressionKernel): + init_arg_names = ("dim", "icomp", "jcomp", "viscosity_mu_name", + "stresslet_vector_name") + + def __init__(self, dim=None, icomp=None, jcomp=None, viscosity_mu_name="mu", + stresslet_vector_name="stresslet_vec"): + """ + :arg viscosity_mu_name: The argument name to use for + dynamic viscosity :math:`\mu` the then generating functions to + evaluate this kernel. + """ + mu = var(viscosity_mu_name) + #Remove mu? (Not needed inside kernel?) + + 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 + ) + scaling = 1/(var("pi")) + + elif dim is None: + expr = None + scaling = None + else: + raise RuntimeError("unsupported dimensionality") + + self.viscosity_mu_name = viscosity_mu_name + self.stresslet_vector_name = stresslet_vector_name + self.icomp = icomp + self.jcomp = jcomp + + ExpressionKernel.__init__( + self, + dim, + expression=expr, + scaling=scaling, + is_complex_valued=False) + + def __getinitargs__(self): + return (self._dim, self.icomp, self.jcomp, self.viscosity_mu_name, + self.stresslet_vector_name) + + def update_persistent_hash(self, key_hash, key_builder): + key_hash.update(type(self).__name__.encode()) + key_builder.rec(key_hash, (self.dim, self.icomp, self.jcomp, self.viscosity_mu_name, self.stresslet_vector_name)) + + def __repr__(self): + return "StressletKnl%dD_%d%d[%s]" % (self.dim, self.icomp, self.jcomp, + self.stresslet_vector_name) + + def get_args(self): + return [ + KernelArgument( + loopy_arg=lp.ValueArg(self.viscosity_mu_name, np.float64), + ), + KernelArgument( + lp.GlobalArg(self.stresslet_vector_name, + None, + shape=(self.dim, "nsources"), + dim_tags="sep,C"), + ) + ] + + def transform_to_code(self, expr): + from sumpy.codegen import VectorComponentRewriter + vcr = VectorComponentRewriter([self.stresslet_vector_name]) + from pymbolic.primitives import Variable + return _VectorIndexAdder(self.stresslet_vector_name, (Variable("isrc"),))( + vcr(expr)) + + mapper_method = "map_stresslet_kernel" + # }}} @@ -638,6 +716,8 @@ class KernelIdentityMapper(KernelMapper): map_laplace_kernel = map_expression_kernel map_helmholtz_kernel = map_expression_kernel + map_stokeslet_kernel = map_expression_kernel + map_stresslet_kernel = map_expression_kernel def map_axis_target_derivative(self, kernel): return AxisTargetDerivative(kernel.axis, self.rec(kernel.inner_kernel)) @@ -669,6 +749,8 @@ class DerivativeCounter(KernelCombineMapper): map_laplace_kernel = map_expression_kernel map_helmholtz_kernel = map_expression_kernel + map_stokeslet_kernel = map_expression_kernel + map_stresslet_kernel = map_expression_kernel def map_axis_target_derivative(self, kernel): return 1 + self.rec(kernel.inner_kernel) @@ -707,6 +789,29 @@ class KernelDimensionSetter(KernelIdentityMapper): helmholtz_k_name=kernel.helmholtz_k_name, allow_evanescent=kernel.allow_evanescent) + def map_stokeslet_kernel(self, kernel): + if kernel._dim is not None and self.dim != kernel.dim: + raise RuntimeError("cannot set kernel dimension to new value (%d)" + "different from existing one (%d)" + % (self.dim, kernel.dim)) + + return StressletKernel(self.dim, + kernel.icomp, + kernel.jcomp, + viscosity_mu_name=kernel.viscosity_mu_name) + + def map_stresslet_kernel(self, kernel): + if kernel._dim is not None and self.dim != kernel.dim: + raise RuntimeError("cannot set kernel dimension to new value (%d)" + "different from existing one (%d)" + % (self.dim, kernel.dim)) + + return StressletKernel(self.dim, + kernel.icomp, + kernel.jcomp, + viscosity_mu_name=kernel.viscosity_mu_name, + stresslet_vector_name=kernel.stresslet_vector_name) + # }}} -- GitLab From e4b8af802079f318ddf22f4e83d710300946fa29 Mon Sep 17 00:00:00 2001 From: Natalie Beams Date: Fri, 15 Apr 2016 15:33:10 -0500 Subject: [PATCH 2/2] Fix vector subscript mapping --- sumpy/symbolic.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py index b2ab466d..338ced14 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -31,6 +31,7 @@ import sympy as sp import numpy as np from pymbolic.mapper import IdentityMapper as IdentityMapperBase from pymbolic.sympy_interface import PymbolicToSympyMapper +import pymbolic.primitives as prim import logging logger = logging.getLogger(__name__) @@ -210,4 +211,11 @@ class PymbolicToSympyMapperWithSymbols(PymbolicToSympyMapper): else: return PymbolicToSympyMapper.map_variable(self, expr) + def map_subscript(self, expr): + if isinstance(expr.aggregate, prim.Variable) and isinstance(expr.index, int): + return sp.Symbol("%s%d" % (expr.aggregate.name, expr.index)) + else: + raise RuntimeError("do not know how to translate '%s' to sympy" + % expr) + # vim: fdm=marker -- GitLab