In this tutorial, we test a very classical lattice Boltzmann scheme \(\DdQq{1}{3}\) on the heat equation.
The problem reads
where \(\mu\) is a constant scalar.
[8]:
from __future__ import print_function, division
from six.moves import range
%matplotlib inline
The numerical simulation of this equation by a lattice Boltzmann scheme consists in the approximatation of the solution on discret points of \((0,1)\) at discret instants.
To simulate this system of equations, we use the \(\DdQq{1}{3}\) scheme given by
three velocities \(v_0=0\), \(v_1=1\), and \(v_2=-1\), with associated distribution functions \(\fk{0}\), \(\fk{1}\), and \(\fk{2}\),
a space step \(\dx\) and a time step \(\dt\), the ration \(\lambda=\dx/\dt\) is called the scheme velocity,
three moments
and their equilibrium values \(\mke{0}\), \(\mke{1}\), and \(\mke{2}\).
two relaxation parameters \(s_1\) and \(s_2\) lying in \([0,2]\).
In order to use the formalism of the package pylbm, we introduce the three polynomials that define the moments: \(P_0 = 1\), \(P_1=X\), and \(P_2=X^2/2\), such that
The transformation \((\fk{0}, \fk{1}, \fk{2})\mapsto(\mk{0},\mk{1}, \mk{2})\) is invertible if, and only if, the polynomials \((P_0,P_1,P_2)\) is a free set over the stencil of velocities.
The lattice Boltzmann method consists to compute the distribution functions \(\fk{0}\), \(\fk{1}\), and \(\fk{2}\) in each point of the lattice \(x\) and at each time \(t^n=n\dt\). A step of the scheme can be read as a splitting between the relaxation phase and the transport phase:
relaxation:
m2f:
transport:
f2m:
The moment of order \(0\), \(\mk{0}\), being conserved during the relaxation phase, a diffusive scaling \(\dt=\dx^2\), yields to the following equivalent equation
if \(\mke{1}=0\). In order to be consistent with the heat equation, the following choice is done:
pylbm uses Python dictionary to describe the simulation. In the following, we will build this dictionary step by step.
In pylbm, the geometry is defined by a box and a label for the boundaries.
[9]:
import pylbm
import numpy as np
xmin, xmax = 0., 1.
dico_geom = {
'box': {'x': [xmin, xmax], 'label':0},
}
geom = pylbm.Geometry(dico_geom)
print(geom)
geom.visualize(viewlabel=True)
Geometry informations
spatial dimension: 1
bounds of the box:
[[0. 1.]]
pylbm provides a class stencil that is used to define the discret velocities of the scheme. In this example, the stencil is composed by the velocities \(v_0=0\), \(v_1=1\) and \(v_2=-1\) numbered by \([0,1,2]\).
[10]:
dico_sten = {
'dim':1,
'schemes':[{'velocities':list(range(3))}],
}
sten = pylbm.Stencil(dico_sten)
print(sten)
sten.visualize()
Stencil informations
* spatial dimension: 1
* maximal velocity in each direction: [1]
* minimal velocity in each direction: [-1]
* Informations for each elementary stencil:
stencil 0
- number of velocities: 3
- velocities: (0: 0), (1: 1), (2: -1),
In order to build the domain of the simulation, the dictionary should contain the space step \(\dx\) and the stencils of the velocities (one for each scheme).
We construct a domain with \(N=10\) points in space.
[11]:
N = 10
dx = (xmax-xmin)/N
dico_dom = {
'box': {'x': [xmin, xmax], 'label':0},
'space_step':dx,
'schemes':[
{
'velocities':list(range(3)),
}
],
}
dom = pylbm.Domain(dico_dom)
print(dom)
dom.visualize()
Domain informations
spatial dimension: 1
space step: dx= 1.000e-01
In pylbm, a simulation can be performed by using several coupled schemes. In this example, a single scheme is used and defined through a list of one single dictionary. This dictionary should contain:
‘velocities’: a list of the velocities
‘conserved_moments’: a list of the conserved moments as sympy variables
‘polynomials’: a list of the polynomials that define the moments
‘equilibrium’: a list of the equilibrium value of all the moments
‘relaxation_parameters’: a list of the relaxation parameters (\(0\) for the conserved moments)
‘init’: a dictionary to initialize the conserved moments
(see the documentation for more details)
The scheme velocity could be taken to \(1/\dx\) and the inital value of \(u\) to
[12]:
import sympy as sp
def solution(x, t):
return np.sin(np.pi*x)*np.exp(-np.pi**2*mu*t)
# parameters
mu = 1.
la = 1./dx
s1 = 2./(1+2*mu)
s2 = 1.
u, X = sp.symbols('u, X')
dico_sch = {
'dim':1,
'scheme_velocity':la,
'schemes':[
{
'velocities':list(range(3)),
'conserved_moments':u,
'polynomials':[1, X, X**2/2],
'equilibrium':[u, 0., .5*u],
'relaxation_parameters':[0., s1, s2],
'init':{u:(solution, (0.,))},
}
],
}
sch = pylbm.Scheme(dico_sch)
print(sch)
[0] WARNING pylbm.scheme in function __init__ line 194
The value 'space_step' is not given or wrong.
The scheme takes default value: dx = 1.
Scheme informations
spatial dimension: dim=1
number of schemes: nscheme=1
number of velocities:
Stencil.nv[0]=3
velocities value:
v[0] = (0: 0), (1: 1), (2: -1),
polynomials:
P[0] = 1, X, X**2/2,
equilibria:
EQ[0] = u, 0.0, 0.5*u,
relaxation parameters:
s[0] = 0.0, 0.666666666666667, 1.00000000000000,
moments matrices
M = Matrix([[1, 1, 1], [0, 10.0000000000000, -10.0000000000000], [0, 50.0000000000000, 50.0000000000000]])
M^(-1) = Matrix([[1.00000000000000, 0, -0.0200000000000000], [0, 0.0500000000000000, 0.0100000000000000], [0, -0.0500000000000000, 0.0100000000000000]])
A simulation is built by defining a correct dictionary.
We combine the previous dictionaries to build a simulation. In order to impose the homogeneous Dirichlet conditions in \(x=0\) and \(x=1\), the dictionary should contain the key ‘boundary_conditions’ (we use pylbm.bc.Anti_bounce_back function).
[13]:
dico = {
'box':{'x':[xmin, xmax], 'label':0},
'space_step':dx,
'scheme_velocity':la,
'schemes':[
{
'velocities':list(range(3)),
'conserved_moments':u,
'polynomials':[1, X, X**2/2],
'equilibrium':[u, 0., .5*u],
'relaxation_parameters':[0., s1, s2],
'init':{u:(solution,(0.,))},
}
],
'boundary_conditions':{
0:{'method':{0:pylbm.bc.anti_bounce_back,}, 'value':None},
},
'generator': 'numpy'
}
sol = pylbm.Simulation(dico)
print(sol)
Simulation informations:
Domain informations
spatial dimension: 1
space step: dx= 1.000e-01
Scheme informations
spatial dimension: dim=1
number of schemes: nscheme=1
number of velocities:
Stencil.nv[0]=3
velocities value:
v[0] = (0: 0), (1: 1), (2: -1),
polynomials:
P[0] = 1, X, X**2/2,
equilibria:
EQ[0] = u, 0.0, 0.5*u,
relaxation parameters:
s[0] = 0.0, 0.666666666666667, 1.00000000000000,
moments matrices
M = Matrix([[1, 1, 1], [0, 10.0000000000000, -10.0000000000000], [0, 50.0000000000000, 50.0000000000000]])
M^(-1) = Matrix([[1.00000000000000, 0, -0.0200000000000000], [0, 0.0500000000000000, 0.0100000000000000], [0, -0.0500000000000000, 0.0100000000000000]])
Once the simulation is initialized, one time step can be performed by using the function one_time_step.
We compute the solution of the heat equation at \(t=0.1\). And, on the same graphic, we plot the initial condition, the exact solution and the numerical solution.
[14]:
import numpy as np
import sympy as sp
import pylab as plt
import pylbm
u, X, LA = sp.symbols('u, X, LA')
def solution(x, t):
return np.sin(np.pi*x)*np.exp(-np.pi**2*mu*t)
xmin, xmax = 0., 1.
N = 128
mu = 1.
Tf = .1
dx = (xmax-xmin)/N # spatial step
la = 1./dx
s1 = 2./(1+2*mu)
s2 = 1.
dico = {
'box':{'x':[xmin,xmax], 'label':0},
'space_step':dx,
'scheme_velocity':la,
'schemes':[
{
'velocities':list(range(3)),
'conserved_moments':u,
'polynomials':[1, X/LA, X**2/(2*LA**2)],
'equilibrium':[u, 0., .5*u],
'relaxation_parameters':[0., s1, s2],
'init':{u:(solution,(0.,))},
}
],
'boundary_conditions':{
0:{'method':{0:pylbm.bc.anti_bounce_back,}, 'value':None},
},
'parameters':{LA: la},
'generator': 'numpy'
}
sol = pylbm.Simulation(dico)
x = sol.domain.x
y = sol.m[u]
plt.figure(1)
plt.plot(x, y,'k', label='initial')
while sol.t < 0.1:
sol.one_time_step()
plt.plot(x, sol.m[u],'b', label=r'$D_1Q_3$')
plt.plot(x, solution(x, sol.t),'r', label='exact')
plt.title('Heat equation t={0:5.3f}'.format(sol.t))
plt.legend()
[14]:
<matplotlib.legend.Legend at 0x7f980defffd0>
[ ]: