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/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 +++++++++++++++++++++++++++
 6 files changed, 83 insertions(+)
 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

(limited to 'boltzgen/lbm/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()
-- 
cgit v1.2.3