diff options
| -rw-r--r-- | symbolic/D2Q9.py | 2 | ||||
| -rw-r--r-- | symbolic/D3Q27.py | 2 | ||||
| -rw-r--r-- | symbolic/characteristics.py | 6 | ||||
| -rw-r--r-- | 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 |