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
|
from sympy import *
from sympy.codegen.ast import Assignment
q = 19
d = 3
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))
f_next = symarray('f_next', q)
f_curr = symarray('f_curr', q)
def moments(f = f_curr, optimize = True):
rho = symbols('rho')
u = Matrix(symarray('u', d))
exprs = [ Assignment(rho, sum(f)) ]
for i, u_i in enumerate(u):
exprs.append(Assignment(u_i, sum([ (c_j*f[j])[i] for j, c_j in enumerate(c) ]) / sum(f)))
if optimize:
return cse(exprs, optimizations='basic', symbols=numbered_symbols(prefix='m'))
else:
return ([], exprs)
def equilibrium():
rho = symbols('rho')
u = Matrix(symarray('u', d))
f_eq = []
for i, c_i in enumerate(c):
f_eq_i = w[i] * rho * ( 1
+ c_i.dot(u) / c_s**2
+ c_i.dot(u)**2 / (2*c_s**4)
- u.dot(u) / (2*c_s**2) )
f_eq.append(f_eq_i)
return f_eq
def bgk(tau, f_eq = equilibrium(), optimize = True):
exprs = [ Assignment(f_next[i], f_curr[i] + 1/tau * ( f_eq_i - f_curr[i] )) for i, f_eq_i in enumerate(f_eq) ]
if optimize:
return cse(exprs, optimizations='basic')
else:
return ([], exprs)
|
