Commit 9437fa74 authored by Markus Holzer's avatar Markus Holzer
Browse files

Finished cumulant tutorial

parent 328deb55
Pipeline #29828 passed with stage
in 23 minutes and 20 seconds
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
```
%% Cell type:markdown id: tags:
# Tutorial 09: 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. For the purpose of Lattice Boltzmann collision, they are superior to moments since they are statistically independent by construction. Moments are definable using the so-called moment generating function, which for the discrete particle distribution of the LB equation is stated as
$$
M( \vec{X} ) =
\sum_i f_i
\exp \left(
\vec{c}_i \cdot \vec{X}
\right)
$$
The raw moments $m_{\alpha \beta \gamma}$ can be expressed as its derivatives, evaluated at zero:
$$
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}
$$
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:
$$
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}
$$
Other than with moments, there is no straightforward way to express cumulants in terms of the populations. However, their generating functions can be related in a way that allows computation of cumulants from both raw and central moments, which in turn are 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 major reasons:
* All cumulants of orders 2 and 3 are equal to their corresponding central moments, up to the density $\rho$ as a proportionality factor
* The conserved modes of first order, which correspond to momentum, are relaxed in central moment space to allow for an especially 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})$. The equation can be recombined with the definition of the cumulant-generating function to obtain
$$
C( \vec{X} ) = \vec{X} \cdot \vec{u} + \log( K( \vec{X} ) ).
$$
Derivatives of $C$ can thus be expressed in terms of derivatives of $K$, directly yielding equations of the cumulants in terms of central moments.
Applying a force in cumulant space is straightforward since all cumulants of the LB forcing term are zero, expect for those related to momentum. The first-order central moments of the force are equal to the force vector. When the velocity entering the collision equilibria is shifted by the half-force shift $F/2$, the first central moments trail behind their frame of reference by exactly this shift. Application of the force only affects the momentum moments:
$$
\kappa_{100} = - \frac{F_x}{2}, \quad
\kappa^{\ast}_{100} = \kappa_{100} + F_x = \frac{F_x}{2}
$$
Basically, 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 simply 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 pre-defined 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. If you'd like to specify a custom set of monomial cumulants, refer to section C.
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.
%% Cell type:code id: tags:
``` python
method_monomial = create_lb_method(method='monomial_cumulant', relaxation_rate=sp.Symbol('omega_v'))
method_monomial
```
%%%% Output: execute_result
<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f45c4f29df0>
%% 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 are assembled into groups selected to enforce rotational invariance (cf. :cite:`Geier2006`, :cite:`geier2015`, :cite:`Coreixas2019`). Default polynomial groups are available for the D2Q9, D3Q19 and D3Q27 stencils. It is again possible to specify only a single relaxation rate which is assigned to the moments governing the shear viscosity. All other relaxation rates are then automatically set to one.
%% Cell type:code id: tags:
``` python
method_polynomial = create_lb_method(method='cumulant', relaxation_rate=sp.Symbol('omega_v'))
method_polynomial
```
%%%% Output: execute_result
<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f454f02bb20>
%% Cell type:markdown id: tags:
### Central Moments and Forcing
In above printouts, the conserved modes are marked with the note *(central moment)*. This is to highlight 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 relaxation rates of the momentum modes 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 specification of a bulk viscosity.
%% Cell type:code id: tags:
``` python
method_params = {
'method' : 'cumulant',
'force' : sp.symbols('F_:3'),
'relaxation_rates' : [sp.Symbol('omega_shear'), sp.Symbol('omega_bulk'), 1, 1]
}
method_with_force = create_lb_method(**method_params)
method_with_force
```
%%%% Output: execute_result
<lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x7f454cea04c0>
%% Cell type:markdown id: tags:
## C) Advanced Options
A lot is going on behind the scenes while the collision equations are derived.
This diff is collapsed.
......@@ -358,7 +358,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
"version": "3.8.2"
}
},
"nbformat": 4,
......
......@@ -393,7 +393,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.9"
"version": "3.8.2"
}
},
"nbformat": 4,
......
......@@ -182,7 +182,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
"version": "3.8.2"
}
},
"nbformat": 4,
......
......@@ -603,7 +603,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
"version": "3.8.2"
}
},
"nbformat": 4,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment