1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
|
from sympy import *
from sympy.codegen.ast import Assignment
import symbolic.optimizations as optimizations
from symbolic.characteristics import weights, c_s
def assign(names, definitions):
return list(map(lambda x: Assignment(*x), zip(names, definitions)))
class LBM:
def __init__(self, descriptor):
self.descriptor = descriptor
self.f_next = symarray('f_next', descriptor.q)
self.f_curr = symarray('f_curr', descriptor.q)
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')
u = Matrix(symarray('u', self.descriptor.d))
exprs = [ Assignment(rho, sum(self.f_curr)) ]
for i, u_i in enumerate(u):
exprs.append(
Assignment(u_i, sum([ (c_j*self.f_curr[j])[i] for j, c_j in enumerate(self.descriptor.c) ]) / sum(self.f_curr)))
if optimize:
return cse(exprs, optimizations=optimizations.custom, symbols=numbered_symbols(prefix='m'))
else:
return ([], exprs)
def equilibrium(self):
rho = symbols('rho')
u = Matrix(symarray('u', self.descriptor.d))
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.append(f_eq_i)
return f_eq
def bgk(self, tau, f_eq, optimize = True):
exprs = [ self.f_curr[i] + 1/tau * (f_eq_i - self.f_curr[i]) for i, f_eq_i in enumerate(f_eq) ]
if optimize:
subexprs, f = cse(exprs, optimizations=optimizations.custom)
return (subexprs, assign(self.f_next, f))
else:
return ([], assign(self.f_next, exprs))
|
