[1]:

from lbmpy.session import *


# Tutorial 03: Defining LB methods in lbmpy¶

## A) General Form¶

The lattice Boltzmann equation in its most general form is:

$f_q(\mathbf{x} + \mathbf{c}_q \delta t, t+\delta t) = K\left( f_q(\mathbf{x}, t) \right)$

with a discrete velocity set $$\mathbf{c}_q$$ (stencil) and a generic collision operator $$K$$.

So a lattice Boltzmann method can be fully defined by picking a stencil and a collision operator. The collision operator $$K$$ has the following structure: - Transformation of particle distribution function $$f$$ into collision space. This transformation has to be invertible and may be nonlinear. - The collision operation is an convex combination of the pdf representation in collision space $$c$$ and some equilibrium vector $$c^{(eq)}$$. This equilibrium can also be defined in physical space, then $c^{(eq)} = C( f^{(eq)} )$. The convex combination is done elementwise using a diagonal matrix $$S$$ where the diagonal entries are the relaxation rates. - After collision, the collided state $$c'$$ is transformed back into physical space

The full collision operator is:

$K(f) = C^{-1}\left( (I-S)C(f) + SC(f^{(eq}) \right)$

or

$K(f) = C^{-1}\left( C(f) - S (C(f) - C(f^{(eq})) \right)$

## B) Moment-based relaxation¶

The most commonly used LBM collision operator is the multi relaxation time (MRT) collision. In MRT methods the collision space is spanned by moments of the distribution function. This is a very natural approach, since the pdf moments are the quantities that go into the Chapman Enskog analysis that is used to show that LB methods can solve the Navier Stokes equations. Also the lower order moments correspond to the macroscopic quantities of interest (density/pressure, velocity, shear rates, heat flux). Furthermore the transformation to collision space is linear in this case, simplifying the collision equations:

$K(f) = C^{-1}\left( C(f) - S (C(f) - C(f^{(eq})) \right)$
$K(f) = f - \underbrace{ C^{-1}SC}_{A}(f - f^{(eq)})$

in lbmpy the following formulation is used, since it is more natural to define the equilibrium in moment-space instead of physical space:

$K(f) = f - C^{-1}S(Cf - c^{(eq)})$

### Use a pre-defined method¶

Lets create a moment-based method in lbmpy and see how the moment transformation $$C$$ and the relaxation rates that comprise the diagonal matrix $$S$$ can be defined. We start with a function that creates a basic MRT model. Don’t use this for real simulations, there orthogonalized MRT methods should be used, as discussed in the next section.

[2]:

from lbmpy.creationfunctions import create_lb_method
lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT_RAW)

method = create_lb_method(lbm_config=lbm_config)
# check also method='srt', 'trt', 'mrt'
method

[2]:

Moment Eq. Value Relaxation Rate
$1$ $\rho$ $\omega$
$x$ $u_{0}$ $\omega$
$y$ $u_{1}$ $\omega$
$x^{2}$ $\frac{\rho}{3} + u_{0}^{2}$ $\omega$
$y^{2}$ $\frac{\rho}{3} + u_{1}^{2}$ $\omega$
$x y$ $u_{0} u_{1}$ $\omega$
$x^{2} y$ $\frac{u_{1}}{3}$ $\omega$
$x y^{2}$ $\frac{u_{0}}{3}$ $\omega$
$x^{2} y^{2}$ $\frac{\rho}{9} + \frac{u_{0}^{2}}{3} + \frac{u_{1}^{2}}{3}$ $\omega$

The first column labeled “Moment” defines the collision space and thus the transformation matrix $$C$$. The remaining columns specify the equilibrium vector in moment space $$c^{(eq)}$$ and the corresponding relaxation rate.

Each row of the “Moment” column defines one row of $$C$$. In the next cells this matrix and the discrete velocity set (stencil) of our method are shown. Check for example the second last row of the table $$x^2 y$$: In the corresponding second last row of the moment matrix $$C$$ where each column stands for a lattice velocity (for ordering visualized stencil below) and each entry is the expression $$x^2 y$$ where $$x$$ and $$y$$ are the components of the lattice velocity.

In general the transformation matrix $$C_{iq}$$ is defined as;

$c_i = C_{iq} f_q = \sum_q m_i(\mathbf{c}_q)$

where $$m_i(\mathbf{c}_q)$$ is the $$i$$’th moment polynomial where $$x$$ and $$y$$ are substituted with the components of the $$q$$’th lattice velocity

[3]:

# Transformation matrix C
method.moment_matrix

[3]:

$\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\0 & 0 & 0 & -1 & 1 & -1 & 1 & -1 & 1\\0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1\\0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1\\0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1\\0 & 0 & 0 & 0 & 0 & -1 & 1 & 1 & -1\\0 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1\\0 & 0 & 0 & 0 & 0 & -1 & 1 & -1 & 1\\0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\end{matrix}\right]$
[4]:

method.stencil.plot()

[5]:

method.stencil

[5]:

Direction Name Direction
$\mathtt{\text{C}}$ $\left( 0, \ 0\right)$
$\mathtt{\text{N}}$ $\left( 0, \ 1\right)$
$\mathtt{\text{S}}$ $\left( 0, \ -1\right)$
$\mathtt{\text{W}}$ $\left( -1, \ 0\right)$
$\mathtt{\text{E}}$ $\left( 1, \ 0\right)$
$\mathtt{\text{NW}}$ $\left( -1, \ 1\right)$
$\mathtt{\text{NE}}$ $\left( 1, \ 1\right)$
$\mathtt{\text{SW}}$ $\left( -1, \ -1\right)$
$\mathtt{\text{SE}}$ $\left( 1, \ -1\right)$

### Orthogonal MRTs¶

For a real MRT method, the moments should be orthogonalized. One can either orthogonalize using the standard scalar product or a scalar product that is weighted with the lattice weights. If unsure, use the weighted version.

The next cell shows how to get both orthogonalized MRT versions in lbmpy.

[6]:

rr = [sp.Symbol('omega_shear'), sp.Symbol('omega_bulk'), sp.Symbol('omega_3'), sp.Symbol('omega_4')]

lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, weighted=True,
relaxation_rates=rr)
weighted_ortho_mrt = create_lb_method(lbm_config=lbm_config)
weighted_ortho_mrt

[6]:

Moment Eq. Value Relaxation Rate
$1$ $\rho$ $0$
$x$ $u_{0}$ $0$
$y$ $u_{1}$ $0$
$x^{2} - y^{2}$ $u_{0}^{2} - u_{1}^{2}$ $\omega_{shear}$
$x y$ $u_{0} u_{1}$ $\omega_{shear}$
$3 x^{2} + 3 y^{2} - 2$ $3 u_{0}^{2} + 3 u_{1}^{2}$ $\omega_{bulk}$
$3 x^{2} y - y$ $0$ $\omega_{3}$
$3 x y^{2} - x$ $0$ $\omega_{3}$
$9 x^{2} y^{2} - 3 x^{2} - 3 y^{2} + 1$ $0$ $\omega_{4}$
[7]:

lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, weighted=False,
relaxation_rates=rr)
ortho_mrt = create_lb_method(lbm_config=lbm_config)
ortho_mrt

[7]:

Moment Eq. Value Relaxation Rate
$1$ $\rho$ $0$
$x$ $u_{0}$ $0$
$y$ $u_{1}$ $0$
$x^{2} - y^{2}$ $u_{0}^{2} - u_{1}^{2}$ $\omega_{shear}$
$x y$ $u_{0} u_{1}$ $\omega_{shear}$
$3 x^{2} + 3 y^{2} - 4$ $- 2 \rho + 3 u_{0}^{2} + 3 u_{1}^{2}$ $\omega_{bulk}$
$3 x^{2} y - 2 y$ $- u_{1}$ $\omega_{3}$
$3 x y^{2} - 2 x$ $- u_{0}$ $\omega_{3}$
$9 x^{2} y^{2} - 6 x^{2} - 6 y^{2} + 4$ $\rho - 3 u_{0}^{2} - 3 u_{1}^{2}$ $\omega_{4}$

One can check if a method is orthogonalized:

[8]:

ortho_mrt.is_orthogonal, weighted_ortho_mrt.is_weighted_orthogonal

[8]:

(True, True)


Note here how the relaxation rate for each moment is set. By default, a single relaxation rate per moment group is assigned. With a moment group, we mean here moment polynomials of the same order. However, the second-order moments are split between the shear and bulk moments (defining the shear and bulk viscosity). Thus the first two relaxation rates in the rr-list are used for these.

If only a single relaxation rate would be defined, it is applied to the shear moments, and all other higher-order moments are set to one. This means the PDFs are directly fixed to the equilibrium in the collision process. Furthermore, the relaxation rate for conserved moments can be chosen freely as they are conserved in the collision. Thus by default, we set it to zero.

It is also possible to define an individual relaxation rate per moment. In this case, Q relaxation rates need to be defined, where Q is the number of moments.

### Central moment lattice Boltzmann methods¶

Another popular method is the cascaded lattice Boltzmann method. The cascaded LBM increases the numerical stability by shifting the collision step to the central moment space. Thus it is applied in the non-moving frame and achieves a better Galilean invariance. Typically the central moment collision operator is derived for compressible LB methods, and a higher-order equilibrium is used. Although incompressible LB methods with a second-order equilibrium can be derived with lbmpy, it violates the Galilean invariance and thus reduces the advantages of the method.

[9]:

lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.CENTRAL_MOMENT, equilibrium_order=4,
compressible=True, relaxation_rates=rr)

central_moment_method = create_lb_method(lbm_config)
central_moment_method

[9]:

Central Moment Eq. Value Relaxation Rate
$1$ $\rho$ $0$
$x$ $0$ $0$
$y$ $0$ $0$
$x y$ $0$ $\omega_{shear}$
$x^{2} - y^{2}$ $0$ $\omega_{shear}$
$x^{2} + y^{2}$ $\frac{2 \rho}{3}$ $\omega_{bulk}$
$x^{2} y$ $0$ $\omega_{3}$
$x y^{2}$ $0$ $\omega_{3}$
$x^{2} y^{2}$ $\frac{\rho}{9}$ $\omega_{4}$

The shift to the central moment space is done by applying a so-called shift matrix. Usually, this introduces a high numerical overhead. This problem is solved with lbmpy because each transformation stage can be specifically optimised individually. Therefore, it is possible to derive a CLBM with only a little numerical overhead.

[10]:

central_moment_method.shift_matrix

[10]:

$\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- u_{0} & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- u_{1} & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\u_{0} u_{1} & - u_{1} & - u_{0} & 1 & 0 & 0 & 0 & 0 & 0\\u_{0}^{2} - u_{1}^{2} & - 2 u_{0} & 2 u_{1} & 0 & 1 & 0 & 0 & 0 & 0\\u_{0}^{2} + u_{1}^{2} & - 2 u_{0} & - 2 u_{1} & 0 & 0 & 1 & 0 & 0 & 0\\- u_{0}^{2} u_{1} & 2 u_{0} u_{1} & u_{0}^{2} & - 2 u_{0} & - \frac{u_{1}}{2} & - \frac{u_{1}}{2} & 1 & 0 & 0\\- u_{0} u_{1}^{2} & u_{1}^{2} & 2 u_{0} u_{1} & - 2 u_{1} & \frac{u_{0}}{2} & - \frac{u_{0}}{2} & 0 & 1 & 0\\u_{0}^{2} u_{1}^{2} & - 2 u_{0} u_{1}^{2} & - 2 u_{0}^{2} u_{1} & 4 u_{0} u_{1} & - \frac{u_{0}^{2}}{2} + \frac{u_{1}^{2}}{2} & \frac{u_{0}^{2}}{2} + \frac{u_{1}^{2}}{2} & - 2 u_{1} & - 2 u_{0} & 1\end{matrix}\right]$

### Define custom MRT method¶

To choose custom values for the left moment column one can pass a nested list of moments. Moments that should be relaxed with the same paramter are grouped together.

lbmpy also comes with a few templates for this list taken from literature:

[11]:

from lbmpy.methods import mrt_orthogonal_modes_literature
from lbmpy.stencils import get_stencil
from lbmpy.moments import MOMENT_SYMBOLS

x, y, z = MOMENT_SYMBOLS

moments = mrt_orthogonal_modes_literature(LBStencil(Stencil.D2Q9), is_weighted=True)
moments

[11]:

$\displaystyle \left[ \left[ 1\right], \ \left[ x, \ y\right], \ \left[ 3 x^{2} + 3 y^{2} - 2, \ x^{2} - y^{2}, \ x y\right], \ \left[ x \left(3 x^{2} + 3 y^{2} - 4\right), \ y \left(3 x^{2} + 3 y^{2} - 4\right)\right], \ \left[ - 15 x^{2} - 15 y^{2} + 9 \left(x^{2} + y^{2}\right)^{2} + 2\right]\right]$

This nested moment list can be passed to create_lb_method:

[12]:

create_lb_method(LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, nested_moments=moments,
relaxation_rates=rr))

[12]:

Moment Eq. Value Relaxation Rate
$1$ $\rho$ $0$
$x$ $u_{0}$ $0$
$y$ $u_{1}$ $0$
$3 x^{2} + 3 y^{2} - 2$ $3 u_{0}^{2} + 3 u_{1}^{2}$ $\omega_{bulk}$
$x^{2} - y^{2}$ $u_{0}^{2} - u_{1}^{2}$ $\omega_{shear}$
$x y$ $u_{0} u_{1}$ $\omega_{shear}$
$x \left(3 x^{2} + 3 y^{2} - 4\right)$ $0$ $\omega_{3}$
$y \left(3 x^{2} + 3 y^{2} - 4\right)$ $0$ $\omega_{3}$
$- 15 x^{2} - 15 y^{2} + 9 \left(x^{2} + y^{2}\right)^{2} + 2$ $9 u_{0}^{2} + 9 u_{1}^{2}$ $\omega_{4}$

If one needs to also specify custom equilibrium moments the following approach can be used

[13]:

rho = sp.symbols("rho")
u = sp.symbols("u_:3")
omega = sp.symbols("omega_:4")

method_table = [
# Conserved moments
(1,    rho,  0 ),
(x,    u[0], 0 ),
(y,    u[1], 0 ),

# Shear moments
(x*y,       u[0]*u[1],                   omega[0]),
(x**2-y**2, u[0]**2 - u[1]**2,           omega[0]),
(x**2+y**2, 2*rho/3 + u[0]**2 + u[1]**2, omega[1]),

# Higher order
(x * y**2,    u[0]/3,                        omega[2]),
(x**2 * y,    u[1]/3,                        omega[2]),
(x**2 * y**2, rho/9 + u[0]**2/3 + u[1]**2/3, omega[3]),
]
method = create_generic_mrt(LBStencil(Stencil.D2Q9), method_table)
method

[13]:

Moment Eq. Value Relaxation Rate
$1$ $\rho$ $0$
$x$ $u_{0}$ $0$
$y$ $u_{1}$ $0$
$x y$ $u_{0} u_{1}$ $\omega_{0}$
$x^{2} - y^{2}$ $u_{0}^{2} - u_{1}^{2}$ $\omega_{0}$
$x^{2} + y^{2}$ $\frac{2 \rho}{3} + u_{0}^{2} + u_{1}^{2}$ $\omega_{1}$
$x y^{2}$ $\frac{u_{0}}{3}$ $\omega_{2}$
$x^{2} y$ $\frac{u_{1}}{3}$ $\omega_{2}$
$x^{2} y^{2}$ $\frac{\rho}{9} + \frac{u_{0}^{2}}{3} + \frac{u_{1}^{2}}{3}$ $\omega_{3}$

Instead of manually defining all entries in the method table, lbmpy has functions to fill the table according to a specific pattern. For example: - for a full stencil (D2Q9, D3Q27) there exist exactly 9 or 27 linearly independent moments. These can either be taken as they are, or orthogonalized using Gram-Schmidt, weighted Gram-Schmidt or a Hermite approach - equilibrium values can be computed from the standard discrete equilibrium of order 1,2 or 3. Alternatively they can also be computed as continuous moments of a Maxwellian distribution

One option is to start with one of lbmpy’s built-in methods and modify it with create_lb_method_from_existing. In the next cell we fix the fourth order relaxation rate to a constant, by writing a function that defines how to alter each row of the collision table. This is for demonstration only, of course we could have done it right away when passing in the collision table.

[14]:

def modification_func(moment, eq, rate):
if rate == omega[3]:
return moment, eq, 1.0
return moment, eq, rate

method = create_lb_method_from_existing(method, modification_func)
method

[14]:

Moment Eq. Value Relaxation Rate
$1$ $\rho$ $0$
$x$ $u_{0}$ $0$
$y$ $u_{1}$ $0$
$x y$ $u_{0} u_{1}$ $\omega_{0}$
$x^{2} - y^{2}$ $u_{0}^{2} - u_{1}^{2}$ $\omega_{0}$
$x^{2} + y^{2}$ $\frac{2 \rho}{3} + u_{0}^{2} + u_{1}^{2}$ $\omega_{1}$
$x y^{2}$ $\frac{u_{0}}{3}$ $\omega_{2}$
$x^{2} y$ $\frac{u_{1}}{3}$ $\omega_{2}$
$x^{2} y^{2}$ $\frac{\rho}{9} + \frac{u_{0}^{2}}{3} + \frac{u_{1}^{2}}{3}$ $1.0$

Our customized method can be directly passed into one of the scenarios. We can for example set up a channel flow with it. Since we used symbols as relaxation rates, we have to pass them in as kernel_params.

[15]:

ch = create_channel(domain_size=(100, 30), lb_method=method, u_max=0.05,
kernel_params={'omega_0': 1.8, 'omega_1': 1.4, 'omega_2': 1.5})
ch.run(500)
plt.figure(dpi=200)
plt.vector_field(ch.velocity[:, :]);


#### Bonus: Automatic analysis¶

Above we have created a non-orthogonal MRT, where the shear viscosity and bulk viscosity can be independently controlled. For moment-based methods, lbmpy also offers an automatic Chapman Enskog analysis that can find the relation between viscosity and relaxation rate(s):

[16]:

from lbmpy.chapman_enskog import ChapmanEnskogAnalysis
analysis = ChapmanEnskogAnalysis(method)
analysis.get_dynamic_viscosity()

[16]:

$\displaystyle - \frac{\omega_{0} - 2}{6 \omega_{0}}$
[17]:

analysis.get_bulk_viscosity()

[17]:

$\displaystyle - \frac{1}{9} - \frac{1}{3 \omega_{1}} + \frac{5}{9 \omega_{0}}$