From 55fb5b62ab3381327d4b04cf1c5ab095874d29cc Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Tue, 2 Jul 2019 20:35:44 +0200
Subject: Determine lattice speed of sound
---
symbolic/D2Q9.py | 2 --
symbolic/D3Q27.py | 2 --
symbolic/characteristics.py | 6 ++++++
symbolic/generator.py | 19 ++++++++++---------
4 files changed, 16 insertions(+), 13 deletions(-)
diff --git a/symbolic/D2Q9.py b/symbolic/D2Q9.py
index 24fb223..9693477 100644
--- a/symbolic/D2Q9.py
+++ b/symbolic/D2Q9.py
@@ -4,5 +4,3 @@ 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)] ]
-
-c_s = sqrt(Rational(1,3))
diff --git a/symbolic/D3Q27.py b/symbolic/D3Q27.py
index 3db45d3..d63aad8 100644
--- a/symbolic/D3Q27.py
+++ b/symbolic/D3Q27.py
@@ -8,5 +8,3 @@ c = [ Matrix(x) for x in [
(-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),
(-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)
]]
-
-c_s = sqrt(Rational(1,3))
diff --git a/symbolic/characteristics.py b/symbolic/characteristics.py
index ca3904e..b68afeb 100644
--- a/symbolic/characteristics.py
+++ b/symbolic/characteristics.py
@@ -19,3 +19,9 @@ def gauss_hermite(n):
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/symbolic/generator.py b/symbolic/generator.py
index 73f3940..629c991 100644
--- a/symbolic/generator.py
+++ b/symbolic/generator.py
@@ -2,7 +2,7 @@ from sympy import *
from sympy.codegen.ast import Assignment
import symbolic.optimizations as optimizations
-from symbolic.characteristics import weights
+from symbolic.characteristics import weights, c_s
def assign(names, definitions):
@@ -14,10 +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)
+ 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 = True):
rho = symbols('rho')
@@ -41,10 +42,10 @@ class LBM:
f_eq = []
for i, c_i in enumerate(self.descriptor.c):
- 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_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
--
cgit v1.2.3