\[\renewcommand{\DdQq}[2]{{\mathrm D}_{#1}{\mathrm Q}_{#2}} \renewcommand{\drondt}{\partial_t} \renewcommand{\drondx}{\partial_x} \renewcommand{\drondtt}{\partial_{tt}} \renewcommand{\drondxx}{\partial_{xx}} \renewcommand{\drondyy}{\partial_{yy}} \renewcommand{\dx}{\Delta x} \renewcommand{\dt}{\Delta t} \renewcommand{\grandO}{{\mathcal O}} \renewcommand{\density}[2]{\,f_{#1}^{#2}} \renewcommand{\fk}[1]{\density{#1}{\vphantom{\star}}} \renewcommand{\fks}[1]{\density{#1}{\star}} \renewcommand{\moment}[2]{\,m_{#1}^{#2}} \renewcommand{\mk}[1]{\moment{#1}{\vphantom{\star}}} \renewcommand{\mke}[1]{\moment{#1}{e}} \renewcommand{\mks}[1]{\moment{#1}{\star}}\]

Analyze your scheme

Two of the biggest problems encountered when starting to use lattice Bolzmann methods are

pylbm tries to give you some ideas to solve them.

For the first one, pylbm can give you the first and second order coefficients of your physical equation (in a next release, it will be possible to have also the third and the fourth order terms). To have better results, it is important to keep SymPy symbols in your scheme as long as possible. Thus, we can see the influence of theses parameters on the physical equations.

For the second one, once you know that you have the good physical equation, pylbm allows to check the linear stability region for these parameters arround a given linearized state. We provide widgets inside a notebook or for a Python script to modify easily these parameters and show the result on the figure interactively.

We believe that these two tools will make it easier for as many people as possible to become familiar with the lattice Boltzmann methods and, in the end, allow them to implement their own schemes.

Let’s take a simple example to illustrate how it works.

Assume that you want to solve the advection equation for 1D problem

\[\begin{split}\begin{aligned} &\drondt u = c \drondx u, && t>0, \quad x\in(0,1),\\ &u(t=0,x) = u_0(x), && x\in(0,1) \\ &u(t,x=0) = u(t,x=1), && t>0. \end{aligned}\end{split}\]

And you already have a lattice Boltzmann scheme that you want to try: the \(\DdQq{1}{3}\) scheme given by

We can write this scheme into pylbm as

[1]:
import sympy as sp

# Symbolic variables definitions
U, X = sp.symbols('u, X')
C, LA, S0, S1 = sp.symbols('c, lambda, s_0, s_1', constants=True)

# The D1Q3 LBM scheme
adv_scheme = {
    'dim': 1,
    'scheme_velocity': LA,
    'schemes': [
        {
            'velocities': list(range(3)),
            'conserved_moments': U,
            'polynomials': [1, X, X**2/2],
            'relaxation_parameters': [0., S0, S1],
            'equilibrium': [U, C*U, C**2*U/2],
        },
    ],
    'parameters': {LA: 1,
                   S0: 1.9,
                   S1: 1.9,
                   C: 1},
}

Let’s create the scheme and look at the information given by pylbm

[2]:
import pylbm

scheme = pylbm.Scheme(adv_scheme)
[3]:
print(scheme)
+--------------------+
| 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  ⎤
            ⎢    ⎥
            ⎢c⋅u ⎥
            ⎢    ⎥
            ⎢ 2  ⎥
            ⎢c ⋅u⎥
            ⎢────⎥
            ⎣ 2  ⎦

        - relaxation parameters

            ⎡0.0⎤
            ⎢   ⎥
            ⎢s₀ ⎥
            ⎢   ⎥
            ⎣s₁ ⎦

    - moments matrices

        ⎡1  1   1 ⎤
        ⎢         ⎥
        ⎢0  λ   -λ⎥
        ⎢         ⎥
        ⎢    2   2⎥
        ⎢   λ   λ ⎥
        ⎢0  ──  ──⎥
        ⎣   2   2 ⎦

    - inverse of moments matrices

        ⎡        -2 ⎤
        ⎢1   0   ───⎥
        ⎢          2⎥
        ⎢         λ ⎥
        ⎢           ⎥
        ⎢    1   1  ⎥
        ⎢0  ───  ── ⎥
        ⎢   2⋅λ   2 ⎥
        ⎢        λ  ⎥
        ⎢           ⎥
        ⎢   -1   1  ⎥
        ⎢0  ───  ── ⎥
        ⎢   2⋅λ   2 ⎥
        ⎣        λ  ⎦


We can see here that we have described one scheme with three 1D velocities. The moment matrix gives how to find the moments from the distribution functions.

Let’s check now if we solve the good physical equations.

[4]:
pde = pylbm.EquivalentEquation(scheme)
[5]:
print(pde)
+----------------------+
| Equivalent Equations |
+----------------------+
    The equation is


    d        d       ∂ ⎛    d    ⎞
    ──(Fx) + ──(U) = ──⎜Bxx⋅──(U)⎟
    dx       dt      ∂x⎝    dx   ⎠

    where


    U = [u]

    Fx = [c⋅u]

    Bxx = [0]

pylbm gives the first and second order terms. In the next release, you will also have access to the third and fourth terms. Our scheme solves the advection equation as exepected.

We can now study the stability of this scheme. Many notions of stability exist and can be used. In pylbm, we focus on a linear notion by computing the eigenvalues of the linear operator corresponding to one time step. The scheme will be considered as stable if all these eigenvalues stay inside the unit circle (as complex values). This notion is sufficient for linear scheme but just gives partial informations for non-linear scheme.

[6]:
stab = pylbm.Stability(scheme)

# linearization around a state
uo = 1.

stab.visualize(
    {
        'linearization': {
            U: uo,
        },
        'parameters': {
            LA: {
                'range': [1, 20],
                'init': 1,
                'step': .1
            },
            U: {
                'range': [0, 20],
                'init': uo,
                'step': .1
            },
            C: {
                'range': [0, 20],
                'init': 1,
                'step': .1
            },
            S0: {
                'range': [0, 2],
                'init': 1.9,
                'step': .1
            },
            S1: {
                'range': [0, 2],
                'init': 1.9,
                'step': .1
            },
        },
        'number_of_wave_vectors': 1024,
    }
)