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.

stencil = LBStencil(Stencil.D2Q9)
pdfs = sp.symbols("f:9")
pdfs

$\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)\).

plt.figure(figsize=(3, 3))
stencil.plot()

discrete_moment(pdfs, (1, 0), stencil)

$\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.

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

$\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

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

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

The same works also for sums:

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

$\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:

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

$\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.

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

$\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}}$

cumulant_as_function_of_raw_moments((2,0))

$\displaystyle \frac{m_{2 0}}{m_{0 0}} - \frac{m_{1 0}^{2}}{m_{0 0}^{2}}$

raw_moment_as_function_of_cumulants((2,0))

$\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}\)

from lbmpy.maxwellian_equilibrium import continuous_maxwellian_equilibrium, discrete_maxwellian_equilibrium

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

$\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:

moments = moments_up_to_component_order(2, dim=dim)
moments

$\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)$

cont_eq_moments = [continuous_moment(maxwellian, m, sp.symbols("v_:2")[:dim]) for m in moments]
cont_eq_moments

$\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:

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

$\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:

M = moment_matrix(moments, stencil)
M

$\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]$

derived_eq = M.inv() * sp.Matrix(cont_eq_moments)
derived_eq

$\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.

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

$\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.

import pytest
pytest.importorskip('ipy_table')

<module 'ipy_table' from '/home/markus/miniconda3/envs/pystencils/lib/python3.8/site-packages/ipy_table/__init__.py'>

moment_equality_table(LBStencil(Stencil.D2Q9), truncate_order=2)

Matched moments 13 - non matched moments 2 - total 15

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

moment_equality_table(LBStencil(Stencil.D3Q19), truncate_order=2)

Matched moments 26 - non matched moments 9 - total 35

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

moment_equality_table(LBStencil(Stencil.D3Q19), truncate_order=2)

Matched moments 26 - non matched moments 9 - total 35

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