From 6682405d0adb9fa6ccbe3e391a61e0317f042532 Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Tue, 27 Mar 2012 17:34:43 -0400
Subject: [PATCH] A zoo of little changes towards enabling the QUACK scheme
code.
---
setup.py | 2 +-
sumpy/symbolic/__init__.py | 55 +++++++++++++++++++++++++++++++++-----
sumpy/symbolic/codegen.py | 4 +++
sumpy/tools.py | 6 +++++
test/test_kernels.py | 4 +++
5 files changed, 63 insertions(+), 8 deletions(-)
diff --git a/setup.py b/setup.py
index c793e73c..58cc30f0 100644
--- a/setup.py
+++ b/setup.py
@@ -37,7 +37,7 @@ setup(name="sumpy",
author="Rio Yokota, Andreas Kloeckner, Matthew Knepley",
author_email="yokota@bu.edu",
license = "MIT",
- packages=["sumpy"],
+ packages=["sumpy", "sumpy.symbolic"],
# 2to3 invocation
cmdclass={'build_py': build_py})
diff --git a/sumpy/symbolic/__init__.py b/sumpy/symbolic/__init__.py
index 9c184b17..58d4f8d5 100644
--- a/sumpy/symbolic/__init__.py
+++ b/sumpy/symbolic/__init__.py
@@ -61,7 +61,7 @@ class IdentityMapper(SympyMapper):
map_Symbol = map_Integer
map_Real = map_Integer
-
+ map_ImaginaryUnit = map_Integer
@@ -129,6 +129,9 @@ class CSEMapper(IdentityMapper):
else:
return IdentityMapper.map_Rational(self, expr)
+ def map_Subs(self, expr):
+ return expr
+
@@ -154,7 +157,7 @@ def eliminate_common_subexpressions(exprs, sym_gen):
else:
iterable = node
- for arg in node.args:
+ for arg in iterable:
gather_use_counts(arg)
for expr in exprs:
@@ -204,20 +207,58 @@ def diff_multi_index(expr, multi_index, var_name):
-def make_coulomb_kernel_in(var_name, dimensions):
- from sumpy.symbolic import make_sym_vector
- dist = make_sym_vector(var_name, dimensions)
+def make_coulomb_kernel(dist_vec):
+ dimensions = len(dist_vec)
+ r = sp.sqrt((dist_vec.T*dist_vec)[0,0])
if dimensions == 2:
- return sp.log(sp.sqrt((dist.T*dist)[0,0]))
+ return sp.log(r)
elif dimensions == 3:
- return 1/sp.sqrt((dist.T*dist)[0,0])
+ return 1/r
else:
raise RuntimeError("unsupported dimensionality")
+
+def make_helmholtz_kernel(dist_vec):
+ dimensions = len(dist_vec)
+ r = sp.sqrt((dist_vec.T*dist_vec)[0,0])
+
+ i = sp.sqrt(-1)
+ k = sp.Symbol("k")
+
+ if dimensions == 2:
+ return i/4 * sp.Function("H1_0")(k*r)
+ elif dimensions == 3:
+ return sp.exp(i*k*r)/r
+ else:
+ raise RuntimeError("unsupported dimensionality")
+
+
+
+
+def make_coulomb_kernel_in(var_name, dimensions):
+ from warnings import warn
+ warn("make_coulomb_kernel_in is deprecated", DeprecationWarning, stacklevel=2)
+
+ from sumpy.symbolic import make_sym_vector
+ return make_coulomb_kernel(make_sym_vector(var_name, dimensions))
+
+
+
+
+def make_helmholtz_kernel_in(var_name, dimensions):
+ from warnings import warn
+ warn("make_helmholtz_kernel_in is deprecated", DeprecationWarning, stacklevel=2)
+
+ from sumpy.symbolic import make_sym_vector
+ return make_helmholtz_kernel(make_sym_vector(var_name, dimensions))
+
+
+
+
def make_coulomb_kernel_ts(dimensions):
old = make_sym_vector("d", dimensions)
new = (make_sym_vector("t", dimensions)
diff --git a/sumpy/symbolic/codegen.py b/sumpy/symbolic/codegen.py
index a30523ca..8b1523f3 100644
--- a/sumpy/symbolic/codegen.py
+++ b/sumpy/symbolic/codegen.py
@@ -88,6 +88,10 @@ class PowRewriter(IdentityMapper):
return IdentityMapper.map_Pow(self, expr)
+ def map_Subs(self, expr):
+ return expr
+
+
diff --git a/sumpy/tools.py b/sumpy/tools.py
index b9d6c8ba..f9f102c2 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -21,6 +21,12 @@ def mi_power(vector, mi):
result *= vec_i**mi_i
return result
+def mi_derivative(expr, vector, mi):
+ for mi_i, vec_i in zip(mi, vector):
+ expr = expr.diff(vec_i, mi_i)
+ return expr
+
+
# }}}
diff --git a/test/test_kernels.py b/test/test_kernels.py
index eb466aea..cb5d585e 100644
--- a/test/test_kernels.py
+++ b/test/test_kernels.py
@@ -156,6 +156,8 @@ def test_p2m2p(ctx_getter):
assert la.norm((potential_dev-potential_dev_direct).get())/res**2 < 1e-3
+ # }}}
+
@@ -263,6 +265,8 @@ def test_p2m2m2p(ctx_getter):
assert la.norm((potential_dev-potential_dev_direct).get())/res**2 < 1e-3
+ # }}}
+
--
GitLab