"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. \n",
"\n",
"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.\n",
"\n",
"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."
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
$$K(f) = C^{-1}\left( C(f) - S (C(f) - C(f^{(eq})) \right)$$
%% Cell type:markdown id: tags:
## 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)})$$
%% Cell type:markdown id: tags:
### 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.
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x10300dfd0>
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x107ac7a30>
%% Cell type:markdown id: tags:
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(c_q)$$
where $m_i(c_q)$ is the $i$'th moment polynomial where $x$ and $y$ are substituted with the components of the $q$'th lattice velocity
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.
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.
%% Cell type:markdown id: tags:
### Central moment lattice Boltzmann methods
%% Cell type:markdown id: tags:
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.
<lbmpy.methods.momentbased.centralmomentbasedmethod.CentralMomentBasedLbMethod at 0x12b422880>
<lbmpy.methods.momentbased.centralmomentbasedmethod.CentralMomentBasedLbMethod at 0x13ad7e460>
%% Cell type:markdown id: tags:
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.
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x12c79daf0>
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x13b22f5b0>
%% Cell type:markdown id: tags:
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.
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x12c5bd670>
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x13b212460>
%% Cell type:markdown id: tags:
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`.
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):
"<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x103b562b0>"
"<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x103b6d610>"
]
},
"execution_count": 9,
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
```
%% Cell type:markdown id: tags:
# Tutorial 04: The Cumulant Lattice Boltzmann Method in lbmpy
## A) Principles of the centered cumulant collision operator
Recently, an advanced Lattice Boltzmann collision operator based on relaxation in cumulant space has gained popularity. Similar to moments, cumulants are statistical quantities inherent to a probability distribution. A significant advantage of the cumulants is that they are statistically independent by construction. Moments can be defined by using the so-called moment generating function, which for the discrete particle distribution of the LB equation is stated as
$$
\begin{align}
M( \vec{X} ) =
\sum_i f_i
\exp \left(
\vec{c}_i \cdot \vec{X}
\right)
\end{align}
$$
The raw moments $m_{\alpha \beta \gamma}$ can be expressed as its derivatives, evaluated at zero:
$$
\begin{align}
m_{\alpha \beta \gamma}
=
\left.
\frac{\partial^{\alpha}}{\partial X^{\alpha}}
\frac{\partial^{\beta}}{\partial Y^{\beta}}
\frac{\partial^{\gamma}}{\partial Z^{\gamma}}
M(X, Y, Z)
\right\vert_{\vec{X}=0}
\end{align}
$$
The cumulant-generating function is defined as the natural logarithm of this moment-generating function, and the cumulants $c_{\alpha \beta \gamma}$ are defined as its derivatives evaluated at zero:
$$
\begin{align}
C(\vec{X}) :=& \, \log ( M(\vec{X}) ) \\
c_{\alpha \beta \gamma}
=& \,
\left.
\frac{\partial^{\alpha}}{\partial X^{\alpha}}
\frac{\partial^{\beta}}{\partial Y^{\beta}}
\frac{\partial^{\gamma}}{\partial Z^{\gamma}}
C(X, Y, Z)
\right\vert_{\vec{X}=0}
\end{align}
$$
Other than with moments, there is no straightforward way to express cumulants in terms of the populations. However, their generating functions can be related to allowing the computation of cumulants from both raw and central moments, computed from populations. In lbmpy, the transformations from populations to cumulants and back are implemented using central moments as intermediaries. This is done for two primary reasons:
1. All cumulants of orders 2 and 3 are equal to their corresponding central moments, up to the density $\rho$ as a proportionality factor.
2. The conserved modes of the first order, which correspond to momentum, are relaxed in central moment space to allow for a more efficient implicit forcing scheme.
The central moment-generating function $K$ can be related to the moment-generating function through $K(\vec{X}) = \exp( - \vec{X} \cdot \vec{u} ) M(\vec{X})$. It is possible to recombine the equation with the definition of the cumulant-generating function
Derivatives of $C$ can thus be expressed in terms of derivatives of $K$, directly yielding equations of the cumulants in terms of central moments.
In the cumulant collision operator, forces can be applied through a simple forcing scheme which in lbmpy is called *implicit forcing*. Ony the force contributions to the first-order (momentum) modes are considered. When the first-order central moments are taken in the frame of reference shifted by $F/2$, they trail behind zero by just that much. Forcing can be applied by first adding half the force *before* collision, and then adding the remaining half *after* collision. Due to this, the first-order central moments entering the cumulant equations are always zero. By applying the force already in central-moment space, the cumulant equations are simplified significantly as no forces need to be taken into account. The pre- and post-collision momentum central moments take the form:
In total, through forcing, the first central moments change sign. This is equivalent to relaxation with a relaxation rate $\omega = 2$. For this reason, lbmpy's implementation of the cumulant LBM calculates the collision of the momentum modes in central moment space, and the default force model overrides their relaxation rate by setting it to 2.
%% Cell type:markdown id: tags:
## B) Method Creation
Using the `create_lb_method` interface, creation of a cumulant-based lattice Boltzmann method in lbmpy is straightforward. Cumulants can either be relaxed in their raw (monomial forms) or in polynomial combinations. Both variants are available as predefined setups.
### Relaxation of Monomial Cumulants
Monomial cumulant relaxation is available through the `method="monomial_cumulant"` parameter setting. This variant requires fewer transformation steps and is a little faster than polynomial relaxation, but it does not allow separation of bulk and shear viscosity. Default monomial cumulant sets are available for the D2Q9, D3Q19 and D3Q27 stencils. It is also possible to define a custom set of monomial cumulants.
When creating a monomial cumulant method, one option is to specify only a single relaxation rate which will then be assigned to all cumulants related to the shear viscosity. In this case, all other non-conserved cumulants will be relaxed to equilibrium. Alternatively, individual relaxation rates for all non-conserved cumulants can be specified. The conserved cumulants are set to zero by default to save computational cost. They can be adjusted with `set_zeroth_moment_relaxation_rate`, `set_first_moment_relaxation_rate` and `set_conserved_moments_relaxation_rate`.
<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x103b14880>
<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x127652c40>
%% Cell type:markdown id: tags:
### Relaxation of Polynomial Cumulants
By setting `method="cumulant"`, a set of default polynomial cumulants is chosen to be relaxed. Those cumulants are taken from literature and assembled into groups selected to enforce rotational invariance (see: `lbmpy.methods.centeredcumulant.centered_cumulants`). Default polynomial groups are available for the D2Q9, D3Q19 and D3Q27 stencils. As before it is possible to specify only a single relaxation rate assigned to the moments governing the shear viscosity. All other relaxation rates are then automatically set to one.
<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x103b14a30>
<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x12745b5e0>
%% Cell type:markdown id: tags:
### Central Moments and Forcing
The conserved modes are marked with the note *(central moment)* in the table above. It highlights the fact that these modes are relaxed in central moment space, other than cumulant space. As described in section A, this is done to enable the implicit forcing scheme. When a force is specified, the momentum modes' relaxation rates are overridden by the force model. In the following cell, a symbolic force is specified. Further, a full list of relaxation rates is passed to allow the specification of bulk viscosity.
<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x11f7fd070>
<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x127652520>
%% Cell type:markdown id: tags:
## C) Exemplary simulation: flow around sphere
To end this tutorial, we show an example simulation of a channel flow with a spherical obstacle. This example is shown in 2D with the D2Q9 stencil but can be adapted easily for a 3D simulation.
%% Cell type:code id: tags:
``` python
from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
```
%% Cell type:markdown id: tags:
To define the channel flow with dimensionless parameters, we first define the reference length in lattice cells. The reference length will be the diameter of the spherical obstacle. Furthermore, we define a maximal velocity which will be set for the inflow later. The Reynolds number is set relatively high to 100000, which will cause a turbulent flow in the channel.
From these definitions, we can calculate the relaxation rate `omega` for the lattice Boltzmann method.
Now the data for the simulation is allocated. We allocate a field `src` for the PDFs and field `dst` used as a temporary field to implement the two grid pull pattern. Additionally, we allocate a velocity field `velField`
We choose a cumulant lattice Boltzmann method, as described above. Here the second-order cumulants are relaxed with the relaxation rate calculated above. All higher-order cumulants are relaxed with one, which means we set them to the equilibrium.
