In [ ]:
from __future__ import print_function, division
from six.moves import range
In [1]:
%matplotlib inline
In this tutorial, we test a very classical lattice Boltzmann scheme \(\DdQq{2}{5}\) on the heat equation.
The problem reads
where \(\mu\) is a constant scalar.
The numerical simulation of this equation by a lattice Boltzmann scheme consists in the approximatation of the solution on discret points of \((0,1)^2\) at discret instants.
To simulate this system of equations, we use the \(\DdQq{2}{5}\) scheme given by
five velocities \(v_0=(0,0)\), \(v_1=(1,0)\), \(v_2=(0,1)\), \(v_3=(-1,0)\), and \(v_4=(0,-1)\) with associated distribution functions \(\fk{i}\), \(0\leq i\leq 4\),
a space step \(\dx\) and a time step \(\dt\), the ration \(\lambda=\dx/\dt\) is called the scheme velocity,
five moments
and their equilibrium values \(\mke{k}\), \(0\leq k\leq 4\). * two relaxation parameters \(s_1\) and \(s_2\) lying in \([0,2]\) (\(s_1\) for the odd moments and \(s_2\) for the odd ones).
In order to use the formalism of the package pyLBM, we introduce the five polynomials that define the moments: \(P_0 = 1\), \(P_1=X\), \(P_2=Y\), \(P_3=(X^2+Y^2)/2\), and \(P_4=(X^2-Y^2)/2\), such that
The transformation \((\fk{0}, \fk{1}, \fk{2}, \fk{3}, \fk{4})\mapsto(\mk{0},\mk{1}, \mk{2}, \mk{3}, \mk{4})\) is invertible if, and only if, the polynomials \((P_0,P_1,P_2,P_3,P_4)\) is a free set over the stencil of velocities.
The lattice Boltzmann method consists to compute the distribution functions \(\fk{i}\), \(0\leq i\leq 4\) 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. We define here a square \((0, 1)^2\).
In [2]:
import pyLBM
import numpy as np
import pylab as plt
xmin, xmax, ymin, ymax = 0., 1., 0., 1.
dico_geom = {
'box': {'x': [xmin, xmax], 'y':[ymin, ymax], 'label':0},
}
geom = pyLBM.Geometry(dico_geom)
print(geom)
geom.visualize(viewlabel=True)
Geometry informations
spatial dimension: 2
bounds of the box:
[[ 0. 1.]
[ 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,0)\), \(v_1=(1,0)\), \(v_2=(-1,0)\), \(v_3=(0,1)\), and \(v_4=(0,-1)\) numbered by \([0,1,2,3,4]\).
In [3]:
dico_sten = {
'dim':2,
'schemes':[{'velocities':list(range(5))}],
}
sten = pyLBM.Stencil(dico_sten)
print(sten)
sten.visualize()
Stencil informations
* spatial dimension: 2
* maximal velocity in each direction: [1 1]
* minimal velocity in each direction: [-1 -1]
* Informations for each elementary stencil:
stencil 0
- number of velocities: 5
- velocities: (0: 0, 0), (1: 1, 0), (2: 0, 1), (3: -1, 0), (4: 0, -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.
In [4]:
N = 10
dx = (xmax-xmin)/N
dico_dom = {
'box': {'x': [xmin, xmax], 'y':[ymin, ymax], 'label':0},
'space_step':dx,
'schemes':[
{
'velocities':list(range(5)),
}
],
}
dom = pyLBM.Domain(dico_dom)
print(dom)
dom.visualize(view_distance=True)
Domain informations
spatial dimension: 2
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:
(see the documentation for more details)
The scheme velocity could be taken to \(1/\dx\) and the inital value of \(u\) to
In [5]:
import sympy as sp
def solution(x, y, t):
return np.sin(np.pi*x)*np.sin(np.pi*y)*np.exp(-2*np.pi**2*mu*t)
# parameters
mu = 1.
la = 1./dx
s1 = 2./(1+4*mu)
s2 = 1.
u, X, Y = sp.symbols('u, X, Y')
dico_sch = {
'dim':2,
'scheme_velocity':la,
'schemes':[
{
'velocities':list(range(5)),
'conserved_moments':u,
'polynomials':[1, X, Y, (X**2+Y**2)/2, (X**2-Y**2)/2],
'equilibrium':[u, 0., 0., .5*u, 0.],
'relaxation_parameters':[0., s1, s1, s2, s2],
'init':{u:(solution, (0.,))},
}
],
}
sch = pyLBM.Scheme(dico_sch)
print(sch)
[0] WARNING pyLBM.scheme in function __init__ line 229
The value 'space_step' is not given or wrong.
The scheme takes default value: dx = 1.
WARNING:pyLBM.scheme:The value 'space_step' is not given or wrong.
The scheme takes default value: dx = 1.
Scheme informations
spatial dimension: dim=2
number of schemes: nscheme=1
number of velocities:
Stencil.nv[0]=5
velocities value:
v[0]=(0: 0, 0), (1: 1, 0), (2: 0, 1), (3: -1, 0), (4: 0, -1),
polynomials:
P[0]=Matrix([[1], [X], [Y], [X**2/2 + Y**2/2], [X**2/2 - Y**2/2]])
equilibria:
EQ[0]=Matrix([[u], [0.0], [0.0], [0.5*u], [0.0]])
relaxation parameters:
s[0]=[0.0, 0.4, 0.4, 1.0, 1.0]
moments matrices
M = [Matrix([
[1, 1, 1, 1, 1],
[0, 1, 0, -1, 0],
[0, 0, 1, 0, -1],
[0, 1/2, 1/2, 1/2, 1/2],
[0, 1/2, -1/2, 1/2, -1/2]])]
invM = [Matrix([
[1, 0, 0, -2, 0],
[0, 1/2, 0, 1/2, 1/2],
[0, 0, 1/2, 1/2, -1/2],
[0, -1/2, 0, 1/2, 1/2],
[0, 0, -1/2, 1/2, -1/2]])]
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\), \(x=1\), \(y=0\), and \(y=1\), the dictionary should contain the key ‘boundary_conditions’ (we use pyLBM.bc.Anti_bounce_back function).
In [6]:
dico = {
'box':{'x':[xmin, xmax], 'y':[ymin, ymax], 'label':0},
'space_step':dx,
'scheme_velocity':la,
'schemes':[
{
'velocities':list(range(5)),
'conserved_moments':u,
'polynomials':[1, X, Y, (X**2+Y**2)/2, (X**2-Y**2)/2],
'equilibrium':[u, 0., 0., .5*u, 0.],
'relaxation_parameters':[0., s1, s1, s2, s2],
'init':{u:(solution,(0.,))},
}
],
'boundary_conditions':{
0:{'method':{0:pyLBM.bc.anti_bounce_back,}, 'value':None},
},
}
sol = pyLBM.Simulation(dico)
print(sol)
Simulation informations:
Domain informations
spatial dimension: 2
space step: dx= 1.000e-01
Scheme informations
spatial dimension: dim=2
number of schemes: nscheme=1
number of velocities:
Stencil.nv[0]=5
velocities value:
v[0]=(0: 0, 0), (1: 1, 0), (2: 0, 1), (3: -1, 0), (4: 0, -1),
polynomials:
P[0]=Matrix([[1], [X], [Y], [X**2/2 + Y**2/2], [X**2/2 - Y**2/2]])
equilibria:
EQ[0]=Matrix([[u], [0.0], [0.0], [0.5*u], [0.0]])
relaxation parameters:
s[0]=[0.0, 0.4, 0.4, 1.0, 1.0]
moments matrices
M = [Matrix([
[1, 1, 1, 1, 1],
[0, 1, 0, -1, 0],
[0, 0, 1, 0, -1],
[0, 1/2, 1/2, 1/2, 1/2],
[0, 1/2, -1/2, 1/2, -1/2]])]
invM = [Matrix([
[1, 0, 0, -2, 0],
[0, 1/2, 0, 1/2, 1/2],
[0, 0, 1/2, 1/2, -1/2],
[0, -1/2, 0, 1/2, 1/2],
[0, 0, -1/2, 1/2, -1/2]])]
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\). On the same graphic, we plot the initial condition, the exact solution and the numerical solution.
In [9]:
import numpy as np
import sympy as sp
import pylab as plt
%matplotlib inline
from mpl_toolkits.axes_grid1 import make_axes_locatable
import pyLBM
u, X, Y = sp.symbols('u, X, Y')
def solution(x, y, t, k, l):
return np.sin(k*np.pi*x)*np.sin(l*np.pi*y)*np.exp(-(k**2+l**2)*np.pi**2*mu*t)
def plot(i, j, z, title):
im = axarr[i,j].imshow(z)
divider = make_axes_locatable(axarr[i, j])
cax = divider.append_axes("right", size="20%", pad=0.05)
cbar = plt.colorbar(im, cax=cax, format='%6.0e')
axarr[i, j].xaxis.set_visible(False)
axarr[i, j].yaxis.set_visible(False)
axarr[i, j].set_title(title)
# parameters
xmin, xmax, ymin, ymax = 0., 1., 0., 1.
N = 128
mu = 1.
Tf = .1
dx = (xmax-xmin)/N # spatial step
la = 1./dx
s1 = 2./(1+4*mu)
s2 = 1.
k, l = 1, 1 # number of the wave
dico = {
'box':{'x':[xmin, xmax], 'y':[ymin, ymax], 'label':0},
'space_step':dx,
'scheme_velocity':la,
'schemes':[
{
'velocities':list(range(5)),
'conserved_moments':u,
'polynomials':[1, X, Y, (X**2+Y**2)/2, (X**2-Y**2)/2],
'equilibrium':[u, 0., 0., .5*u, 0.],
'relaxation_parameters':[0., s1, s1, s2, s2],
'init':{u:(solution,(0.,k,l))},
}
],
'boundary_conditions':{
0:{'method':{0:pyLBM.bc.anti_bounce_back,}, 'value':None},
},
'generator':pyLBM.generator.CythonGenerator,
}
sol = pyLBM.Simulation(dico)
x = sol.domain.x
y = sol.domain.y
f, axarr = plt.subplots(2, 2)
f.suptitle('Heat equation', fontsize=20)
plot(0, 0, sol.m[u].copy(), 'initial')
while sol.t < Tf:
sol.one_time_step()
sol.f2m()
z = sol.m[u]
ze = solution(x[:,np.newaxis], y[np.newaxis,:], sol.t, k, l)
plot(1, 0, z, 'final')
plot(0, 1, ze, 'exact')
plot(1, 1, z-ze, 'error')
plt.show()
In [ ]: