The Scheme

With pylbm, elementary schemes can be gathered and coupled through the equilibrium in order to simplify the implementation of the vectorial schemes. Of course, the user can implement a single elementary scheme and then recover the classical framework of the d’Humières schemes.

For pylbm, the scheme is performed through a dictionary. The generalized d’Humières framework for vectorial schemes is used [dH92], [G14]. In the first section, we describe how build an elementary scheme. Then the vectorial schemes are introduced as coupled elementary schemes.

The elementary schemes

Let us first consider a regular lattice \(L\) in dimension \(d\) with a typical mesh size \(dx\), and the time step \(dt\). The scheme velocity \(\lambda\) is then defined by \(\lambda = dx/dt\). We introduce a set of \(q\) velocities adapted to this lattice \(\{v_0, \ldots, v_{q-1}\}\), that is to say that, if \(x\) is a point of the lattice \(L\), the point \(x+v_jdt\) is on the lattice for every \(j\in\{0,\ldots,q{-}1\}\).

The aim of the \(DdQq\) scheme is to compute a distribution function vector \({\boldsymbol f} = (f_0,\ldots,f_{q-1})\) on the lattice \(L\) at discret values of time. The scheme splits into two phases: the relaxation and the transport. That is, the passage from the time \(t\) to the time \(t+dt\) consists in the succession of these two phases.

  • the relaxation phase

    This phase, also called collision, is local in space: on every site \(x\) of the lattice, the values of the vector \({\boldsymbol f}\) are modified, the result after the collision being denoted by \({\boldsymbol f}^\star\). The operator of collision is a linear operator of relaxation toward an equilibrium value denoted \({\boldsymbol f}^{\textrm{eq}}\).

    pylbm uses the framework of d’Humières: the linear operator of the collision is diagonal in a special basis called moments denoted by \({\boldsymbol m} = (m_0,\ldots,m_{q-1})\). The change-of-basis matrix \(M\) is such that \({\boldsymbol m} = M{\boldsymbol f}\) and \({\boldsymbol f} = M^{-1}{\boldsymbol m}\). In the basis of the moments, the collision operator then just reads

    \begin{equation*} m_k^\star = m_k - s_k (m_k - m_k^{\textrm{eq}}), \qquad 0\leqslant k\leqslant q{-}1, \end{equation*}

    where \(s_k\) is the relaxation parameter associated to the kth moment. The kth moment is said conserved during the collision if the associated relaxation parameter \(s_k=0\).

    By analogy with the kinetic theory, the change-of-basis matrix \(M\) is defined by a set of polynomials with \(d\) variables \((P_0,\ldots,P_{q-1})\) by

    \begin{equation*} M_{ij} = P_i(v_j). \end{equation*}
  • the transport phase

    This phase just consists in a shift of the indices and reads

    \begin{equation*} f_j(x, t+dt) = f_j^\star(x-v_jdt, t), \qquad 0\leqslant j\leqslant q{-}1, \end{equation*}

Notations

The scheme is defined and build through a dictionary in pylbm. Let us first list the several key words of this dictionary:

  • dim: the spatial dimension. This argument is optional if the geometry is known, that is if the dimension can be computed through the list of the variables;

  • scheme_velocity: the velocity of the scheme denoted by \(\lambda\) in the previous section and defined as the spatial step over the time step (\(\lambda=dx/dt\)) ;

  • schemes: the list of the schemes. In pylbm, several coupled schemes can be used, the coupling being done through the equilibrium values of the moments. Some examples with only one scheme and with more than one schemes are given in the next sections. Each element of the list should be a dictionay with the following key words:

    • velocities: the list of the velocity indices,

    • conserved_moments: the list of the conserved moments (list of symbolic variables),

    • polynomials: the list of the polynomials that define the moments, the polynomials are built with the symbolic variables X, Y, and Z,

    • equilibrium: the list of the equilibrium value of the moments,

    • relaxation_parameters: the list of the relaxation parameters, (by convention, the relaxation parameter of a conserved moment is taken to 0).

Examples in 1D

script

\(D_1Q_2\) for the advection

A velocity \(c\in{\mathbb R}\) being given, the advection equation reads

\begin{equation*} \partial_t u(t,x) + c \partial_x u(t,x) = 0, \qquad t>0, x\in{\mathbb R}. \end{equation*}

Taken for instance \(c=0.5\), the following scheme can be used:

import pylbm

u, X = sp.symbols("u, X")

d = {
    "dim": 1,
    "scheme_velocity": 1.0,
    "schemes": [
        {
            "velocities": [1, 2],
            "conserved_moments": u,
            "polynomials": [1, X],
            "equilibrium": [u, 0.5 * u],
            "relaxation_parameters": [0.0, 1.9],
        },
    ],
}
s = pylbm.Scheme(d)
print(s)

The dictionary d is used to set the dimension to 1, the scheme velocity to 1. The used scheme has two velocities: the first one \(v_0=1\) and the second one \(v_1=-1\). The polynomials that define the moments are \(P_0 = 1\) and \(P_1 = X\) so that the matrix of the moments is

\begin{equation*} M = \begin{pmatrix} 1&1\\ 1&-1 \end{pmatrix} \end{equation*}

with the convention \(M_{ij} = P_i(v_j)\). Then, there are two distribution functions \(f_0\) and \(f_1\) that move at the velocities \(v_0\) and \(v_1\), and two moments \(m_0 = f_0+f_1\) and \(m_1 = f_0-f_1\). The first moment \(m_0\) is conserved during the relaxation phase (as the associated relaxation parameter is set to 0), while the second moment \(m_1\) relaxes to its equilibrium value \(0.5m_0\) with a relaxation parameter \(1.9\) by the relation

\begin{equation*} m_1^\star = m_1 - 1.9 (m_1 - 0.5m_0). \end{equation*}

script

\(D_1Q_2\) for Burger’s

The Burger’s equation reads

\begin{equation*} \partial_t u(t,x) + \tfrac{1}{2}\partial_x u^2(t,x) = 0, \qquad t>0, x\in{\mathbb R}. \end{equation*}

The following scheme can be used:

import pylbm

u, X = sp.symbols("u, X")

d = {
    "dim": 1,
    "scheme_velocity": 1.0,
    "schemes": [
        {
            "velocities": [1, 2],
            "conserved_moments": u,
            "polynomials": [1, X],
            "equilibrium": [u, 0.5 * u**2],
            "relaxation_parameters": [0.0, 1.9],
        },
    ],
}
s = pylbm.Scheme(d)
print(s)

The same dictionary has been used for this application with only one modification: the equilibrium value of the second moment \(m_1^{\textrm{eq}}\) is taken to \(\tfrac{1}{2}m_0^2\).

script

\(D_1Q_3\) for the wave equation

The wave equation is rewritten into the system of two partial differential equations

\begin{equation*} \left\{ \begin{aligned} &\partial_t u(t, x) + \partial_x v(t, x) = 0, & t>0, x\in{\mathbb R},\\ &\partial_t v(t, x) + c^2\partial_x u(t, x) = 0, & t>0, x\in{\mathbb R}. \end{aligned} \right. \end{equation*}

The following scheme can be used:

import pylbm

u, v, X = sp.symbols("u, v, X")

c = 0.5
d = {
    "dim": 1,
    "scheme_velocity": 1.0,
    "schemes": [
        {
            "velocities": [0, 1, 2],
            "conserved_moments": [u, v],
            "polynomials": [1, X, 0.5 * X**2],
            "equilibrium": [u, v, 0.5 * c**2 * u],
            "relaxation_parameters": [0.0, 0.0, 1.9],
        },
    ],
}
s = pylbm.Scheme(d)
print(s)

Examples in 2D

script

\(D_2Q_4\) for the advection

A velocity \((c_x, c_y)\in{\mathbb R}^2\) being given, the advection equation reads

\begin{equation*} \partial_t u(t, x, y) + c_x \partial_x u(t, x, y) + c_y \partial_y u(t, x, y) = 0, \qquad t>0, x, y \in{\mathbb R}. \end{equation*}

Taken for instance \(c_x=0.1, c_y=0.2\), the following scheme can be used:

import pylbm

u, X, Y = sp.symbols("u, X, Y")

d = {
    "dim": 2,
    "scheme_velocity": 1.0,
    "schemes": [
        {
            "velocities": [1, 2, 3, 4],
            "conserved_moments": u,
            "polynomials": [1, X, Y, X**2 - Y**2],
            "equilibrium": [u, 0.1 * u, 0.2 * u, 0.0],
            "relaxation_parameters": [0.0, 1.9, 1.9, 1.4],
        },
    ],
}
s = pylbm.Scheme(d)
print(s)

The dictionary d is used to set the dimension to 2, the scheme velocity to 1. The used scheme has four velocities: \(v_0=(1,0)\), \(v_1=(0,1)\), \(v_2=(-1,0)\), and \(v_3=(0,-1)\). The polynomials that define the moments are \(P_0 = 1\), \(P_1 = X\), \(P_2 = Y\), and \(P_3 = X^2-Y^2\) so that the matrix of the moments is

\begin{equation*} M = \begin{pmatrix} 1&1&1&1\\ 1&0&-1&0\\ 0&1&0&-1\\ 1&-1&1&-1 \end{pmatrix} \end{equation*}

with the convention \(M_{ij} = P_i(v_j)\). Then, there are four distribution functions \(f_j, 0\leq j\leq 3\) that move at the velocities \(v_j\), and four moments \(m_k = \sum_{j=0}^3 M_{kj}f_j\). The first moment \(m_0\) is conserved during the relaxation phase (as the associated relaxation parameter is set to 0), while the other moments \(m_k, 1\leq k\leq 3\) relaxe to their equilibrium values by the relations

\begin{equation*} \left\{ \begin{aligned} m_1^\star &= m_1 - 1.9 (m_1 - 0.1m_0),\\ m_2^\star &= m_2 - 1.9 (m_2 - 0.2m_0),\\ m_3^\star &= (1-1.4)m_3. \end{aligned} \right. \end{equation*}

script

\(D_2Q_9\) for Navier-Stokes

The system of the compressible Navier-Stokes equations reads

\begin{equation*} \left\{ \begin{aligned} &\partial_t \rho + \nabla{\cdot}(\rho {\boldsymbol u}) = 0,\\ &\partial_t (\rho {\boldsymbol u}) + \nabla{\cdot}(\rho {\boldsymbol u}{\otimes}{\boldsymbol u}) + \nabla p = \kappa \nabla (\nabla{\cdot}{\boldsymbol u}) + \eta \nabla{\cdot} \bigl(\nabla{\boldsymbol u} + (\nabla{\boldsymbol u})^T - \nabla{\cdot}{\boldsymbol u}{\mathbb I}\bigr), \end{aligned} \right. \end{equation*}

where we removed the dependency of all unknown on the variables \((t, x, y)\). The vector \({\boldsymbol x}\) stands for \((x, y)\) and, if \(\psi\) is a scalar function of \({\boldsymbol x}\) and \({\boldsymbol\phi}=(\phi_x,\phi_y)\) is a vectorial function of \({\boldsymbol x}\), the usual partial differential operators read

\begin{equation*} \begin{aligned} &\nabla\psi = (\partial_x\psi, \partial_y\psi),\\ &\nabla{\cdot}{\boldsymbol\phi} = \partial_x\phi_x + \partial_y\phi_y,\\ &\nabla{\cdot}({\boldsymbol\phi}{\otimes}{\boldsymbol\phi}) = (\nabla{\cdot}(\phi_x{\boldsymbol\phi}), \nabla{\cdot}(\phi_y{\boldsymbol\phi})). \end{aligned} \end{equation*}

The coefficients \(\kappa\) and \(\eta\) are the bulk and the shear viscosities.

The following dictionary can be used to simulate the system of Navier-Stokes equations in the limit of small velocities:

import sympy as sp
import pylbm

rho, qx, qy, X, Y = sp.symbols("rho, qx, qy, X, Y")

dx = 1.0 / 256  # space step
eta = 1.25e-5  # shear viscosity
kappa = 10 * eta  # bulk viscosity
sb = 1.0 / (0.5 + kappa * 3.0 / dx)
ss = 1.0 / (0.5 + eta * 3.0 / dx)
d = {
    "dim": 2,
    "scheme_velocity": 1.0,
    "schemes": [
        {
            "velocities": list(range(9)),
            "conserved_moments": [rho, qx, qy],
            "polynomials": [
                1,
                X,
                Y,
                3 * (X**2 + Y**2) - 4,
                (9 * (X**2 + Y**2) ** 2 - 21 * (X**2 + Y**2) + 8) / 2,
                3 * X * (X**2 + Y**2) - 5 * X,
                3 * Y * (X**2 + Y**2) - 5 * Y,
                X**2 - Y**2,
                X * Y,
            ],
            "relaxation_parameters": [0.0, 0.0, 0.0, sb, sb, sb, sb, ss, ss],
            "equilibrium": [
                rho,
                qx,
                qy,
                -2 * rho + 3 * qx**2 + 3 * qy**2,
                rho + 3 / 2 * qx**2 + 3 / 2 * qy**2,
                -qx,
                -qy,
                qx**2 - qy**2,
                qx * qy,
            ],
        },
    ],
}
s = pylbm.Scheme(d)
print(s)

The scheme generated by the dictionary is the 9 velocities scheme with orthogonal moments introduced in [QdHL92]

Examples in 3D

script

\(D_3Q_6\) for the advection

A velocity \((c_x, c_y, c_z)\in{\mathbb R}^2\) being given, the advection equation reads

\begin{equation*} \partial_t u(t, x, y, z) + c_x \partial_x u(t, x, y, z) + c_y \partial_y u(t, x, y, z) + c_z \partial_z u(t, x, y, z) = 0, \quad t>0, x, y, z \in{\mathbb R}. \end{equation*}

Taken for instance \(c_x=0.1, c_y=-0.1, c_z=0.2\), the following scheme can be used:

import sympy as sp
import pylbm

u, X, Y, Z = sp.symbols("u, X, Y, Z")

cx, cy, cz = 0.1, -0.1, 0.2
d = {
    "dim": 3,
    "scheme_velocity": 1.0,
    "schemes": [
        {
            "velocities": list(range(1, 7)),
            "conserved_moments": u,
            "polynomials": [1, X, Y, Z, X**2 - Y**2, X**2 - Z**2],
            "equilibrium": [u, cx * u, cy * u, cz * u, 0.0, 0.0],
            "relaxation_parameters": [0.0, 1.5, 1.5, 1.5, 1.5, 1.5],
        },
    ],
}
s = pylbm.Scheme(d)
print(s)

The dictionary d is used to set the dimension to 3, the scheme velocity to 1. The used scheme has six velocities: \(v_0=(0,0,1)\), \(v_1=(0,0,-1)\), \(v_2=(0,1,0)\), \(v_3=(0,-1,0)\), \(v_4=(1,0,0)\), and \(v_5=(-1,0,0)\). The polynomials that define the moments are \(P_0 = 1\), \(P_1 = X\), \(P_2 = Y\), \(P_3 = Z\), \(P_4 = X^2-Y^2\), and \(P_5 = X^2-Z^2\) so that the matrix of the moments is

\begin{equation*} M = \begin{pmatrix} 1&1&1&1&1&1\\ 0&0&0&0&1&-1\\ 0&0&1&-1&0&0\\ 1&-1&0&0&0&0\\ 0&0&-1&-1&1&1\\ -1&-1&0&0&1&1 \end{pmatrix} \end{equation*}

with the convention \(M_{ij} = P_i(v_j)\). Then, there are six distribution functions \(f_j, 0\leq j\leq 5\) that move at the velocities \(v_j\), and six moments \(m_k = \sum_{j=0}^5 M_{kj}f_j\). The first moment \(m_0\) is conserved during the relaxation phase (as the associated relaxation parameter is set to 0), while the other moments \(m_k, 1\leq k\leq 5\) relaxe to their equilibrium values by the relations

\begin{equation*} \left\{ \begin{aligned} m_1^\star &= m_1 - 1.5 (m_1 - 0.1m_0),\\ m_2^\star &= m_2 - 1.5 (m_2 + 0.1m_0),\\ m_3^\star &= m_3 - 1.5 (m_3 - 0.2m_0),\\ m_4^\star &= (1-1.5)m_4,\\ m_5^\star &= (1-1.5)m_5. \end{aligned} \right. \end{equation*}

The vectorial schemes

With pylbm, vectorial schemes can be built easily by using a list of elementary schemes. Each elementary scheme is given by a dictionary as in the previous section. The conserved moments of all the elementary schemes can be used in the equilibrium values of the non conserved moments, in order to couple the schemes. For more details on the vectorial schemes, the reader can refer to [G14].

Examples in 1D

script

\(D_1Q_{2,2}\) for the shallow water equation

A constant \(g\in{\mathbb R}\) being given, the shallow water system reads

\begin{align*} &\partial_t h(t,x) + \partial_x q(t,x) = 0, &\qquad t>0, x\in{\mathbb R},\\ &\partial_t q(t,x) + \partial_x \bigl(q^2(t,x)/h(t,x) + gh^2(t,x)/2\bigr) = 0, &\qquad t>0, x\in{\mathbb R}. \end{align*}

Taken for instance \(g=1\), the following scheme can be used:

import pylbm

# parameters
h, q, X, LA, g = sp.symbols("h, q, X, LA, g")
la = 2.0  # velocity of the scheme
s_h, s_q = 1.7, 1.5  # relaxation parameters

d = {
    "dim": 1,
    "scheme_velocity": la,
    "schemes": [
        {
            "velocities": [1, 2],
            "conserved_moments": h,
            "polynomials": [1, LA * X],
            "relaxation_parameters": [0, s_h],
            "equilibrium": [h, q],
        },
        {
            "velocities": [1, 2],
            "conserved_moments": q,
            "polynomials": [1, LA * X],
            "relaxation_parameters": [0, s_q],
            "equilibrium": [q, q**2 / h + 0.5 * g * h**2],
        },
    ],
    "parameters": {LA: la, g: 1.0},
}
s = pylbm.Scheme(d)
print(s)

Two elementary schemes have been built, these two schemes are identical except for the equilibrium values of the non conserved moment and of the relaxation parameter: The first one is used to simulate the equation on \(h\) and the second one to simulate the equation on \(q\). For each scheme, the equilibrium value of the non conserved moment is equal to the flux of the corresponding equation: the equilibrium value of the kth scheme can so depend on all the conserved moments (and not only on those of the kth scheme).

Examples in 2D

script

\(D_2Q_{4,4,4}\) for the shallow water equation

A constant \(g\in{\mathbb R}\) being given, the shallow water system reads

\begin{align*} &\partial_t h(t,x,y) + \partial_x q_x(t,x,y) + \partial_y q_y(t,x,y) = 0, &\qquad t>0, x,y\in{\mathbb R},\\ &\partial_t q_x(t,x,y) + \partial_x \bigl(q_x^2(t,x,y)/h(t,x,y) + gh^2(t,x,y)/2\bigr) \\ & \phantom{\partial_t q_x(t,x,y)} + \partial_y \bigl( q_x(t,x,y)q_y(t,x,y)/h(t,x,y)\bigr) = 0, &\qquad t>0, x,y\in{\mathbb R},\\ &\partial_t q_y(t,x,y) + \partial_x \bigl( q_x(t,x,y)q_y(t,x,y)/h(t,x,y)\bigr)\\ & \phantom{\partial_t q_y(t,x,y)} + \partial_y \bigl(q_y^2(t,x,y)/h(t,x,y) + gh^2(t,x,y)/2\bigr) = 0, &\qquad t>0, x,y\in{\mathbb R}. \end{align*}

Taken for instance \(g=1\), the following scheme can be used:

import pylbm

X, Y, LA, g = sp.symbols("X, Y, LA, g")
h, qx, qy = sp.symbols("h, qx, qy")

# parameters
la = 4  # velocity of the scheme
s_h = [0.0, 2.0, 2.0, 1.5]
s_q = [0.0, 1.5, 1.5, 1.2]

vitesse = [1, 2, 3, 4]
polynomes = [1, LA * X, LA * Y, X**2 - Y**2]

d = {
    "dim": 2,
    "scheme_velocity": la,
    "schemes": [
        {
            "velocities": vitesse,
            "conserved_moments": h,
            "polynomials": polynomes,
            "relaxation_parameters": s_h,
            "equilibrium": [h, qx, qy, 0.0],
        },
        {
            "velocities": vitesse,
            "conserved_moments": qx,
            "polynomials": polynomes,
            "relaxation_parameters": s_q,
            "equilibrium": [qx, qx**2 / h + 0.5 * g * h**2, qx * qy / h, 0.0],
        },
        {
            "velocities": vitesse,
            "conserved_moments": qy,
            "polynomials": polynomes,
            "relaxation_parameters": s_q,
            "equilibrium": [qy, qy * qx / h, qy**2 / h + 0.5 * g * h**2, 0.0],
        },
    ],
    "parameters": {LA: la, g: 1.0},
}

s = pylbm.Scheme(d)
print(s)

Three elementary schemes have been built, these three schemes are identical except for the equilibrium values of the non conserved moment and of the relaxation parameter: The first one is used to simulate the equation on \(h\) and the others to simulate the equation on \(q_x\) and \(q_y\). For each scheme, the equilibrium value of the non conserved moment is equal to the flux of the corresponding equation: the equilibrium value of the kth scheme can so depend on all the conserved moments (and not only on those of the kth scheme).