aboutsummaryrefslogtreecommitdiff
path: root/symbolic
diff options
context:
space:
mode:
authorAdrian Kummerlaender2019-07-02 14:58:36 +0200
committerAdrian Kummerlaender2019-07-02 14:59:46 +0200
commit7ee1ca337f25ab6ead86ead28c10df520cff1b63 (patch)
treec73d4f9ae857b39c45a99149a3f08fabba82e905 /symbolic
parentf62fcf3e47919b1c893ac2335b9fdfb842f9595f (diff)
downloadsymlbm_playground-7ee1ca337f25ab6ead86ead28c10df520cff1b63.tar
symlbm_playground-7ee1ca337f25ab6ead86ead28c10df520cff1b63.tar.gz
symlbm_playground-7ee1ca337f25ab6ead86ead28c10df520cff1b63.tar.bz2
symlbm_playground-7ee1ca337f25ab6ead86ead28c10df520cff1b63.tar.lz
symlbm_playground-7ee1ca337f25ab6ead86ead28c10df520cff1b63.tar.xz
symlbm_playground-7ee1ca337f25ab6ead86ead28c10df520cff1b63.tar.zst
symlbm_playground-7ee1ca337f25ab6ead86ead28c10df520cff1b63.zip
Determine weights using Gauss-Hermite quadrature
Diffstat (limited to 'symbolic')
-rw-r--r--symbolic/D2Q9.py1
-rw-r--r--symbolic/D3Q27.py6
-rw-r--r--symbolic/characteristics.py21
-rw-r--r--symbolic/generator.py14
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