Installation¶

The Burger’s equation¶

The problem reads

The previous \(\DdQq{1}{2}\) scheme can simulate the Burger’s equation by modifying the equilibrium value of the moment of order \(1\) \(\mke{1}\). It now reads \(\mke{1} = {\mk{0}}^2/2\).

More generaly, the simulated equation is into the conservative form

the equilibrium has to be taken to \(\mke{1}=\varphi(\mk{0})\).

We just have to modify the equilibrium and the initialization of the previous example to simulate the Burger’s equation. The initial condition can be a discontinuous function in order to simulate Riemann problems. Note that the function f2m, m2f, relaxation, and transport are unchanged.

[7]:

def equilibrium(m0):
    return .5*m0**2

def initialize(mesh, la):
    ug, ud = 0.25, -0.15
    xmin, xmax = .5*np.sum(mesh[:2]), .5*np.sum(mesh[-2:])
    xc = xmin + .5*(xmax-xmin)
    m0 = ug*(mesh<xc) + ud*(mesh>xc) + .5*(ug+ud)*(mesh==xc)
    m1 = equilibrium(m0)
    f0 = np.empty(m0.shape)
    f1 = np.empty(m0.shape)
    return f0, f1, m0, m1

def relaxation(m0, m1, s):
    m1[:] = (1-s)*m1 + s*equilibrium(m0)

# parameters
Tf = 1. # final time
N = 128 # number of points in space
la = 1. # scheme velocity
s = 1.8 # relaxation parameter
# initialization
x = mesh(N)     # mesh
dx = x[1]-x[0]  # space step
dt = dx/la      # time step
f0, f1, m0, m1 = initialize(x, la)
plt.figure(1)
plt.plot(x[1:-1], m0[1:-1], 'b', label='initial')
# time loops
t = 0.
while (t<Tf):
    t += dt
    relaxation(m0, m1, s)
    m2f(f0, f1, m0, m1, la)
    transport(f0, f1)
    f2m(f0, f1, m0, m1, la)
plt.plot(x[1:-1], m0[1:-1], 'r', label='final')
plt.title('Burgers equation')
plt.legend(loc='best')
plt.show()

We can test different values of the relaxation parameter \(s\). In particular, we observe that the scheme remains stable if \(s\in[0,2]\). More \(s\) is small, more the numerical diffusion is important and if \(s\) is close to \(2\), oscillations appear behind the shock.

In order to simulate a Riemann problem, the boundary conditions have to be modified. A classical way is to impose entry conditions for hyperbolic problems. The lattice Boltzmann methods lend themselves very well to that conditions: the scheme only needs the distributions corresponding to a velocity that goes inside the domain. Nevertheless, on a physical edge where the flux is going outside, a non physical distribution that goes inside has to be imposed. A first simple way is to leave the initial value: this is correct while the discontinuity does not reach the edge. A second way is to impose Neumann condition by repeating the inner value.

We modify the previous script to take into account these new boundary conditions.

[8]:

def transport(f0, f1):
    # Neumann boundary conditions
    f0[-1] = f0[-2]
    f1[0] = f1[1]
    # transport
    f0[1:-1] = f0[2:]
    f1[1:-1] = f1[:-2]

# parameters
Tf = 1. # final time
N = 128 # number of points in space
la = 1. # scheme velocity
s = 1.8 # relaxation parameter

# initialization
x = mesh(N)     # mesh
dx = x[1]-x[0]  # space step
dt = dx/la      # time step
f0, f1, m0, m1 = initialize(x, la)
plt.figure(1)
plt.plot(x[1:-1], m0[1:-1], 'b', label='initial')
# time loops
t = 0.
while (t<Tf):
    t += dt
    relaxation(m0, m1, s)
    m2f(f0, f1, m0, m1, la)
    transport(f0, f1)
    f2m(f0, f1, m0, m1, la)
plt.plot(x[1:-1], m0[1:-1], 'r', label='final')
plt.title('Burgers equation')
plt.legend(loc='best')
plt.show()

The scheme \(\DdQq{1}{3}\)¶

The numerical simulation of this equation by a lattice Boltzmann scheme consists in the approximation of the solution on discret points of \((0,2\pi)\) at discret instants.

The spatial mesh is defined by using a numpy array. To simplify, the mesh is supposed to be uniform.

First, we import the package numpy and we create the spatial mesh. One phantom cell has to be added at each bound for the treatment of the boundary conditions.

[1]:

%matplotlib inline

[2]:

import numpy as np
import pylab as plt

def mesh(N):
    xmin, xmax = 0., 2.*np.pi
    dx = (xmax-xmin)/N
    x = np.linspace(xmin-.5*dx, xmax+.5*dx, N+2)
    return x

x = mesh(10)
plt.plot(x, 0.*x, 'sk')
plt.show()

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

  \[\mk{0}=\sum_{i=0}^{2} \fk{i}, \quad \mk{1}= \lambda \sum_{i=0}^{2} v_i \fk{i}, \quad \mk{2}= \frac{\lambda^2}{2} \sum_{i=0}^{2} v_i^2 \fk{i},\]

and their equilibrium values \(\mke{0} = \mk{0}\), \(\mke{1} = \mk{1}\), and \(\mke{2} = c^2/2 \mk{0}\).

• a relaxation parameter \(s\) lying in \([0,2]\).

In order to prepare the formalism of the package pylbm, we introduce the three polynomials that define the moments: \(P_0 = 1\), \(P_1=\lambda X\), and \(P_2=\lambda^2/2 X^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:

  \[\mks{2}(t,x)=(1-s)\mk{2}(t,x)+s\mke{2}(t,x).\]

• m2f:

  \[\begin{split}\begin{aligned}\fks{0}(t,x)&\;=\mk{0}(t,x)-2\mks{2}(t,x)/\lambda^2, \\ \fks{1}(t,x)&\;=\mk{1}(t,x)/(2\lambda)+\mks{2}(t,x)/\lambda^2, \\ \fks{2}(t,x)&\;=-\mk{1}(t,x)/(2\lambda)+\mks{2}(t,x)/\lambda^2.\end{aligned}\end{split}\]

• transport:

  \[\begin{split}\begin{aligned} \fk{0}(t+\dt, x)&\;=\fks{0}(t,x), \\ \fk{1}(t+\dt, x)&\;=\fks{1}(t,x-\dx), \\ \fk{2}(t+\dt, x)&\;=\fks{2}(t,x+\dx). \end{aligned}\end{split}\]

• f2m:

  \[\begin{split}\begin{aligned}\mk{0}(t+\dt,x)&\;=\fk{0}(t+\dt,x)+\fk{1}(t+\dt,x)+\fk{2}(t+\dt,x), \\ \mk{1}(t+\dt,x)&\;=\lambda\fk{1}(t+\dt,x)-\lambda\fk{2}(t+\dt,x), \\ \mk{2}(t+\dt,x)&\;=\tfrac{1}{2}\lambda^2\fk{1}(t+\dt,x)+\tfrac{1}{2}\lambda^2\fk{2}(t+\dt,x).\end{aligned}\end{split}\]

The moments of order \(0\), \(\mk{0}\), and of order \(1\), \(\mk{1}\), being conserved during the relaxation phase, the equivalent equations of this scheme read at first order

We implement a function equilibrium that computes the equilibrium value \(\mke{2}\), the moment of order \(0\), \(\mk{0}\), and the velocity \(c\) being given in argument.

[3]:

def equilibrium(m0, c):
    return .5*c**2*m0

We create three vectors \(\mk{0}\), \(\mk{1}\), and \(\mk{2}\) with shape the shape of the mesh and initialize them. The moments of order \(0\) and \(1\) should contain the initial value of the unknowns \(\rho\) and \(q\), and the moment of order \(2\) the corresponding equilibrium value.

We create also three vectors \(\fk{0}\), \(\fk{1}\) and \(\fk{2}\).

[4]:

def initialize(mesh, c, la):
    m0 = np.sin(mesh)
    m1 = np.zeros(mesh.shape)
    m2 = equilibrium(m0, c)
    f0 = np.empty(m0.shape)
    f1 = np.empty(m0.shape)
    f2 = np.empty(m0.shape)
    return f0, f1, f2, m0, m1, m2

## Periodic boundary conditions

We implement the four elementary functions f2m, relaxation, m2f, and transport. In the transport function, the boundary conditions should be implemented: we will use periodic conditions by copying the informations in the phantom cells.

[5]:

def f2m(f0, f1, f2, m0, m1, m2, la):
    m0[:] = f0 + f1 + f2
    m1[:] = la * (f2 - f1)
    m2[:] = .5* la**2 * (f1 + f2)

def m2f(f0, f1, f2, m0, m1, m2, la):
    f0[:] = m0 - 2./la**2 * m2
    f1[:] = -.5/la * m1 + 1/la**2 * m2
    f2[:] = .5/la * m1 + 1/la**2 * m2

def relaxation(m0, m1, m2, c, s):
    m2[:] *= (1-s)
    m2[:] += s*equilibrium(m0, c)

def transport(f0, f1, f2):
    # periodic boundary conditions
    f1[-1] = f1[1]
    f2[0] = f2[-2]
    # transport
    f1[1:-1] = f1[2:]
    f2[1:-1] = f2[:-2]

We compute and we plot the numerical solution at time \(T_f=2\pi\).

[6]:

# parameters
c = 1   # velocity for the transport equation
Tf = 2.*np.pi # final time
N = 128 # number of points in space
la = 1. # scheme velocity
s = 2.  # relaxation parameter
# initialization
x = mesh(N)     # mesh
dx = x[1]-x[0]  # space step
dt = dx/la      # time step
f0, f1, f2, m0, m1, m2 = initialize(x, c, la)
plt.figure(1)
plt.plot(x[1:-1], m0[1:-1], 'r', label=r'$\rho$')
plt.plot(x[1:-1], m1[1:-1], 'b', label=r'$q$')
plt.title('Initial moments')
plt.legend(loc='best')
# time loops
nt = int(Tf/dt)
m2f(f0, f1, f2, m0, m1, m2, la)
for k in range(nt):
    transport(f0, f1, f2)
    f2m(f0, f1, f2, m0, m1, m2, la)
    relaxation(m0, m1, m2, c, s)
    m2f(f0, f1, f2, m0, m1, m2, la)
plt.figure(2)
plt.plot(x[1:-1], m0[1:-1], 'r', label=r'$\rho$')
plt.plot(x[1:-1], m1[1:-1], 'b', label=r'$q$')
plt.title('Final moments')
plt.legend(loc='best')
plt.show()

Anti bounce back conditions¶

In order to take into account homogenous Dirichlet conditions over \(\rho\), we introduce the bounce back conditions. At edge \(x=0\), two points are involved: \(x_0=-\dx/2\) and \(x_1=\dx/2\). We impose \(\fk{1}(x_0)=-\fk{2}(x_1)\). And at edge \(x=2\pi\), the two involved points are \(x_{N}\) and \(x_{N+1}\). We impose \(\fk{2}(x_{N+1})=-\fk{1}(x_{N})\).

We modify the transport function to impose anti bounce back conditions. We can compare the solutions obtained with the two different boundary conditions.

[7]:

def transport(f0, f1, f2):
    # anti bounce back boundary conditions
    f1[-1] = -f2[-2]
    f2[0] = -f1[1]
    # transport
    f1[1:-1] = f1[2:]
    f2[1:-1] = f2[:-2]


# parameters
c = 1   # velocity for the transport equation
Tf = 2*np.pi # final time
N = 128 # number of points in space
la = 1. # scheme velocity
s = 2.  # relaxation parameter
# initialization
x = mesh(N)     # mesh
dx = x[1]-x[0]  # space step
dt = dx/la      # time step
f0, f1, f2, m0, m1, m2 = initialize(x, c, la)
plt.figure(1)
plt.plot(x[1:-1], m0[1:-1], 'r', label=r'$\rho$')
plt.plot(x[1:-1], m1[1:-1], 'b', label=r'$q$')
plt.title('Initial moments')
plt.legend(loc='best')
# time loops
nt = int(Tf/dt)
m2f(f0, f1, f2, m0, m1, m2, la)
for k in range(nt):
    transport(f0, f1, f2)
    f2m(f0, f1, f2, m0, m1, m2, la)
    relaxation(m0, m1, m2, c, s)
    m2f(f0, f1, f2, m0, m1, m2, la)
plt.figure(2)
plt.plot(x[1:-1], m0[1:-1], 'r', label=r'$\rho$')
plt.plot(x[1:-1], m1[1:-1], 'b', label=r'$q$')
plt.title('Final moments')
plt.legend(loc='best')
plt.show()

[ ]:

The scheme \(\DdQq{1}{3}\)¶

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

  \[\mk{0}=\sum_{i=0}^{2} \fk{i}, \quad \mk{1}= \sum_{i=0}^{2} v_i \fk{i}, \quad \mk{2}= \frac{1}{2} \sum_{i=0}^{2} v_i^2 \fk{i},\]

  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:

  \[\begin{split}\begin{aligned}\mks{1}(t,x)&=(1-s_1)\mk{1}(t,x)+s_1\mke{1}(t,x),\\ \mks{2}(t,x)&=(1-s_2)\mk{2}(t,x)+s_2\mke{2}(t,x).\end{aligned}\end{split}\]

• m2f:

  \[\begin{split}\begin{aligned}\fks{0}(t,x)&\;=\mk{0}(t,x)-2\mks{2}(t,x), \\ \fks{1}(t,x)&\;=\mks{1}(t,x)/2+\mks{2}(t,x), \\ \fks{2}(t,x)&\;=-\mks{1}(t,x)/2+\mks{2}(t,x).\end{aligned}\end{split}\]

• transport:

  \[\begin{split}\begin{aligned} \fk{0}(t+\dt, x)&\;=\fks{0}(t,x), \\ \fk{1}(t+\dt, x)&\;=\fks{1}(t,x-\dx), \\ \fk{2}(t+\dt, x)&\;=\fks{2}(t,x+\dx). \end{aligned}\end{split}\]

• f2m:

  \[\begin{split}\begin{aligned}\mk{0}(t+\dt,x)&\;=\fk{0}(t+\dt,x)+\fk{1}(t+\dt,x)+\fk{2}(t+\dt,x), \\ \mk{1}(t+\dt,x)&\;=\fk{1}(t+\dt,x)-\fk{2}(t+\dt,x), \\ \mk{2}(t+\dt,x)&\;=\tfrac{1}{2}\fk{1}(t+\dt,x)+\tfrac{1}{2}\fk{2}(t+\dt,x).\end{aligned}\end{split}\]

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:

Using pylbm¶

pylbm uses Python dictionary to describe the simulation. In the following, we will build this dictionary step by step.

The geometry¶

In pylbm, the geometry is defined by a box and a label for the boundaries.

[2]:

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 information |
+----------------------+
    - spatial dimension: 1
    - bounds of the box: [0. 1.]
    - labels: [0, 0]

The stencil¶

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]\).

[3]:

dico_sten = {
    'dim': 1,
    'schemes':[{'velocities': list(range(3))}],
}
sten = pylbm.Stencil(dico_sten)
print(sten)
sten.visualize();

+---------------------+
| Stencil information |
+---------------------+
    - spatial dimension: 1
    - minimal velocity in each direction: [-1]
    - maximal velocity in each direction: [1]
    - information for each elementary stencil:
        stencil 0
            - number of velocities: 3
            - velocities
                (0: 0)
                (1: 1)
                (2: -1)

The domain¶

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.

[4]:

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 information |
+--------------------+
    - spatial dimension: 1
    - space step: 0.1
    - with halo:
        bounds of the box: [-0.05] x [1.05]
        number of points: [12]
    - without halo:
        bounds of the box: [0.05] x [0.95]
        number of points: [10]

    +----------------------+
    | Geometry information |
    +----------------------+
        - spatial dimension: 1
        - bounds of the box: [0. 1.]
        - labels: [0, 0]

[5]:

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],
        }
    ],
}

sch = pylbm.Scheme(dico_sch)
print(sch)

+--------------------+
| Scheme information |
+--------------------+
    - spatial dimension: 1
    - number of schemes: 1
    - number of velocities: 3
    - conserved moments: [u]

    +----------+
    | Scheme 0 |
    +----------+
        - velocities
            (0: 0)
            (1: 1)
            (2: -1)

        - polynomials

            ⎡1 ⎤
            ⎢  ⎥
            ⎢X ⎥
            ⎢  ⎥
            ⎢ 2⎥
            ⎢X ⎥
            ⎢──⎥
            ⎣2 ⎦

        - equilibria

            ⎡  u  ⎤
            ⎢     ⎥
            ⎢ 0.0 ⎥
            ⎢     ⎥
            ⎣0.5⋅u⎦

        - relaxation parameters

            ⎡       0.0       ⎤
            ⎢                 ⎥
            ⎢0.666666666666667⎥
            ⎢                 ⎥
            ⎣       1.0       ⎦

    - moments matrices

        ⎡1  1    1 ⎤
        ⎢          ⎥
        ⎢0  10  -10⎥
        ⎢          ⎥
        ⎣0  50  50 ⎦

    - inverse of moments matrices

        ⎡1    0    -1/50⎤
        ⎢               ⎥
        ⎢0  1/20   1/100⎥
        ⎢               ⎥
        ⎣0  -1/20  1/100⎦

[6]:

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.AntiBounceBack,}},
    },
    'generator': 'numpy'
}

sol = pylbm.Simulation(dico)
print(sol)

+------------------------+
| Simulation information |
+------------------------+

    +--------------------+
    | Domain information |
    +--------------------+
        - spatial dimension: 1
        - space step: 0.1
        - with halo:
            bounds of the box: [-0.05] x [1.05]
            number of points: [12]
        - without halo:
            bounds of the box: [0.05] x [0.95]
            number of points: [10]

        +----------------------+
        | Geometry information |
        +----------------------+
            - spatial dimension: 1
            - bounds of the box: [0. 1.]
            - labels: [0, 0]

    +--------------------+
    | Scheme information |
    +--------------------+
        - spatial dimension: 1
        - number of schemes: 1
        - number of velocities: 3
        - conserved moments: [u]

        +----------+
        | Scheme 0 |
        +----------+
            - velocities
                (0: 0)
                (1: 1)
                (2: -1)

            - polynomials

                ⎡1 ⎤
                ⎢  ⎥
                ⎢X ⎥
                ⎢  ⎥
                ⎢ 2⎥
                ⎢X ⎥
                ⎢──⎥
                ⎣2 ⎦

            - equilibria

                ⎡  u  ⎤
                ⎢     ⎥
                ⎢ 0.0 ⎥
                ⎢     ⎥
                ⎣0.5⋅u⎦

            - relaxation parameters

                ⎡       0.0       ⎤
                ⎢                 ⎥
                ⎢0.666666666666667⎥
                ⎢                 ⎥
                ⎣       1.0       ⎦

        - moments matrices

            ⎡1  1    1 ⎤
            ⎢          ⎥
            ⎢0  10  -10⎥
            ⎢          ⎥
            ⎣0  50  50 ⎦

        - inverse of moments matrices

            ⎡1    0    -1/50⎤
            ⎢               ⎥
            ⎢0  1/20   1/100⎥
            ⎢               ⎥
            ⎣0  -1/20  1/100⎦

Run a simulation¶

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.

[7]:

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.AntiBounceBack,}},
    },
    '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();

The scheme \(\DdQq{2}{5}\)¶

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

  \[\mk{0}=\sum_{i=0}^{4} \fk{i}, \quad \mk{1}= \sum_{i=0}^{4} v_{ix} \fk{i}, \quad \mk{2}= \sum_{i=0}^{4} v_{iy} \fk{i}, \quad \mk{3}= \frac{1}{2} \sum_{i=0}^{5} (v_{ix}^2+v_{iy}^2) \fk{i}, \quad \mk{4}= \frac{1}{2} \sum_{i=0}^{5} (v_{ix}^2-v_{iy}^2) \fk{i},\]

  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:

  \[\begin{split}\begin{aligned}\mks{1}(t,x,y)&=(1-s_1)\mk{1}(t,x,y)+s_1\mke{1}(t,x,y),\\ \mks{2}(t,x,y)&=(1-s_1)\mk{2}(t,x,y)+s_1\mke{2}(t,x,y),\\ \mks{3}(t,x,y)&=(1-s_2)\mk{3}(t,x,y)+s_2\mke{3}(t,x,y),\\ \mks{4}(t,x,y)&=(1-s_2)\mk{4}(t,x,y)+s_2\mke{4}(t,x,y).\end{aligned}\end{split}\]

• m2f:

  \[\begin{split}\begin{aligned}\fks{0}(t,x,y)&\;=\mk{0}(t,x,y)-2\mks{3}(t,x,y), \\ \fks{1}(t,x,y)&\;=\tfrac{1}{2}(\phantom{-}\mks{1}(t,x,y)+\mks{3}(t,x,y)+\mks{4}(t,x,y)), \\ \fks{2}(t,x,y)&\;=\tfrac{1}{2}(\phantom{-}\mks{2}(t,x,y)+\mks{3}(t,x,y)-\mks{4}(t,x,y)), \\ \fks{3}(t,x,y)&\;=\tfrac{1}{2}(-\mks{1}(t,x,y)+\mks{3}(t,x,y)+\mks{4}(t,x,y)), \\ \fks{4}(t,x,y)&\;=\tfrac{1}{2}(-\mks{2}(t,x,y)+\mks{3}(t,x,y)-\mks{4}(t,x,y)).\end{aligned}\end{split}\]

• transport:

  \[\begin{split}\begin{aligned} \fk{0}(t+\dt, x,y)&\;=\fks{0}(t,x,y), \\ \fk{1}(t+\dt, x,y)&\;=\fks{1}(t,x-\dx,y), \\ \fk{2}(t+\dt, x,y)&\;=\fks{2}(t,x,y-\dx), \\ \fk{3}(t+\dt, x,y)&\;=\fks{3}(t,x+\dx,y), \\ \fk{4}(t+\dt, x,y)&\;=\fks{4}(t,x,y+\dx). \end{aligned}\end{split}\]

• f2m:

  \[\begin{split}\begin{aligned}\mk{0}(t+\dt,x,y)&\;=\fk{0}(t+\dt,x,y)+\fk{1}(t+\dt,x,y)+\fk{2}(t+\dt,x,y)\\&\;\phantom{=}+\fk{3}(t+\dt,x,y)+\fk{4}(t+\dt,x,y), \\ \mk{1}(t+\dt,x,y)&\;=\fk{1}(t+\dt,x,y)-\fk{3}(t+\dt,x,y), \\ \mk{2}(t+\dt,x,y)&\;=\fk{2}(t+\dt,x,y)-\fk{4}(t+\dt,x,y), \\ \mk{3}(t+\dt,x,y)&\;=\tfrac{1}{2}(\fk{1}(t+\dt,x,y)+\fk{2}(t+\dt,x,y)+\fk{3}(t+\dt,x,y)+\fk{4}(t+\dt,x,y)), \\ \mk{4}(t+\dt,x,y)&\;=\tfrac{1}{2}(\fk{1}(t+\dt,x,y)-\fk{2}(t+\dt,x,y)+\fk{3}(t+\dt,x,y)-\fk{4}(t+\dt,x,y)).\end{aligned}\end{split}\]

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:

Using pylbm¶

pylbm uses Python dictionary to describe the simulation. In the following, we will build this dictionary step by step.

The geometry¶

In pylbm, the geometry is defined by a box and a label for the boundaries. We define here a square \((0, 1)^2\).

[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 information |
+----------------------+
    - spatial dimension: 2
    - bounds of the box: [0. 1.] x [0. 1.]
    - labels: [0, 0, 0, 0]

[3]:

dico_sten = {
    'dim':2,
    'schemes': [
        {'velocities': list(range(5))}
    ],
}
sten = pylbm.Stencil(dico_sten)
print(sten)
sten.visualize();