In :
from lbmpy.session import *
from lbmpy.cumulants import *
from lbmpy.moments import *
from lbmpy.continuous_distribution_measures import continuous_moment
from lbmpy.stencils import get_stencil


# Demo: Moments, Cumulants and Maxwellian Equilibrium¶

## 1) Moments & Cumulants¶

### Moments¶

The moments and cumulants modules contain functions to calculate moments and cumulants of functions on a discrete velocity space defined by a stencil.

In :
stencil = get_stencil("D2Q9")
pdfs = sp.symbols("f:9")
pdfs

Out:
$\displaystyle \left( f_{0}, \ f_{1}, \ f_{2}, \ f_{3}, \ f_{4}, \ f_{5}, \ f_{6}, \ f_{7}, \ f_{8}\right)$

Discrete moments are computed by following formula:

$$\sum_{d \in S} d_1^{m_1} d_2^{m_2} \; f_d$$

with $S$ being the stencil, $d_i$ the direction components, $f_d$ the function values for each direction and $m_j$ the components of the moment tuple.

Lets compute the first moment in the first direction, i.e. $(m_1, m_2) = (1,0)$.

In :
plt.figure(figsize=(3, 3))
ps.stencil.plot(stencil) In :
discrete_moment(pdfs, (1, 0), stencil)

Out:
$\displaystyle - f_{3} + f_{4} - f_{5} + f_{6} - f_{7} + f_{8}$

We get contributions of all directions that have a non-zero x-component, weighted with the sign of the direction.

For the second order moment, the direction components are squared, so all contributions come in with a positive sign.

In :
discrete_moment(pdfs, (2, 0), stencil)

Out:
$\displaystyle f_{3} + f_{4} + f_{5} + f_{6} + f_{7} + f_{8}$

We can specify the moments not only as exponent tuples but also as polynomials in the $d_i$'s. The symbols for the $d_i$ are provided by the moments module as MOMENT_SYMBOLS

In :
x, y, z = MOMENT_SYMBOLS
discrete_moment(pdfs, x, stencil)

Out:
$\displaystyle - f_{3} + f_{4} - f_{5} + f_{6} - f_{7} + f_{8}$

The same works also for sums:

In :
discrete_moment(pdfs, x**2 * y + y**2, stencil)

Out:
$\displaystyle f_{1} + f_{2} + 2 f_{5} + 2 f_{6}$

Here the advantage of the polynomial representation becomes visible. To compute the same moment with exponent tuples takes two calls:

In :
discrete_moment(pdfs, (2, 1), stencil) + discrete_moment(pdfs, (0, 2), stencil)

Out:
$\displaystyle f_{1} + f_{2} + 2 f_{5} + 2 f_{6}$

### Cumulants¶

Cumulants are an alternative to the moment represenation, see the Wikipedia article. Cumulants can be calculated directly using the cumulant generating function, or can be calculated by moments.

In :
discrete_cumulant(pdfs, (2, 0), stencil)

Out:
$\displaystyle \frac{f_{3} + f_{4} + f_{5} + f_{6} + f_{7} + f_{8} - \frac{\left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8}\right)^{2}}{f_{0} + f_{1} + f_{2} + f_{3} + f_{4} + f_{5} + f_{6} + f_{7} + f_{8}}}{f_{0} + f_{1} + f_{2} + f_{3} + f_{4} + f_{5} + f_{6} + f_{7} + f_{8}}$
In :
cumulant_as_function_of_raw_moments((2,0))

Out:
$\displaystyle \frac{m_{2 0}}{m_{0 0}} - \frac{m_{1 0}^{2}}{m_{0 0}^{2}}$
In :
raw_moment_as_function_of_cumulants((2,0))

Out:
$\displaystyle c_{1 0}^{2} e^{c_{0 0}} + c_{2 0} e^{c_{0 0}}$

## 2) Using Moments to derive discrete LBM equilibrium¶

For full stencils i.e. D2Q9 and D3Q27 a LBM equilibrium can be derived using the following strategy:

• calculate all moments or cumulants up to second order of the continuous Maxwell-Boltzmann distribution
• for full stencils there are as many moments/cumulants as stencil directions, so we can determine the pdf values from them

First we obtain the continuous Maxwellian equilibrium as sympy expression. We fix the speed of sound to $\frac{1}{3}$

In :
from lbmpy.maxwellian_equilibrium import continuous_maxwellian_equilibrium, discrete_maxwellian_equilibrium

dim = len(stencil)
maxwellian = continuous_maxwellian_equilibrium(dim=dim, c_s_sq=sp.Rational(1,3))
maxwellian

Out:
$\displaystyle \frac{3 \rho e^{- \frac{3 \left(- u_{0} + v_{0}\right)^{2}}{2} - \frac{3 \left(- u_{1} + v_{1}\right)^{2}}{2}}}{2 \pi}$

Then we get all moment exponent tuples up to order 2 and calculate these moments of the continuous Maxwellian equilibrium:

In :
moments = moments_up_to_component_order(2, dim=dim)
moments

Out:
$\displaystyle \left( \left( 0, \ 0\right), \ \left( 0, \ 1\right), \ \left( 0, \ 2\right), \ \left( 1, \ 0\right), \ \left( 1, \ 1\right), \ \left( 1, \ 2\right), \ \left( 2, \ 0\right), \ \left( 2, \ 1\right), \ \left( 2, \ 2\right)\right)$
In :
cont_eq_moments = [continuous_moment(maxwellian, m, sp.symbols("v_:2")[:dim]) for m in moments]
cont_eq_moments

Out:
$\displaystyle \left[ \rho, \ \rho u_{1}, \ \frac{\rho \left(9 u_{1}^{2} + 3\right)}{9}, \ \rho u_{0}, \ \rho u_{0} u_{1}, \ \frac{\rho u_{0} \left(9 u_{1}^{2} + 3\right)}{9}, \ \frac{\rho \left(9 u_{0}^{2} + 3\right)}{9}, \ \frac{\rho u_{1} \left(9 u_{0}^{2} + 3\right)}{9}, \ \frac{\rho \left(9 u_{0}^{2} + 3\right) \left(9 u_{1}^{2} + 3\right)}{81}\right]$

To obtain the same equilibrium as in the LBM literature, we have to drop all terms that have order 3 or higher:

In :
cont_eq_moments = [remove_higher_order_terms(m, order=3, symbols=sp.symbols("u_:3"))
for m in cont_eq_moments]
cont_eq_moments

Out:
$\displaystyle \left[ \rho, \ \rho u_{1}, \ \rho u_{1}^{2} + \frac{\rho}{3}, \ \rho u_{0}, \ \rho u_{0} u_{1}, \ \rho u_{0} u_{1}^{2} + \frac{\rho u_{0}}{3}, \ \rho u_{0}^{2} + \frac{\rho}{3}, \ \rho u_{0}^{2} u_{1} + \frac{\rho u_{1}}{3}, \ \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right]$

Then we take these equilibrium moments and determine pdf values, which would lead to these moment values. The moment matrix transforms pdfs into moment space. To obtain the equilibrium pdfs we only have to transform the equilibrium moments with the inverse of this matrix:

In :
M = moment_matrix(moments, stencil)
M

Out:
$\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1\\0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1\\0 & 0 & 0 & -1 & 1 & -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 & 1 & 1 & 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]$
In :
derived_eq = M.inv() * sp.Matrix(cont_eq_moments)
derived_eq

Out:
$\displaystyle \left[\begin{matrix}- \frac{2 \rho u_{0}^{2}}{3} - \frac{2 \rho u_{1}^{2}}{3} + \frac{4 \rho}{9}\\- \frac{\rho u_{0}^{2} u_{1}}{2} - \frac{\rho u_{0}^{2}}{6} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho u_{1}}{3} + \frac{\rho}{9}\\\frac{\rho u_{0}^{2} u_{1}}{2} - \frac{\rho u_{0}^{2}}{6} + \frac{\rho u_{1}^{2}}{3} - \frac{\rho u_{1}}{3} + \frac{\rho}{9}\\\frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{0} u_{1}^{2}}{2} - \frac{\rho u_{0}}{3} - \frac{\rho u_{1}^{2}}{6} + \frac{\rho}{9}\\\frac{\rho u_{0}^{2}}{3} - \frac{\rho u_{0} u_{1}^{2}}{2} + \frac{\rho u_{0}}{3} - \frac{\rho u_{1}^{2}}{6} + \frac{\rho}{9}\\\frac{\rho u_{0}^{2} u_{1}}{4} + \frac{\rho u_{0}^{2}}{12} - \frac{\rho u_{0} u_{1}^{2}}{4} - \frac{\rho u_{0} u_{1}}{4} - \frac{\rho u_{0}}{12} + \frac{\rho u_{1}^{2}}{12} + \frac{\rho u_{1}}{12} + \frac{\rho}{36}\\\frac{\rho u_{0}^{2} u_{1}}{4} + \frac{\rho u_{0}^{2}}{12} + \frac{\rho u_{0} u_{1}^{2}}{4} + \frac{\rho u_{0} u_{1}}{4} + \frac{\rho u_{0}}{12} + \frac{\rho u_{1}^{2}}{12} + \frac{\rho u_{1}}{12} + \frac{\rho}{36}\\- \frac{\rho u_{0}^{2} u_{1}}{4} + \frac{\rho u_{0}^{2}}{12} - \frac{\rho u_{0} u_{1}^{2}}{4} + \frac{\rho u_{0} u_{1}}{4} - \frac{\rho u_{0}}{12} + \frac{\rho u_{1}^{2}}{12} - \frac{\rho u_{1}}{12} + \frac{\rho}{36}\\- \frac{\rho u_{0}^{2} u_{1}}{4} + \frac{\rho u_{0}^{2}}{12} + \frac{\rho u_{0} u_{1}^{2}}{4} - \frac{\rho u_{0} u_{1}}{4} + \frac{\rho u_{0}}{12} + \frac{\rho u_{1}^{2}}{12} - \frac{\rho u_{1}}{12} + \frac{\rho}{36}\end{matrix}\right]$

This is the same as the standard discrete equilibrium found in LBM literature.

In :
literature_version = sp.Matrix(discrete_maxwellian_equilibrium(stencil, c_s_sq=sp.Rational(1,3), order=3))
literature_version

Out:
$\displaystyle \left[\begin{matrix}- \frac{2 \rho u_{0}^{2}}{3} - \frac{2 \rho u_{1}^{2}}{3} + \frac{4 \rho}{9}\\- \frac{\rho u_{0}^{2} u_{1}}{2} - \frac{\rho u_{0}^{2}}{6} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho u_{1}}{3} + \frac{\rho}{9}\\\frac{\rho u_{0}^{2} u_{1}}{2} - \frac{\rho u_{0}^{2}}{6} + \frac{\rho u_{1}^{2}}{3} - \frac{\rho u_{1}}{3} + \frac{\rho}{9}\\\frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{0} u_{1}^{2}}{2} - \frac{\rho u_{0}}{3} - \frac{\rho u_{1}^{2}}{6} + \frac{\rho}{9}\\\frac{\rho u_{0}^{2}}{3} - \frac{\rho u_{0} u_{1}^{2}}{2} + \frac{\rho u_{0}}{3} - \frac{\rho u_{1}^{2}}{6} + \frac{\rho}{9}\\\frac{\rho u_{0}^{2} u_{1}}{4} + \frac{\rho u_{0}^{2}}{12} - \frac{\rho u_{0} u_{1}^{2}}{4} - \frac{\rho u_{0} u_{1}}{4} - \frac{\rho u_{0}}{12} + \frac{\rho u_{1}^{2}}{12} + \frac{\rho u_{1}}{12} + \frac{\rho}{36}\\\frac{\rho u_{0}^{2} u_{1}}{4} + \frac{\rho u_{0}^{2}}{12} + \frac{\rho u_{0} u_{1}^{2}}{4} + \frac{\rho u_{0} u_{1}}{4} + \frac{\rho u_{0}}{12} + \frac{\rho u_{1}^{2}}{12} + \frac{\rho u_{1}}{12} + \frac{\rho}{36}\\- \frac{\rho u_{0}^{2} u_{1}}{4} + \frac{\rho u_{0}^{2}}{12} - \frac{\rho u_{0} u_{1}^{2}}{4} + \frac{\rho u_{0} u_{1}}{4} - \frac{\rho u_{0}}{12} + \frac{\rho u_{1}^{2}}{12} - \frac{\rho u_{1}}{12} + \frac{\rho}{36}\\- \frac{\rho u_{0}^{2} u_{1}}{4} + \frac{\rho u_{0}^{2}}{12} + \frac{\rho u_{0} u_{1}^{2}}{4} - \frac{\rho u_{0} u_{1}}{4} + \frac{\rho u_{0}}{12} + \frac{\rho u_{1}^{2}}{12} - \frac{\rho u_{1}}{12} + \frac{\rho}{36}\end{matrix}\right]$

## 3) Reduced Stencils¶

This method for deriving a discrete equilibrium works well for "full stencils" i.e. with $3^d$ directions, where $d$ is the dimension.

In [ ]:
import pytest
pytest.importorskip('ipy_table')

In :
moment_equality_table(get_stencil("D2Q9"), truncate_order=2)

Matched moments 13 - non matched moments 2 - total 15

Out:
 order 0 (0, 0)  x 1 1 (1, 0)  x 2 2 (1, 1)  x 1 (2, 0)  x 2 3 (2, 1)  x 2 (3, 0)  x 2 4 (2, 2)  x 1 (3, 1)  x 2 (4, 0)  x 2
In :
moment_equality_table(get_stencil("D3Q27"), truncate_order=2)

Matched moments 32 - non matched moments 3 - total 35

Out:
 order 0 (0, 0, 0)  x 1 1 (1, 0, 0)  x 3 2 (1, 1, 0)  x 3 (2, 0, 0)  x 3 3 (1, 1, 1)  x 1 (2, 1, 0)  x 6 (3, 0, 0)  x 3 4 (2, 1, 1)  x 3 (2, 2, 0)  x 3 (3, 1, 0)  x 6 (4, 0, 0)  x 3
In :
moment_equality_table(get_stencil("D3Q19"), truncate_order=2)

Matched moments 26 - non matched moments 9 - total 35

Out:
 order 0 (0, 0, 0)  x 1 1 (1, 0, 0)  x 3 2 (1, 1, 0)  x 3 (2, 0, 0)  x 3 3 (1, 1, 1)  x 1 (2, 1, 0)  x 6 (3, 0, 0)  x 3 4 (2, 1, 1)  x 3 (2, 2, 0)  x 3 (3, 1, 0)  x 6 (4, 0, 0)  x 3