Determine weights using Gauss-Hermite quadrature

4 files changed, 31 insertions, 11 deletions

diff --git a/symbolic/D2Q9.py b/symbolic/D2Q9.py
index 22f7ed5..24fb223 100644
--- a/symbolic/D2Q9.py
+++ b/symbolic/D2Q9.py
@@ -4,6 +4,5 @@ q = 9
 d = 2
 
 c = [ Matrix(x) for x in [(-1, 1), ( 0, 1), ( 1, 1), (-1, 0), ( 0, 0), ( 1, 0), (-1,-1), ( 0, -1), ( 1, -1)] ]
-w = [ Rational(*x) for x in [(1,36), (1,9), (1,36), (1,9), (4,9), (1,9), (1,36), (1,9), (1,36)] ]
 
 c_s = sqrt(Rational(1,3))
diff --git a/symbolic/D3Q27.py b/symbolic/D3Q27.py
index 6a81de5..3db45d3 100644
--- a/symbolic/D3Q27.py
+++ b/symbolic/D3Q27.py
@@ -9,10 +9,4 @@ c = [ Matrix(x) for x in [
     (-1, 1,-1), ( 0, 1,-1), ( 1, 1,-1), (-1, 0,-1), ( 0, 0,-1), ( 1, 0,-1), (-1,-1,-1), ( 0, -1,-1), ( 1, -1,-1)
 ]]
 
-w = [Rational(*x) for x in [
-    (1, 216), (1,54), (1,216), (1,54), (2,27), (1,54), (1,216), (1,54), (1,216),
-    (1,  54), (2,27), (1, 54), (2,27), (8,27), (2,27), (1, 54), (2,27), (1, 54),
-    (1, 216), (1,54), (1,216), (1,54), (2,27), (1,54), (1,216), (1,54), (1,216)
-]]
-
 c_s = sqrt(Rational(1,3))
diff --git a/symbolic/characteristics.py b/symbolic/characteristics.py
new file mode 100644
index 0000000..ca3904e
--- /dev/null
+++ b/symbolic/characteristics.py
@@ -0,0 +1,21 @@
+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))
diff --git a/symbolic/generator.py b/symbolic/generator.py
index f94031f..73f3940 100644
--- a/symbolic/generator.py
+++ b/symbolic/generator.py
@@ -2,6 +2,7 @@ from sympy import *
 from sympy.codegen.ast import Assignment
 
 import symbolic.optimizations as optimizations
+from symbolic.characteristics import weights
 
 
 def assign(names, definitions):
@@ -13,6 +14,11 @@ class LBM:
         self.f_next = symarray('f_next', descriptor.q)
         self.f_curr = symarray('f_curr', descriptor.q)
 
+        if hasattr(descriptor, 'w'):
+            self.w = descriptor.w
+        else:
+            self.w = weights(descriptor.d, descriptor.c)
+
     def moments(self, optimize = True):
         rho = symbols('rho')
         u   = Matrix(symarray('u', self.descriptor.d))
@@ -35,10 +41,10 @@ class LBM:
         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_i = self.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