In [323]:
from itertools import combinations, combinations_with_replacement
In [326]:
import numpy as np
In [7]:
A_R__1 = AffineSpace(1,RR,names='t')
In [8]:
A_R__2.<x,y> = AffineSpace(2,RR)
In [9]:
A_R__3.<x,y,z> = AffineSpace(3,RR)
In [14]:
[_ for _ in A_R__3.points_of_bounded_height(3.0)]
In [16]:
A_R__1.<t> = AffineSpace(1,RealField())
In [17]:
A_R__2.<x,y> = AffineSpace(2,RealField())
In [18]:
A_R__3.<x,y,z> = AffineSpace(3,RealField())
In [21]:
A_R__1.coordinate_ring()
Out[21]:
In [22]:
A_R__1.coordinate_ring().gens()
Out[22]:
In [23]:
A_R__2.coordinate_ring()
Out[23]:
In [24]:
A_R__2.coordinate_ring().gens()
Out[24]:
In [25]:
A_R__3.coordinate_ring()
Out[25]:
In [26]:
A_R__3.coordinate_ring().gens()
Out[26]:
In [29]:
A_R__1.coordinate_ring().monomial()
Out[29]:
In [30]:
A_R__2.coordinate_ring().monomial()
Out[30]:
In [31]:
A_R__3.coordinate_ring().monomial()
Out[31]:
In [35]:
cf. http://staff.polito.it/pietro.asinari/publications/preprint_Asinari_PA_2010a.pdf, I. Karlin and P. Asinari, Factorization symmetry in the lattice Boltzmann method. Physica A 389, 1530 (2010). The prepaper that this seemd to be based upon and had some more calculation details is
In [47]:
u = var("u",domain="real")
#T_0 = var("T_0",domain="real")
T_0 =var("T_0",domain="positive")
In [268]:
v = var("v",domain="real")
In [40]:
phi_v = sqrt( pi/(2*T_0))*exp( - (v-u)**2/(2*T_0))
In [48]:
integrate( phi_v * 1, v)
Out[48]:
In [54]:
integrate( phi_v * 1, (v,-oo,oo)).expand().simplify()
Out[54]:
In [57]:
print( integrate( phi_v * v**1, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**2, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**3, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**4, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**5, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**6, (v,-oo,oo)).expand().simplify() )
In [71]:
Q=5
q=Q//2
In [271]:
a_lst = var( ",".join(map(str,["a"+str(k+1) for k in range(q)])) )
v_i = var("v_i",domain="integer")
#v_lst = var( ",".join(map(str,["v"+str(k+1) for k in range(q)])) )
# variables named or labeled from 1,2,...q, Python counts from 0,1..q-1
In [76]:
k=var('k')
sum( a_lst[int(k)]*v_lst[k], k,0,q-1)
In [95]:
closurerel = sum( [a_lst[k]*v_i**(2*k+1) for k in range(q)] ) == v_i**Q # closure relation RHS
closurerelRHS = sum( [a_lst[k]*v_i**(2*k+1) for k in range(q)] ) # closure relation RHS
In [100]:
v_ipos = [1,3]
[closurerel.subs(v_i== v_ipos[k]) for k in range(len(v_ipos))]
#[closurerelRHS.subs((v_i, v_ipos[k])) for k in range(len(v_ipos))]
Out[100]:
In [101]:
solve( [closurerel.subs(v_i== v_ipos[k]) for k in range(len(v_ipos))], a_lst )
Out[101]:
In [105]:
solve( [closurerel.subs(v_i== v_ipos[k]) for k in range(len(v_ipos))], a_lst )[0]
Out[105]:
In [106]:
closurerel.subs( solve( [closurerel.subs(v_i== v_ipos[k]) for k in range(len(v_ipos))], a_lst )[0] )
Out[106]:
In [117]:
asolved=[aeqn.rhs() for aeqn in solve( [closurerel.subs(v_i== v_ipos[k]) for k in range(len(v_ipos))], a_lst )[0]]
print(asolved)
Try v_ipos=[1,2]
In [107]:
v_ipos = [1,2]
solve( [closurerel.subs(v_i== v_ipos[k]) for k in range(len(v_ipos))], a_lst )
Out[107]:
In [108]:
1.multifactorial(2)
Out[108]:
In [109]:
3.multifactorial(2)
Out[109]:
In [110]:
5.multifactorial(2)
Out[110]:
In [111]:
7.multifactorial(2)
Out[111]:
In [112]:
T_0=var('T_0',domain="positive")
In [122]:
sum( [asolved[k]*ZZ(2*k+1).multifactorial(2)*T_0**(k) for k in range(q)] )
Out[122]:
In [123]:
ZZ(Q).multifactorial(2)*T_0**((Q-1)/2)-sum( [asolved[k]*ZZ(2*k+1).multifactorial(2)*T_0**(k) for k in range(q)] )
Out[123]:
In [124]:
solve( ZZ(Q).multifactorial(2)*T_0**((Q-1)/2)-sum( [asolved[k]*ZZ(2*k+1).multifactorial(2)*T_0**(k) for k in range(q)] ) , T_0)
Out[124]:
In [125]:
solve( ZZ(Q).multifactorial(2)*T_0**((Q-1)/2)-sum( [asolved[k]*ZZ(2*k+1).multifactorial(2)*T_0**(k) for k in range(q)] )==0 , T_0)
Out[125]:
For
In [126]:
v_ipos = [1,3]
asolved=[aeqn.rhs() for aeqn in solve( [closurerel.subs(v_i== v_ipos[k]) for k in range(len(v_ipos))], a_lst )[0]]
In [127]:
solve( ZZ(Q).multifactorial(2)*T_0**((Q-1)/2)-sum( [asolved[k]*ZZ(2*k+1).multifactorial(2)*T_0**(k) for k in range(q)] )==0 , T_0)
Out[127]:
In [131]:
solve( ZZ(Q).multifactorial(2)*T_0**((Q-1)/2)-sum( [asolved[k]*ZZ(2*k+1).multifactorial(2)*T_0**(k) for k in range(q)] )==0 , T_0)[0].rhs().N >0
Out[131]:
In [132]:
solve( ZZ(Q).multifactorial(2)*T_0**((Q-1)/2)-sum( [asolved[k]*ZZ(2*k+1).multifactorial(2)*T_0**(k) for k in range(q)] )==0 , T_0)[1].rhs().N >0
Out[132]:
In [234]:
def make_closure_rel(Q, v_ipos=None,v_i=None):
""" make_closure_rel - make closure relation
returns, as Python tuple, the closure relation and reference temperature T_0, as (closure relation, T_0)
INPUT/ARGUMENT
==============
@type Q : positive integer
@param Q : Q is the total number of discrete velocities
"""
q=Q//2
a_lst = var( ",".join(map(str,["a"+str(k+1) for k in range(q)])) )
if v_i is None:
v_i = var("v_i",domain="integer")
if v_ipos is None:
v_ipos = range(1,q+1)
closurerel_gen = sum( [a_lst[k]*v_i**(2*k+1) for k in range(q)] ) == v_i**Q # closure relation, in general form
sysofpolys = [closurerel_gen.subs(v_i== v_ipos[k]) for k in range(len(v_ipos))] # q different nonzero values for velocities pluuged in
solved_as = solve( sysofpolys, a_lst ) # solved for a_{Q-2k},k=1...q, for given different velocities v_i
solved_as = solved_as[0]
closurerel = closurerel_gen.subs(solved_as)
# try:
# check_closurerel_unique = closurerel.rhs()-closurerel.lhs() # check the closure relation to see if it's nontrivial
# if check_closurerel_unique.N() == 0:
# print "check_closurerel_unique failed, equals 0, i.e. no unique solution for v_i's; lattice invalid"
# return closurerel
# except TypeError as err:
# print err
# print "so far closure relation is unique, proceed"
solved_as = [aeqn.rhs() for aeqn in solved_as] # [a_1,a_2,...a_q] as Python list
T_0=var('T_0',domain="positive")
matchingcond = -sum( [solved_as[k]*ZZ(2*k+1).multifactorial(2)*T_0**(k) for k in range(q)]) # \sum_{j=1}^q a_{Q-2k} b_{Q-2k} T_0^{(Q-2k-1)/2}
matchingcond += ZZ(Q).multifactorial(2)*T_0**((Q-1)/2) # b_QT_0^{(Q-1)/2}
# try:
# check_matchingcond = matchingcond.N()
# if check_matchingcond == 0:
# print "check_matchingcond failed, i.e. algebraic eqn. for ref. T_0 is 0"
# return closurerel
# except TypeError as err:
# print err
# print "We're actually safe, so far lattice exists; proceed"
T_0roots = solve( matchingcond ==0, T_0) # solve for the roots of T_0
T_0roots_passed = [] # T_0 roots s.t. T_0 >= 0
for Tsoln in T_0roots:
if (Tsoln.rhs().N >= 0):
T_0roots_passed.append( Tsoln.rhs())
return closurerel, T_0roots_passed
In [238]:
def make_closure_rel(Q, v_ipos=None,v_i=None):
""" make_closure_rel - make closure relation
returns, as Python tuple, the closure relation and reference temperature T_0, as (closure relation, T_0)
INPUT/ARGUMENT
==============
@type Q : positive integer
@param Q : Q is the total number of discrete velocities
"""
q=Q//2
a_lst = var( ",".join(map(str,["a"+str(k+1) for k in range(q)])) )
if v_i is None:
v_i = var("v_i",domain="integer")
if v_ipos is None:
v_ipos = range(1,q+1)
closurerel_gen = sum( [a_lst[k]*v_i**(2*k+1) for k in range(q)] ) == v_i**Q # closure relation, in general form
sysofpolys = [closurerel_gen.subs(v_i== v_ipos[k]) for k in range(len(v_ipos))] # q different nonzero values for velocities pluuged in
solved_as = solve( sysofpolys, a_lst ) # solved for a_{Q-2k},k=1...q, for given different velocities v_i
solved_as = solved_as[0]
closurerel = closurerel_gen.subs(solved_as)
solved_as = [aeqn.rhs() for aeqn in solved_as] # [a_1,a_2,...a_q] as Python list
T_0=var('T_0',domain="positive")
matchingcond = -sum( [solved_as[k]*ZZ(2*k+1).multifactorial(2)*T_0**(k) for k in range(q)]) # \sum_{j=1}^q a_{Q-2k} b_{Q-2k} T_0^{(Q-2k-1)/2}
matchingcond += ZZ(Q).multifactorial(2)*T_0**((Q-1)/2) # b_QT_0^{(Q-1)/2}
T_0roots = solve( matchingcond ==0, T_0) # solve for the roots of T_0
T_0roots_passed = [] # T_0 roots s.t. T_0 >= 0
for Tsoln in T_0roots:
if (Tsoln.rhs().N >= 0):
T_0roots_passed.append( Tsoln.rhs())
return closurerel, T_0roots_passed
In [229]:
make_closure_rel(5)
Out[229]:
In [230]:
make_closure_rel(5,[1,3])
Out[230]:
In [231]:
make_closure_rel(7)
Out[231]:
In [214]:
q=9//2
a_lst = var( ",".join(map(str,["a"+str(k+1) for k in range(q)])) )
v_i = var("v_i",domain="integer")
v_ipos = range(1,q+1)
closurerel_gen = sum( [a_lst[k]*v_i**(2*k+1) for k in range(q)] ) == v_i**Q
sysofpolys = [closurerel_gen.subs(v_i== v_ipos[k]) for k in range(len(v_ipos))]
solved_as_init = solve( sysofpolys, a_lst )
solved_as = solved_as_init[0]
closurerel = closurerel_gen.subs(solved_as)
solved_as = [aeqn.rhs() for aeqn in solved_as]
T_0=var('T_0',domain="positive")
matchingcond = -sum( [solved_as[k]*ZZ(2*k+1).multifactorial(2)*T_0**(k) for k in range(q)])
matchingcond += ZZ(Q).multifactorial(2)*T_0**((Q-1)/2)
matchingcond
Out[214]:
In [212]:
matchingcond.N() == 0
Out[212]:
In [227]:
( closurerel.lhs() - closurerel.rhs() ).N() == 0
Out[227]:
In [232]:
make_closure_rel(9)
In [239]:
make_closure_rel(9,[1,2,3,5])
Out[239]:
In [246]:
print( len(make_closure_rel(9,[1,2,3,5])[1]) )
print(make_closure_rel(9,[1,2,3,5])[1][0].N() )
make_closure_rel(9,[1,2,3,5])[1][1].N()
Out[246]:
In [253]:
print( make_closure_rel(11)[0] )
print( make_closure_rel(11)[1] )
find_root( make_closure_rel(11)[1][0], 0,2)
Out[253]:
In [254]:
make_closure_rel(13)
Out[254]:
In [259]:
find_root( make_closure_rel(13)[1][0], 0,10)
In [260]:
make_closure_rel(13,[1,2,3,4,5,7])
Out[260]:
In [261]:
find_root( make_closure_rel(13,[1,2,3,4,5,7])[1][0],0,2 ) # find a root between 0 and 2 for T_0
Out[261]:
In [262]:
make_closure_rel(15)
Out[262]:
In [264]:
find_root( make_closure_rel(15)[1][0],0,2 ) # find a root between 0 and 2 for T_0
Out[264]:
| $Q$ | $V$ | Closure | $T_0$ |
|---|---|---|---|
| 3 | $\lbrace 0, \pm 1 \rbrace$ | $v_{(i)}^3=v_{(i)}$ | $1/3$ |
| 5 | $\lbrace 0, \pm 1, \pm 3 \rbrace$ | $v_{(i)}^5=10v_{(i)}^3 - 9v_{(i)}$ | $1 \pm \sqrt{2/5}$ |
Table 1: 1-dim. Maxwell lattices with odd number of integer-valued velocities, $Q=3,5,7,9,11,13,15$
In [274]:
phi_v = sqrt( 1/(2*T_0*pi))*exp( - (v-u)**2/(2*T_0))
For $Q=3$
In [275]:
phi_v.subs(T_0==1/3).subs(v==v_i)
Out[275]:
In [276]:
phi_v.subs(T_0==1/3).subs(v==0)
Out[276]:
In [277]:
phi_v.subs(T_0==1/3).subs(v==1)
Out[277]:
In [278]:
phi_v.subs(T_0==1/3).subs(v==-1)
Out[278]:
In [282]:
( 0 + phi_v.subs(T_0==1/3).subs(v==1) + phi_v.subs(T_0==1/3).subs(v==-1)*(-1) ).expand().simplify()
Out[282]:
In [283]:
v_i = var("v_i",domain="integer")
make_closure_rel(5,[1,3])
Out[283]:
In [ ]:
In [ ]:
In [284]:
u = var("u",domain="real")
T_0 =var("T_0",domain="positive")
v = var("v",domain="real")
In [285]:
phi_v = sqrt( 1/(2*pi*T_0))*exp( - (v-u)**2/(2*T_0))
In [286]:
print( integrate( phi_v * v**0, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**1, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**2, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**3, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**4, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**5, (v,-oo,oo)).expand().simplify() )
print( integrate( phi_v * v**6, (v,-oo,oo)).expand().simplify() )
$d\equiv D=2$, 2-dim. case
In [288]:
u1,u2=var("u1,u2",domain="real")
v1,v2=var("v1,v2",domain="positive")
phi_v = 1/(2*pi*T_0)*exp( - ((v1-u1)**2+(v2-u2)**2)/(2*T_0))
In [298]:
print( integrate( integrate( phi_v * v1**0*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**1*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**0*v2**1, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**2*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**0*v2**2, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**1*v2**1, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**3*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**2*v2**1, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**1*v2**2, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**0*v2**3, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**4*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**3*v2**1, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**2*v2**2, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**1*v2**3, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**0*v2**4, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**5*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**4*v2**1, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**3*v2**2, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**2*v2**3, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**1*v2**4, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**0*v2**5, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
In [299]:
Q=9
W_lst = var(",".join(map(str,["W"+str(k) for k in range(Q)])) ) # Python list of weights W_i, i=0,1,...Q-1
In [300]:
c_0 = (0,0)
c_1,c_3 = (1,0),(-1,0); c_2,c_4=(0,1),(0,-1)
c_5,c_6,c_7,c_8 =(1,1),(-1,1),(1,-1),(-1,-1)
cs = [c_0,c_1,c_2,c_3,c_4,c_5,c_6,c_7,c_8]
In [302]:
sum(W_lst)
Out[302]:
In [303]:
sum( [W_lst[k]*cs[k][0] for k in range(Q)] )
Out[303]:
In [304]:
sum( [W_lst[k]*cs[k][1] for k in range(Q)] )
Out[304]:
In [306]:
print( sum( [W_lst[k]*cs[k][0]*cs[k][0] for k in range(Q)] ) )
print( sum( [W_lst[k]*cs[k][0]*cs[k][1] for k in range(Q)] ) )
print( sum( [W_lst[k]*cs[k][1]*cs[k][1] for k in range(Q)] ) )
In [307]:
print( sum( [W_lst[k]*cs[k][0]*cs[k][0]*cs[k][0] for k in range(Q)] ) )
print( sum( [W_lst[k]*cs[k][1]*cs[k][0]*cs[k][0] for k in range(Q)] ) )
print( sum( [W_lst[k]*cs[k][1]*cs[k][1]*cs[k][0] for k in range(Q)] ) )
# print( sum( [W_lst[k]*cs[k][1]*cs[k][1]*cs[k][1] for k in range(Q)] ) ) # the same
In [312]:
# print( sum( [W_lst[k]*cs[k][0]*cs[k][0]*cs[k][0]*cs[k][0] for k in range(Q)] ) ) # the same
# print( sum( [W_lst[k]*cs[k][1]*cs[k][0]*cs[k][0]*cs[k][0] for k in range(Q)] ) ) # the same
print( sum( [W_lst[k]*cs[k][0]*cs[k][0]*cs[k][1]*cs[k][1] for k in range(Q)] ) )
#print( sum( [W_lst[k]*cs[k][0]*cs[k][1]*cs[k][1]*cs[k][1] for k in range(Q)] ) ) # the same
#print( sum( [W_lst[k]*cs[k][1]*cs[k][1]*cs[k][1]*cs[k][1] for k in range(Q)] ) )
Is there a way to generate all the discrete velocities, $\mathbf{c}_i \equiv \mathbf{v}_i$?
I will use the mirror symmetry so we only need to consider the "positive values"
In [315]:
d=2 # D=2
v_ipos = [1,] # v_i positive values
dimidx = range(d) # indices for the components
In [317]:
[zeroidx for zeroidx in combinations(dimidx,1)]
Out[317]:
In [318]:
[_ for _ in combinations(v_ipos,1)]
Out[318]:
If $Q= $ total number of discrete vectors is given and known beforehand:
In [320]:
v_is =[[0 for _ in range(d)] for _ in range(Q)]
In [322]:
[_ for _ in combinations(dimidx,d)]
Out[322]:
In [328]:
[_ for _ in combinations(range(3),2)]
Out[328]:
In [324]:
[_ for _ in combinations_with_replacement(v_ipos,1) ]
Out[324]:
In [325]:
[_ for _ in combinations_with_replacement(v_ipos,2) ]
Out[325]:
In [335]:
d=3 # D=3
v_ipos = [1,2] # v_i positive values
dimidx = range(d) # indices for the components
In [336]:
c_d_nonzeros = [np.zeros(d) for _ in range( len([_ for _ in combinations(dimidx,2)]) )]
for nonzeros in combinations(dimidx,2): # d_nonzero = 0,1,..d-1
print( [_ for _ in combinations_with_replacement(v_ipos,2)] )
In [339]:
[_ for _ in combinations_with_replacement([0,1],1)]
Out[339]:
In [340]:
[_ for _ in combinations_with_replacement([0,1],2)]
Out[340]:
In [341]:
[_ for _ in combinations_with_replacement([0,1],3)]
Out[341]:
In [342]:
[_ for _ in combinations_with_replacement([0,1],4)]
Out[342]:
In [347]:
d=2 # D=2
v_ipos = [1,] # v_i positive values
dimidx = range(d) # indices for the components
In [343]:
unique_W_lst = []
In [345]:
if sum(W_lst) not in unique_W_lst:
unique_W_lst.append( sum(W_lst))
In [348]:
for alpha in combinations_with_replacement(dimidx,1):
W_to_lst = sum( [W_lst[k]*cs[k][alpha[0]] for k in range(Q)] )
if W_to_lst not in unique_W_lst:
unique_W_lst.append( W_to_lst )
print(unique_W_lst)
In [349]:
for alpha in combinations_with_replacement(dimidx,2):
W_to_lst = sum( [W_lst[k]*cs[k][alpha[0]]*cs[k][alpha[1]] for k in range(Q)] )
if W_to_lst not in unique_W_lst:
unique_W_lst.append( W_to_lst )
print(unique_W_lst)
In [350]:
for alpha in combinations_with_replacement(dimidx,3):
W_to_lst = sum( [W_lst[k]*cs[k][alpha[0]]*cs[k][alpha[1]]*cs[k][alpha[2]] for k in range(Q)] )
if W_to_lst not in unique_W_lst:
unique_W_lst.append( W_to_lst )
print(unique_W_lst)
In [351]:
for alpha in combinations_with_replacement(dimidx,4):
W_to_lst = sum( [W_lst[k]*cs[k][alpha[0]]*cs[k][alpha[1]]*cs[k][alpha[2]]*cs[k][alpha[3]] for k in range(Q)] )
if W_to_lst not in unique_W_lst:
unique_W_lst.append( W_to_lst )
print(unique_W_lst)
In [352]:
len(unique_W_lst)
Out[352]:
In [354]:
for alpha in combinations_with_replacement(dimidx,3):
todolst = [ W_lst[k]*prod( [cs[k][alpha[order]] for order in range(3)]) for k in range(Q)]
print sum(todolst)
In [376]:
def make_W_lin_sys(d,Q,cs, ORDER=4):
""" make_W_lin_sys - make the linear system of W_i's
NOTES
=====
Remember that W_lst corresponds directly with cs, i.e. W_i's with c_i \equiv v_i's (discrete velocity i)
"""
W_lst = var(",".join(map(str,["W"+str(k) for k in range(Q)])) ) # Python list of weights W_i, i=0,1,...Q-1
dimidx = range(d)
unique_W_lst = []
if sum(W_lst) not in unique_W_lst:
unique_W_lst.append( sum(W_lst))
for order in range(1,ORDER+1): # order = 1,2,..ORDER=4
for alpha in combinations_with_replacement(dimidx,order):
W_to_lst = [W_lst[k]*prod( [cs[k][alpha[orderidx]] for orderidx in range(order)]) for k in range(Q)]
W_to_lst = sum( W_to_lst )
if W_to_lst not in unique_W_lst:
unique_W_lst.append(W_to_lst)
return unique_W_lst
In [378]:
test_d2Q9_W_is = make_W_lin_sys(2,9,cs)
print(len(test_d2Q9_W_is))
In [380]:
print(len(test_d2Q9_W_is))
test_d2Q9_W_is
Out[380]:
In [359]:
range(1)
Out[359]:
Now consider the "right-hand side (RHS)", the "moments" or the velocity moments-over Maxwell distribution (what Wolf calls it)
In [364]:
u1,u2=var("u1,u2",domain="real")
v1,v2=var("v1,v2",domain="positive")
phi_v = 1/(2*pi*T_0)*exp( - ((v1-u1)**2+(v2-u2)**2)/(2*T_0))
In [365]:
print( integrate( integrate( phi_v * v1**0*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**1*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**0*v2**1, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**2*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**0*v2**2, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**1*v2**1, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**3*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**2*v2**1, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**1*v2**2, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**0*v2**3, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**4*v2**0, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**3*v2**1, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**2*v2**2, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**1*v2**3, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
print( integrate( integrate( phi_v * v1**0*v2**4, (v1,-oo,oo)),(v2,-oo,oo)).expand().simplify() )
In [366]:
for alpha in combinations_with_replacement(dimidx,1):
print alpha
In [367]:
for alpha in combinations_with_replacement(dimidx,2):
print alpha
In [368]:
for alpha in combinations_with_replacement(dimidx,3):
print alpha
In [496]:
def make_moment(d,alphas, u_vec=None,T_0=None):
""" make_moment - compute moment of Maxwell distribution
INPUTS/ARGUMENTS
================
@type alphas : Python iterable, most likely a Python tuple or Python list of length order, order of moment we want to calculate
@param alphas : alpha exponents
"""
v_i_vec = var(",".join(map(str,["v"+str(k+1) for k in range(d)])) )
if u_vec is None:
u_vec = var(",".join(map(str,["u"+str(k+1) for k in range(d)])) )
if T_0 is None:
T_0 =var("T_0",domain="positive")
if (d>=2):
norm_uminusv = sum( [ (v_i_vec[k]-u_vec[k])**2 for k in range(d)])
elif (d==1):
norm_uminusv = (v_i_vec - u_vec)**2
phi_v = 1/(2*pi*T_0)**(d/2)*exp( - norm_uminusv/(2*T_0))
if (d>=2):
tointegrate = phi_v*prod([v_i_vec[i] for i in alphas])
elif (d==1):
tointegrate = phi_v*prod([v_i_vec for i in alphas])
if (d>=2):
for v_i in v_i_vec:
tointegrate = integrate( tointegrate, (v_i,-oo,oo))
elif (d==1):
tointegrate = integrate( tointegrate, (v_i_vec,-oo,oo))
return tointegrate.expand().simplify()
In [484]:
test_M3 = make_moment(2,(0,1,1))
In [485]:
test_M3
Out[485]:
In [497]:
def make_W_lin_sys(d,Q,cs, ORDER=4,u_vec=None,T_0=None):
""" make_W_lin_sys - make the linear system of W_i's
NOTES
=====
Remember that W_lst corresponds directly with cs, i.e. W_i's with c_i \equiv v_i's (discrete velocity i)
"""
v_i_vec = var(",".join(map(str,["v"+str(k+1) for k in range(d)])) )
if u_vec is None:
u_vec = var(",".join(map(str,["u"+str(k+1) for k in range(d)])) )
if T_0 is None:
T_0 =var("T_0",domain="positive")
W_lst = var(",".join(map(str,["W"+str(k) for k in range(Q)])) ) # Python list of weights W_i, i=0,1,...Q-1
dimidx = range(d)
unique_W_lst = [] # Use this Python list to keep out redundancies, check for redundancies
unique_W_dict_lst = []
if sum(W_lst) not in unique_W_lst:
unique_W_lst.append( sum(W_lst))
rhs_out = make_moment(d,[],u_vec,T_0)
Weqn_out = {"lhs": sum(W_lst),"rhs":rhs_out}
unique_W_dict_lst.append(Weqn_out)
for order in range(1,ORDER+1): # order = 1,2,..ORDER=4
for alpha in combinations_with_replacement(dimidx,order):
W_to_lst = [W_lst[k]*prod( [cs[k][alpha[orderidx]] for orderidx in range(order)]) for k in range(Q)]
W_to_lst = sum( W_to_lst )
if W_to_lst not in unique_W_lst:
unique_W_lst.append(W_to_lst)
rhs_out = make_moment(d,alpha,u_vec,T_0)
Weqn_out = { "lhs":W_to_lst,"rhs":rhs_out}
unique_W_dict_lst.append( Weqn_out)
return unique_W_dict_lst, W_lst
In [498]:
test_W_sys, test_W_is = make_W_lin_sys(2,9,cs,4)
print(len(test_W_sys))
In [499]:
test_W_sys
Out[499]:
In [462]:
prod( [range(3)[i] for i in []] )
Out[462]:
The $W_i$ for directions with identical speeds are equal for reasons of symmetry.
In [500]:
#test_W_sys[0]["lhs"].subs(W2==W1).subs(W3==W1).subs(W4==W1).subs(W5==W2).subs(W6==W2).subs(W7==W2).subs(W8==W2)
for i in range(Q):
print( test_W_sys[i]["lhs"].subs(W2==W1).subs(W3==W1).subs(W4==W1).subs(W5==W2).subs(W6==W2).subs(W7==W2).subs(W8==W2) )
for i in range(Q):
print( test_W_sys[i]["rhs"].subs(u1==0).subs(u2==0))
In [501]:
test_W_sys_temp = []
for eqn in test_W_sys:
lhs =eqn["lhs"].subs(W2==W1).subs(W3==W1).subs(W4==W1).subs(W5==W2).subs(W6==W2).subs(W7==W2).subs(W8==W2)
rhs =eqn["rhs"].subs(u1==0).subs(u2==0)
test_W_sys_temp.append( lhs==rhs)
In [502]:
test_W_sys_temp
Out[502]:
In [503]:
solve( test_W_sys_temp, [W0,W1,W2])
Out[503]:
In [504]:
print( solve( test_W_sys_temp, [W0,W1,W2])[0][0].subs(T_0==1/3) )
print( solve( test_W_sys_temp, [W0,W1,W2])[0][1].subs(T_0==1/3) )
print( solve( test_W_sys_temp, [W0,W1,W2])[0][2].subs(T_0==1/3) )
So pp. 169 solution in Section 5.2.1of Wolf, and (same thing) on pp. 166, for the D2Q9 lattice (5.2.1), the $W_i$ are obtained.
In [505]:
make_closure_rel(5,[1,3])
Out[505]:
In [506]:
make_closure_rel(7)
Out[506]:
In [537]:
make_closure_rel(7)[1][0].N()
Out[537]:
In [415]:
make_closure_rel(9,[1,2,3,5])
Out[415]:
In [422]:
print( make_closure_rel(9,[1,2,3,5])[1][0].N() )
print( make_closure_rel(9,[1,2,3,5])[1][1].N() )
In [416]:
make_closure_rel(11)
Out[416]:
In [425]:
print( find_root( make_closure_rel(11)[1][0], 0,3) )
#print( find_root( make_closure_rel(11)[1][0], 3,6) )
In [442]:
make_closure_rel(13,[1,2,3,4,5,7])
Out[442]:
In [446]:
print( find_root( make_closure_rel(13,[1,2,3,4,5,7])[1][0], 0,3) )
#print( find_root( make_closure_rel(13,[1,2,3,4,5,7])[1][0], 3,6) )
In [418]:
make_closure_rel(15)
Out[418]:
In [436]:
print( find_root( make_closure_rel(15)[1][0], 0,2) )
In [447]:
v_0 = (0,)
v_1 = (1,)
v_m1=(-1,)
V = [v_0,v_1,v_m1]
In [508]:
#make_W_lin_sys(1,3,V, ORDER=2,u_vec=None,T_0=None)
D1Q3_W_sys, D1Q3_W_is = make_W_lin_sys(1,3,V,ORDER=2)
In [510]:
def W_sys_solver( sys_from_make_W_lin_sys, W_is):
eqns_to_solve = []
for eqn in sys_from_make_W_lin_sys:
eqns_to_solve.append( eqn['lhs'] == eqn['rhs'] )
return solve( eqns_to_solve,W_is)
In [511]:
W_sys_solver( D1Q3_W_sys, D1Q3_W_is)
Out[511]:
In [520]:
print( W_sys_solver( D1Q3_W_sys, D1Q3_W_is)[0][0].subs(u1==0).subs(T_0==1/3) )
print( W_sys_solver( D1Q3_W_sys, D1Q3_W_is)[0][1].subs(u1==0).subs(T_0==1/3) )
print( W_sys_solver( D1Q3_W_sys, D1Q3_W_is)[0][2].subs(u1==0).subs(T_0==1/3) )
In [457]:
len(v_i_vec)
Out[457]:
In [512]:
V = [(0,),(1,),(-1,),(3,),(-3,)]
D1Q5_W_sys, D1Q5_W_is = make_W_lin_sys(1,5,V,ORDER=4)
In [521]:
W_sys_solver( D1Q5_W_sys, D1Q5_W_is)
Out[521]:
In [531]:
print( W_sys_solver( D1Q5_W_sys, D1Q5_W_is)[0][0].subs(u1==0).subs(T_0==-1/5*sqrt(10) + 1).rhs().expand().simplify() )
print( W_sys_solver( D1Q5_W_sys, D1Q5_W_is)[0][1].subs(u1==0).subs(T_0==-1/5*sqrt(10) + 1).rhs().expand().simplify() )
print( W_sys_solver( D1Q5_W_sys, D1Q5_W_is)[0][2].subs(u1==0).subs(T_0==-1/5*sqrt(10) + 1).rhs().expand().simplify() )
print( W_sys_solver( D1Q5_W_sys, D1Q5_W_is)[0][3].subs(u1==0).subs(T_0==-1/5*sqrt(10) + 1).rhs().expand().simplify() )
print( W_sys_solver( D1Q5_W_sys, D1Q5_W_is)[0][4].subs(u1==0).subs(T_0==-1/5*sqrt(10) + 1).rhs().expand().simplify() )
In [532]:
print( W_sys_solver( D1Q5_W_sys, D1Q5_W_is)[0][0].subs(u1==0).subs(T_0==1/5*sqrt(10) + 1).rhs().expand().simplify() )
print( W_sys_solver( D1Q5_W_sys, D1Q5_W_is)[0][1].subs(u1==0).subs(T_0==1/5*sqrt(10) + 1).rhs().expand().simplify() )
print( W_sys_solver( D1Q5_W_sys, D1Q5_W_is)[0][2].subs(u1==0).subs(T_0==1/5*sqrt(10) + 1).rhs().expand().simplify() )
print( W_sys_solver( D1Q5_W_sys, D1Q5_W_is)[0][3].subs(u1==0).subs(T_0==1/5*sqrt(10) + 1).rhs().expand().simplify() )
print( W_sys_solver( D1Q5_W_sys, D1Q5_W_is)[0][4].subs(u1==0).subs(T_0==1/5*sqrt(10) + 1).rhs().expand().simplify() )
In [533]:
V = [(0,),(1,),(-1,),(2,),(-2,),(3,),(-3,)]
D1Q7_W_sys, D1Q7_W_is = make_W_lin_sys(1,7,V,ORDER=6)
In [535]:
D1Q7_W_solved= W_sys_solver( D1Q7_W_sys, D1Q7_W_is)
In [542]:
print( D1Q7_W_solved[0][0].subs(u1==0).subs(T_0==make_closure_rel(7)[1][0]).rhs().N() )
print( D1Q7_W_solved[0][1].subs(u1==0).subs(T_0==make_closure_rel(7)[1][0]).rhs().N() )
print( D1Q7_W_solved[0][2].subs(u1==0).subs(T_0==make_closure_rel(7)[1][0]).rhs().N() )
print( D1Q7_W_solved[0][3].subs(u1==0).subs(T_0==make_closure_rel(7)[1][0]).rhs().N() )
print( D1Q7_W_solved[0][4].subs(u1==0).subs(T_0==make_closure_rel(7)[1][0]).rhs().N() )
print( D1Q7_W_solved[0][5].subs(u1==0).subs(T_0==make_closure_rel(7)[1][0]).rhs().N() )
print( D1Q7_W_solved[0][6].subs(u1==0).subs(T_0==make_closure_rel(7)[1][0]).rhs().N() )
In [555]:
for i in range(7):
print( D1Q7_W_solved[0][i].subs(u1==0) )
Note that Eqns.(7-10) of Physical Review E 93, 063302 (2016), "Entropic Lattice Boltzmann Model for Gas..." is computed!
In [543]:
V = [(0,),(1,),(-1,),(2,),(-2,),(3,),(-3,),(5,),(-5,)]
D1Q9_W_sys, D1Q9_W_is = make_W_lin_sys(1,9,V,ORDER=8)
In [544]:
D1Q9_W_solved= W_sys_solver( D1Q9_W_sys, D1Q9_W_is)
In [547]:
print( D1Q9_W_solved[0][0].subs(u1==0).subs(T_0==0.756080852594268 ) )
print( D1Q9_W_solved[0][1].subs(u1==0).subs(T_0==0.756080852594268 ) )
print( D1Q9_W_solved[0][2].subs(u1==0).subs(T_0==0.756080852594268 ) )
print( D1Q9_W_solved[0][3].subs(u1==0).subs(T_0==0.756080852594268 ) )
print( D1Q9_W_solved[0][4].subs(u1==0).subs(T_0==0.756080852594268 ) )
print( D1Q9_W_solved[0][5].subs(u1==0).subs(T_0==0.756080852594268 ) )
print( D1Q9_W_solved[0][6].subs(u1==0).subs(T_0==0.756080852594268 ) )
print( D1Q9_W_solved[0][7].subs(u1==0).subs(T_0==0.756080852594268 ) )
print( D1Q9_W_solved[0][8].subs(u1==0).subs(T_0==0.756080852594268 ) )
In [548]:
V = [(0,),(1,),(-1,),(2,),(-2,),(3,),(-3,),(4,),(-4,),(5,),(-5,)]
D1Q11_W_sys, D1Q11_W_is = make_W_lin_sys(1,11,V,ORDER=10)
In [549]:
D1Q11_W_solved= W_sys_solver( D1Q11_W_sys, D1Q11_W_is)
In [550]:
for i in range(11):
print( D1Q11_W_solved[0][i].subs(u1==0).subs(T_0==1.06279357749 ) )
In [551]:
V = [(0,),(1,),(-1,),(2,),(-2,),(3,),(-3,),(4,),(-4,),(5,),(-5,),(7,),(-7,)]
D1Q13_W_sys, D1Q13_W_is = make_W_lin_sys(1,13,V,ORDER=12)
In [552]:
D1Q13_W_solved= W_sys_solver( D1Q13_W_sys, D1Q13_W_is)
In [554]:
for i in range(13):
print( D1Q13_W_solved[0][i].subs(u1==0).subs(T_0==2.67428803233 ) )
In [ ]: