Determine weights using Gauss-Hermite quadrature

1 files changed, 21 insertions, 0 deletions

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))