From 24847cbb2567f508a7c30b39c6fb7ba6379d1adc Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Sat, 2 Nov 2019 17:18:32 +0100
Subject: Restructure LBM model / lattice distinction

---
 boltzgen.py                                        | 11 ++-
 boltzgen/__init__.py                               |  2 +-
 .../kernel/template/bounce_back_boundary.cl.mako   |  2 +-
 .../kernel/template/bounce_back_boundary.cpp.mako  |  2 +-
 .../kernel/template/collide_and_stream.cl.mako     |  2 +-
 .../kernel/template/collide_and_stream.cpp.mako    |  2 +-
 boltzgen/kernel/template/momenta_boundary.cl.mako  |  4 +-
 boltzgen/kernel/template/momenta_boundary.cpp.mako |  2 +-
 boltzgen/lbm/__init__.py                           | 80 +---------------------
 boltzgen/lbm/lattice/D2Q9.py                       |  7 ++
 boltzgen/lbm/lattice/D3Q19.py                      | 18 +++++
 boltzgen/lbm/lattice/D3Q27.py                      |  7 ++
 boltzgen/lbm/lattice/D3Q7.py                       | 18 +++++
 boltzgen/lbm/lattice/__init__.py                   |  6 ++
 boltzgen/lbm/lattice/characteristics.py            | 27 ++++++++
 boltzgen/lbm/model/BGK.py                          | 76 ++++++++++++++++++++
 boltzgen/lbm/model/D2Q9.py                         |  7 --
 boltzgen/lbm/model/D3Q19.py                        | 18 -----
 boltzgen/lbm/model/D3Q27.py                        |  7 --
 boltzgen/lbm/model/D3Q7.py                         | 18 -----
 boltzgen/lbm/model/__init__.py                     |  7 +-
 boltzgen/lbm/model/characteristics.py              | 27 --------
 boltzgen/utility/__init__.py                       |  5 ++
 23 files changed, 183 insertions(+), 172 deletions(-)
 create mode 100644 boltzgen/lbm/lattice/D2Q9.py
 create mode 100644 boltzgen/lbm/lattice/D3Q19.py
 create mode 100644 boltzgen/lbm/lattice/D3Q27.py
 create mode 100644 boltzgen/lbm/lattice/D3Q7.py
 create mode 100644 boltzgen/lbm/lattice/__init__.py
 create mode 100644 boltzgen/lbm/lattice/characteristics.py
 create mode 100644 boltzgen/lbm/model/BGK.py
 delete mode 100644 boltzgen/lbm/model/D2Q9.py
 delete mode 100644 boltzgen/lbm/model/D3Q19.py
 delete mode 100644 boltzgen/lbm/model/D3Q27.py
 delete mode 100644 boltzgen/lbm/model/D3Q7.py
 delete mode 100644 boltzgen/lbm/model/characteristics.py

diff --git a/boltzgen.py b/boltzgen.py
index 82adafe..4317c68 100755
--- a/boltzgen.py
+++ b/boltzgen.py
@@ -8,7 +8,8 @@ argparser = argparse.ArgumentParser(
 
 argparser.add_argument('target', help = 'Target language (currently either "cl" or "cpp")')
 
-argparser.add_argument('--lattice',   required = True,  help = 'Lattice type (D2Q9, D3Q7, D3Q19, D3Q27)')
+argparser.add_argument('--lattice',   required = True,  help = 'Lattice type ("D2Q9", "D3Q7", "D3Q19", "D3Q27")')
+argparser.add_argument('--model',     required = False, help = 'LBM model (currently only "BGK")')
 argparser.add_argument('--layout',    required = True,  help = 'Memory layout ("AOS" or "SOA")')
 argparser.add_argument('--index',     required = False, help = 'Cell indexing ("XYZ" or "ZYX")')
 argparser.add_argument('--precision', required = True,  help = 'Floating precision ("single" or "double")')
@@ -21,13 +22,17 @@ argparser.add_argument('--extras',      action = 'append', nargs = '+', default
 
 args = argparser.parse_args()
 
-lattice = eval("lbm.model.%s" % args.lattice)
+if args.model is None:
+    args.model = "BGK"
+
+lattice = eval("lbm.lattice.%s" % args.lattice)
+model   = eval("lbm.model.%s"   % args.model)
 
 if args.index is None:
     args.index = 'XYZ'
 
 generator = Generator(
-    model      = LBM(lattice, tau = float(args.tau), optimize = not args.disable_cse),
+    model      = model(lattice, tau = float(args.tau), optimize = not args.disable_cse),
     target     = args.target,
     precision  = args.precision,
     index      = args.index,
diff --git a/boltzgen/__init__.py b/boltzgen/__init__.py
index 01fc32c..d319fdc 100644
--- a/boltzgen/__init__.py
+++ b/boltzgen/__init__.py
@@ -1,5 +1,5 @@
+import boltzgen.lbm.lattice
 import boltzgen.lbm.model
 
-from boltzgen.lbm              import LBM
 from boltzgen.geometry         import Geometry
 from boltzgen.kernel.generator import Generator
diff --git a/boltzgen/kernel/template/bounce_back_boundary.cl.mako b/boltzgen/kernel/template/bounce_back_boundary.cl.mako
index 8bb4f6c..4766bba 100644
--- a/boltzgen/kernel/template/bounce_back_boundary.cl.mako
+++ b/boltzgen/kernel/template/bounce_back_boundary.cl.mako
@@ -10,7 +10,7 @@ __kernel void bounce_back_boundary_gid(__global ${float_type}* f_next,
 % endfor
 
 <%
-    subexpr, assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = True))
+    subexpr, assignment = model.collision(f_eq = model.equilibrium(resolve_moments = True))
 %>
 
 % for i, expr in enumerate(subexpr):
diff --git a/boltzgen/kernel/template/bounce_back_boundary.cpp.mako b/boltzgen/kernel/template/bounce_back_boundary.cpp.mako
index 91dcfa0..c37233e 100644
--- a/boltzgen/kernel/template/bounce_back_boundary.cpp.mako
+++ b/boltzgen/kernel/template/bounce_back_boundary.cpp.mako
@@ -10,7 +10,7 @@ void bounce_back_boundary(      ${float_type}* f_next,
 % endfor
 
 <%
-    subexpr, assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = True))
+    subexpr, assignment = model.collision(f_eq = model.equilibrium(resolve_moments = True))
 %>
 
 % for i, expr in enumerate(subexpr):
diff --git a/boltzgen/kernel/template/collide_and_stream.cl.mako b/boltzgen/kernel/template/collide_and_stream.cl.mako
index 303b21b..6b55cd2 100644
--- a/boltzgen/kernel/template/collide_and_stream.cl.mako
+++ b/boltzgen/kernel/template/collide_and_stream.cl.mako
@@ -10,7 +10,7 @@ __kernel void collide_and_stream_gid(__global ${float_type}* f_next,
 % endfor
 
 <%
-    subexpr, assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = True))
+    subexpr, assignment = model.collision(f_eq = model.equilibrium(resolve_moments = True))
 %>
 
 % for i, expr in enumerate(subexpr):
diff --git a/boltzgen/kernel/template/collide_and_stream.cpp.mako b/boltzgen/kernel/template/collide_and_stream.cpp.mako
index 47dd6dd..2c52b54 100644
--- a/boltzgen/kernel/template/collide_and_stream.cpp.mako
+++ b/boltzgen/kernel/template/collide_and_stream.cpp.mako
@@ -10,7 +10,7 @@ void collide_and_stream(      ${float_type}* f_next,
 % endfor
 
 <%
-    subexpr, assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = True))
+    subexpr, assignment = model.collision(f_eq = model.equilibrium(resolve_moments = True))
 %>
 
 % for i, expr in enumerate(subexpr):
diff --git a/boltzgen/kernel/template/momenta_boundary.cl.mako b/boltzgen/kernel/template/momenta_boundary.cl.mako
index 4e29601..7c1e3df 100644
--- a/boltzgen/kernel/template/momenta_boundary.cl.mako
+++ b/boltzgen/kernel/template/momenta_boundary.cl.mako
@@ -1,6 +1,6 @@
 <%
 moments_subexpr, moments_assignment = model.moments()
-collision_subexpr, collision_assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = False))
+collision_subexpr, collision_assignment = model.collision(f_eq = model.equilibrium(resolve_moments = False))
 %>
 
 <%def name="momenta_boundary(name, param)">
@@ -30,7 +30,7 @@ __kernel void ${name}_momenta_boundary_gid(
     const ${float_type} ${ccode(expr)}
 % endfor
 
-% for i in range(0,descriptor.q):
+% for i, expr in enumerate(collision_assignment):
     preshifted_f_next[${layout.pop_offset(i)}] = f_next_${i};
 % endfor
 }
diff --git a/boltzgen/kernel/template/momenta_boundary.cpp.mako b/boltzgen/kernel/template/momenta_boundary.cpp.mako
index 6d05288..c41d07e 100644
--- a/boltzgen/kernel/template/momenta_boundary.cpp.mako
+++ b/boltzgen/kernel/template/momenta_boundary.cpp.mako
@@ -1,6 +1,6 @@
 <%
 moments_subexpr, moments_assignment = model.moments()
-collision_subexpr, collision_assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = False))
+collision_subexpr, collision_assignment = model.collision(f_eq = model.equilibrium(resolve_moments = False))
 %>
 
 <%def name="momenta_boundary(name, param)">
diff --git a/boltzgen/lbm/__init__.py b/boltzgen/lbm/__init__.py
index 3f608cb..276c0bb 100644
--- a/boltzgen/lbm/__init__.py
+++ b/boltzgen/lbm/__init__.py
@@ -1,78 +1,2 @@
-from sympy import *
-from sympy.codegen.ast import Assignment
-
-import boltzgen.utility.optimizations as optimizations
-from boltzgen.lbm.model.characteristics import weights, c_s
-
-
-def assign(names, definitions):
-    return list(map(lambda x: Assignment(*x), zip(names, definitions)))
-
-class LBM:
-    def __init__(self, descriptor, tau, optimize = True):
-        self.descriptor = descriptor
-        self.tau        = tau
-        self.optimize   = optimize
-
-        if self.tau <= 0.5:
-            raise Exception('Relaxation time must be larger than 0.5')
-
-        self.f_next = symarray('f_next', descriptor.q)
-        self.f_curr = symarray('f_curr', descriptor.q)
-
-        if not hasattr(descriptor, 'w'):
-            self.descriptor.w = weights(descriptor.d, descriptor.c)
-
-        if not hasattr(descriptor, 'c_s'):
-            self.descriptor.c_s = c_s(descriptor.d, descriptor.c, self.descriptor.w)
-
-    def moments(self, optimize = None):
-        if optimize is None:
-            optimize = self.optimize
-
-        rho = symbols('rho')
-        u   = Matrix(symarray('u', self.descriptor.d))
-
-        exprs = [ Assignment(rho, sum(self.f_curr)) ]
-
-        for i, u_i in enumerate(u):
-            exprs.append(
-                Assignment(u_i, sum([ (c_j*self.f_curr[j])[i] for j, c_j in enumerate(self.descriptor.c) ]) / sum(self.f_curr)))
-
-        if optimize:
-            return cse(exprs, optimizations=optimizations.custom, symbols=numbered_symbols(prefix='m'))
-        else:
-            return ([], exprs)
-
-    def equilibrium(self, resolve_moments = False):
-        rho = symbols('rho')
-        u   = Matrix(symarray('u', self.descriptor.d))
-
-        if resolve_moments:
-            moments = self.moments(optimize = False)[1]
-            rho = moments[0].rhs
-            for i, m in enumerate(moments[1:]):
-                u[i] = m.rhs
-
-        f_eq = []
-
-        for i, c_i in enumerate(self.descriptor.c):
-            f_eq_i = self.descriptor.w[i] * rho * ( 1
-                                                  + c_i.dot(u)    /    self.descriptor.c_s**2
-                                                  + c_i.dot(u)**2 / (2*self.descriptor.c_s**4)
-                                                  - u.dot(u)      / (2*self.descriptor.c_s**2) )
-            f_eq.append(f_eq_i)
-
-        return f_eq
-
-    def bgk(self, f_eq, optimize = None):
-        if optimize is None:
-            optimize = self.optimize
-
-        exprs = [ self.f_curr[i] + 1/self.tau * (f_eq_i - self.f_curr[i]) for i, f_eq_i in enumerate(f_eq) ]
-
-        if optimize:
-            subexprs, f = cse(exprs, optimizations=optimizations.custom)
-            return (subexprs, assign(self.f_next, f))
-        else:
-            return ([], assign(self.f_next, exprs))
+from . import model
+from . import lattice
diff --git a/boltzgen/lbm/lattice/D2Q9.py b/boltzgen/lbm/lattice/D2Q9.py
new file mode 100644
index 0000000..e3ac9de
--- /dev/null
+++ b/boltzgen/lbm/lattice/D2Q9.py
@@ -0,0 +1,7 @@
+from sympy import Matrix
+from itertools import product
+
+d = 2
+q = 9
+
+c = [ Matrix(x) for x in product([-1,0,1], repeat=d) ]
diff --git a/boltzgen/lbm/lattice/D3Q19.py b/boltzgen/lbm/lattice/D3Q19.py
new file mode 100644
index 0000000..e9e6eec
--- /dev/null
+++ b/boltzgen/lbm/lattice/D3Q19.py
@@ -0,0 +1,18 @@
+from sympy import Matrix, Rational, sqrt
+
+d = 3
+q = 19
+
+c = [ Matrix(x) for x in [
+    ( 0, 1, 1), (-1, 0, 1), ( 0, 0, 1), ( 1, 0, 1), ( 0, -1, 1),
+    (-1, 1, 0), ( 0, 1, 0), ( 1, 1, 0), (-1, 0, 0), ( 0, 0, 0), ( 1, 0, 0), (-1,-1, 0), ( 0, -1, 0), ( 1, -1, 0),
+    ( 0, 1,-1), (-1, 0,-1), ( 0, 0,-1), ( 1, 0,-1), ( 0, -1,-1)
+]]
+
+w = [Rational(*x) for x in [
+    (1,36), (1,36), (1,18), (1,36), (1,36),
+    (1,36), (1,18), (1,36), (1,18), (1,3), (1,18), (1,36), (1,18), (1,36),
+    (1,36), (1,36), (1,18), (1,36), (1,36)
+]]
+
+c_s = sqrt(Rational(1,3))
diff --git a/boltzgen/lbm/lattice/D3Q27.py b/boltzgen/lbm/lattice/D3Q27.py
new file mode 100644
index 0000000..6fb1f80
--- /dev/null
+++ b/boltzgen/lbm/lattice/D3Q27.py
@@ -0,0 +1,7 @@
+from sympy import Matrix
+from itertools import product
+
+d = 3
+q = 27
+
+c = [ Matrix(x) for x in product([-1,0,1], repeat=d) ]
diff --git a/boltzgen/lbm/lattice/D3Q7.py b/boltzgen/lbm/lattice/D3Q7.py
new file mode 100644
index 0000000..04e16a3
--- /dev/null
+++ b/boltzgen/lbm/lattice/D3Q7.py
@@ -0,0 +1,18 @@
+from sympy import *
+
+q = 7
+d = 3
+
+c = [ Matrix(x) for x in [
+    ( 0, 0, 1),
+    ( 0, 1, 0), (-1, 0, 0), ( 0, 0, 0), ( 1, 0, 0), ( 0, -1, 0),
+    ( 0, 0,-1)
+]]
+
+w = [Rational(*x) for x in [
+    (1,8),
+    (1,8), (1,8), (1,4), (1,8), (1,8),
+    (1,8)
+]]
+
+c_s = sqrt(Rational(1,4))
diff --git a/boltzgen/lbm/lattice/__init__.py b/boltzgen/lbm/lattice/__init__.py
new file mode 100644
index 0000000..c956759
--- /dev/null
+++ b/boltzgen/lbm/lattice/__init__.py
@@ -0,0 +1,6 @@
+from . import D2Q9
+from . import D3Q7
+from . import D3Q19
+from . import D3Q27
+
+from . import characteristics
diff --git a/boltzgen/lbm/lattice/characteristics.py b/boltzgen/lbm/lattice/characteristics.py
new file mode 100644
index 0000000..b68afeb
--- /dev/null
+++ b/boltzgen/lbm/lattice/characteristics.py
@@ -0,0 +1,27 @@
+from sympy import *
+
+# copy of `sympy.integrals.quadrature.gauss_hermite` sans evaluation
+def gauss_hermite(n):
+    x = Dummy("x")
+    p  = hermite_poly(n, x, polys=True)
+    p1 = hermite_poly(n-1, x, polys=True)
+    xi = []
+    w  = []
+    for r in p.real_roots():
+        xi.append(r)
+        w.append(((2**(n-1) * factorial(n) * sqrt(pi))/(n**2 * p1.subs(x, r)**2)))
+    return xi, w
+
+# determine weights of a d-dimensional LBM model on velocity set c
+# (only works for velocity sets that result into NSE-recovering LB models when
+#  plugged into Gauss-Hermite quadrature without any additional arguments
+#  i.e. D2Q9 and D3Q27 but not D3Q19)
+def weights(d, c):
+    _, omegas = gauss_hermite(3)
+    return list(map(lambda c_i: Mul(*[ omegas[1+c_i[iDim]] for iDim in range(0,d) ]) / pi**(d/2), c))
+
+# determine lattice speed of sound using directions and their weights
+def c_s(d, c, w):
+    speeds = set([ sqrt(sum([ w[i] * c_i[j]**2 for i, c_i in enumerate(c) ])) for j in range(0,d) ])
+    assert len(speeds) == 1 # verify isotropy
+    return speeds.pop()
diff --git a/boltzgen/lbm/model/BGK.py b/boltzgen/lbm/model/BGK.py
new file mode 100644
index 0000000..f5a34e2
--- /dev/null
+++ b/boltzgen/lbm/model/BGK.py
@@ -0,0 +1,76 @@
+from sympy import *
+from sympy.codegen.ast import Assignment
+
+from boltzgen.utility import assign
+import boltzgen.utility.optimizations as optimizations
+
+from boltzgen.lbm.lattice.characteristics import weights, c_s
+
+class BGK:
+    def __init__(self, descriptor, tau, optimize = True):
+        self.descriptor = descriptor
+        self.tau        = tau
+        self.optimize   = optimize
+
+        if self.tau <= 0.5:
+            raise Exception('Relaxation time must be larger than 0.5')
+
+        self.f_next = symarray('f_next', descriptor.q)
+        self.f_curr = symarray('f_curr', descriptor.q)
+
+        if not hasattr(descriptor, 'w'):
+            self.descriptor.w = weights(descriptor.d, descriptor.c)
+
+        if not hasattr(descriptor, 'c_s'):
+            self.descriptor.c_s = c_s(descriptor.d, descriptor.c, self.descriptor.w)
+
+    def moments(self, optimize = None):
+        if optimize is None:
+            optimize = self.optimize
+
+        rho = symbols('rho')
+        u   = Matrix(symarray('u', self.descriptor.d))
+
+        exprs = [ Assignment(rho, sum(self.f_curr)) ]
+
+        for i, u_i in enumerate(u):
+            exprs.append(
+                Assignment(u_i, sum([ (c_j*self.f_curr[j])[i] for j, c_j in enumerate(self.descriptor.c) ]) / sum(self.f_curr)))
+
+        if optimize:
+            return cse(exprs, optimizations=optimizations.custom, symbols=numbered_symbols(prefix='m'))
+        else:
+            return ([], exprs)
+
+    def equilibrium(self, resolve_moments = False):
+        rho = symbols('rho')
+        u   = Matrix(symarray('u', self.descriptor.d))
+
+        if resolve_moments:
+            moments = self.moments(optimize = False)[1]
+            rho = moments[0].rhs
+            for i, m in enumerate(moments[1:]):
+                u[i] = m.rhs
+
+        f_eq = []
+
+        for i, c_i in enumerate(self.descriptor.c):
+            f_eq_i = self.descriptor.w[i] * rho * ( 1
+                                                  + c_i.dot(u)    /    self.descriptor.c_s**2
+                                                  + c_i.dot(u)**2 / (2*self.descriptor.c_s**4)
+                                                  - u.dot(u)      / (2*self.descriptor.c_s**2) )
+            f_eq.append(f_eq_i)
+
+        return f_eq
+
+    def collision(self, f_eq, optimize = None):
+        if optimize is None:
+            optimize = self.optimize
+
+        exprs = [ self.f_curr[i] + 1/self.tau * (f_eq_i - self.f_curr[i]) for i, f_eq_i in enumerate(f_eq) ]
+
+        if optimize:
+            subexprs, f = cse(exprs, optimizations=optimizations.custom)
+            return (subexprs, assign(self.f_next, f))
+        else:
+            return ([], assign(self.f_next, exprs))
diff --git a/boltzgen/lbm/model/D2Q9.py b/boltzgen/lbm/model/D2Q9.py
deleted file mode 100644
index e3ac9de..0000000
--- a/boltzgen/lbm/model/D2Q9.py
+++ /dev/null
@@ -1,7 +0,0 @@
-from sympy import Matrix
-from itertools import product
-
-d = 2
-q = 9
-
-c = [ Matrix(x) for x in product([-1,0,1], repeat=d) ]
diff --git a/boltzgen/lbm/model/D3Q19.py b/boltzgen/lbm/model/D3Q19.py
deleted file mode 100644
index e9e6eec..0000000
--- a/boltzgen/lbm/model/D3Q19.py
+++ /dev/null
@@ -1,18 +0,0 @@
-from sympy import Matrix, Rational, sqrt
-
-d = 3
-q = 19
-
-c = [ Matrix(x) for x in [
-    ( 0, 1, 1), (-1, 0, 1), ( 0, 0, 1), ( 1, 0, 1), ( 0, -1, 1),
-    (-1, 1, 0), ( 0, 1, 0), ( 1, 1, 0), (-1, 0, 0), ( 0, 0, 0), ( 1, 0, 0), (-1,-1, 0), ( 0, -1, 0), ( 1, -1, 0),
-    ( 0, 1,-1), (-1, 0,-1), ( 0, 0,-1), ( 1, 0,-1), ( 0, -1,-1)
-]]
-
-w = [Rational(*x) for x in [
-    (1,36), (1,36), (1,18), (1,36), (1,36),
-    (1,36), (1,18), (1,36), (1,18), (1,3), (1,18), (1,36), (1,18), (1,36),
-    (1,36), (1,36), (1,18), (1,36), (1,36)
-]]
-
-c_s = sqrt(Rational(1,3))
diff --git a/boltzgen/lbm/model/D3Q27.py b/boltzgen/lbm/model/D3Q27.py
deleted file mode 100644
index 6fb1f80..0000000
--- a/boltzgen/lbm/model/D3Q27.py
+++ /dev/null
@@ -1,7 +0,0 @@
-from sympy import Matrix
-from itertools import product
-
-d = 3
-q = 27
-
-c = [ Matrix(x) for x in product([-1,0,1], repeat=d) ]
diff --git a/boltzgen/lbm/model/D3Q7.py b/boltzgen/lbm/model/D3Q7.py
deleted file mode 100644
index 04e16a3..0000000
--- a/boltzgen/lbm/model/D3Q7.py
+++ /dev/null
@@ -1,18 +0,0 @@
-from sympy import *
-
-q = 7
-d = 3
-
-c = [ Matrix(x) for x in [
-    ( 0, 0, 1),
-    ( 0, 1, 0), (-1, 0, 0), ( 0, 0, 0), ( 1, 0, 0), ( 0, -1, 0),
-    ( 0, 0,-1)
-]]
-
-w = [Rational(*x) for x in [
-    (1,8),
-    (1,8), (1,8), (1,4), (1,8), (1,8),
-    (1,8)
-]]
-
-c_s = sqrt(Rational(1,4))
diff --git a/boltzgen/lbm/model/__init__.py b/boltzgen/lbm/model/__init__.py
index c956759..fb8f12b 100644
--- a/boltzgen/lbm/model/__init__.py
+++ b/boltzgen/lbm/model/__init__.py
@@ -1,6 +1 @@
-from . import D2Q9
-from . import D3Q7
-from . import D3Q19
-from . import D3Q27
-
-from . import characteristics
+from .BGK import BGK
diff --git a/boltzgen/lbm/model/characteristics.py b/boltzgen/lbm/model/characteristics.py
deleted file mode 100644
index b68afeb..0000000
--- a/boltzgen/lbm/model/characteristics.py
+++ /dev/null
@@ -1,27 +0,0 @@
-from sympy import *
-
-# copy of `sympy.integrals.quadrature.gauss_hermite` sans evaluation
-def gauss_hermite(n):
-    x = Dummy("x")
-    p  = hermite_poly(n, x, polys=True)
-    p1 = hermite_poly(n-1, x, polys=True)
-    xi = []
-    w  = []
-    for r in p.real_roots():
-        xi.append(r)
-        w.append(((2**(n-1) * factorial(n) * sqrt(pi))/(n**2 * p1.subs(x, r)**2)))
-    return xi, w
-
-# determine weights of a d-dimensional LBM model on velocity set c
-# (only works for velocity sets that result into NSE-recovering LB models when
-#  plugged into Gauss-Hermite quadrature without any additional arguments
-#  i.e. D2Q9 and D3Q27 but not D3Q19)
-def weights(d, c):
-    _, omegas = gauss_hermite(3)
-    return list(map(lambda c_i: Mul(*[ omegas[1+c_i[iDim]] for iDim in range(0,d) ]) / pi**(d/2), c))
-
-# determine lattice speed of sound using directions and their weights
-def c_s(d, c, w):
-    speeds = set([ sqrt(sum([ w[i] * c_i[j]**2 for i, c_i in enumerate(c) ])) for j in range(0,d) ])
-    assert len(speeds) == 1 # verify isotropy
-    return speeds.pop()
diff --git a/boltzgen/utility/__init__.py b/boltzgen/utility/__init__.py
index 8927868..fa9c760 100644
--- a/boltzgen/utility/__init__.py
+++ b/boltzgen/utility/__init__.py
@@ -1,2 +1,7 @@
 from . import optimizations
 from . import ndindex
+
+from sympy.codegen.ast import Assignment
+
+def assign(names, definitions):
+    return list(map(lambda x: Assignment(*x), zip(names, definitions)))
-- 
cgit v1.2.3