Lattice Boltzmann Method



In [1]:

    
import sympy



In [2]:

    
from sympy import exp, integrate, pi, sqrt, Symbol, symbols, oo



In [3]:

    
from sympy import Abs, Q, periodic_argument, polar_lift, refine



In [4]:

    
from sympy import cos



In [5]:

    
from sympy import Rational as Rat

$d2q9$

cf.Jonas Tölke. Implementation of a Lattice Boltzmann kernel using the Compute Unified Device Architecture developed by nVIDIA. Comput. Visual Sci. DOI 10.1007/s00791-008-0120-2

Affine Spaces

Lattices, that are "sufficiently" Galilean invariant, through non-perturbative algebraic theory

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

Maxwell Lattices in 1-dim.

Maxwell's (M) moment relations



In [19]:

    
#u = Symbol("u",assume="real")
u = Symbol("u",real=True)
#T_0 =Symbol("T_0",assume="positive")
T_0 =Symbol("T_0",real=True,positive=True)



In [30]:

    
#v = Symbol("v",assume="real")
v = Symbol("v",real=True)



In [37]:

    
#phi_v = sqrt( pi/(Rat(2)*T_0))*exp( - (v-u)**2/(Rat(2)*T_0))
phi_v = sqrt( pi/(2*T_0))*exp( - (v-u)**2/(2*T_0))



In [38]:

    
integrate(phi_v,v)









    Out[38]:





pi*sqrt(T_0)*sqrt(1/T_0)*erf(sqrt(2)*(-u + v)/(2*sqrt(T_0)))/2



In [39]:

    
integrate(phi_v,(v,-oo,oo))









    Out[39]:





Piecewise((pi*(-erfc(sqrt(2)*u/(2*sqrt(T_0))) + 2)/2 + pi*erfc(sqrt(2)*u/(2*sqrt(T_0)))/2, And(Or(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(exp_polar(2*I*pi)*polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(exp_polar(2*I*pi)*polar_lift(u)**2/polar_lift(T_0)**2, oo)) <= pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) <= pi/2, Abs(periodic_argument(exp_polar(2*I*pi)*polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi)), Or(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(polar_lift(u)**2/polar_lift(T_0)**2, oo)) <= pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) <= pi/2, Abs(periodic_argument(polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi)))), (Integral(sqrt(2)*sqrt(pi)*sqrt(1/T_0)*exp(-(-u + v)**2/(2*T_0))/2, (v, -oo, oo)), True))



In [40]:

    
(integrate(phi_v,v).subs(u,oo)- integrate(phi_v,v).subs(u,-oo)).expand()









    Out[40]:





pi*sqrt(T_0)*sqrt(1/T_0)*erf(sqrt(2)*v/(2*sqrt(T_0)) - oo/sqrt(T_0))/2 - pi*sqrt(T_0)*sqrt(1/T_0)*erf(sqrt(2)*v/(2*sqrt(T_0)) + oo/sqrt(T_0))/2



In [41]:

    
integrate(phi_v,(v,-oo,oo),conds='none')









    Out[41]:





pi*(-erfc(sqrt(2)*u/(2*sqrt(T_0))) + 2)/2 + pi*erfc(sqrt(2)*u/(2*sqrt(T_0)))/2



In [42]:

    
integrate(phi_v,(v,-oo,oo),conds='separate')









    Out[42]:





(pi*(-erfc(sqrt(2)*u/(2*sqrt(T_0))) + 2)/2 + pi*erfc(sqrt(2)*u/(2*sqrt(T_0)))/2,
 And(Or(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(exp_polar(2*I*pi)*polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(exp_polar(2*I*pi)*polar_lift(u)**2/polar_lift(T_0)**2, oo)) <= pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) <= pi/2, Abs(periodic_argument(exp_polar(2*I*pi)*polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi)), Or(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(polar_lift(u)**2/polar_lift(T_0)**2, oo)) <= pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) <= pi/2, Abs(periodic_argument(polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi))))

cf. http://stackoverflow.com/questions/16599325/simplify-conditional-integrals-in-sympy



In [51]:

    
refine(integrate(phi_v,(v,-oo,oo)), Q.is_true(Abs(periodic_argument(1/polar_lift(sqrt(T_0))**2, oo)) <= pi/2))









    Out[51]:





Piecewise((pi*(-erfc(sqrt(2)*u/(2*sqrt(T_0))) + 2)/2 + pi*erfc(sqrt(2)*u/(2*sqrt(T_0)))/2, And(Or(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(exp_polar(2*I*pi)*polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(exp_polar(2*I*pi)*polar_lift(u)**2/polar_lift(T_0)**2, oo)) <= pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) <= pi/2, Abs(periodic_argument(exp_polar(2*I*pi)*polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi)), Or(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) < pi/2, Abs(periodic_argument(polar_lift(u)**2/polar_lift(T_0)**2, oo)) <= pi), And(Abs(periodic_argument(1/polar_lift(T_0), oo)) <= pi/2, Abs(periodic_argument(polar_lift(u)**2/polar_lift(T_0)**2, oo)) < pi)))), (Integral(sqrt(2)*sqrt(pi)*sqrt(1/T_0)*exp(-(-u + v)**2/(2*T_0))/2, (v, -oo, oo)), True))



In [ ]: