[@ggruszczynski](https://github.com/ggruszczynski)

# LBM - some theory

In this tutorial you will get the idea of:

* symbolic code  generation
* a bit more advanced LBM models 


## Moments of Equilibrium Distribution Function

We are going to show how the moments of equilibrium distribution function can be calculated.

The formulas for the discrete equilibrium distribution function $ f^{eq}_i $
comes from a discretization of the continous Maxwell-Boltzmann distribution function.
The Maxwell-Boltzmann equilibrium distribution function in a continuous, velocity space is known as:

$$
\Psi^{\textit{M-B, eq}} = 
\Psi^{\textit{M-B, eq}}(\psi, \boldsymbol{\xi}, \boldsymbol{u}) =
\dfrac{\psi}{(2 \pi c_s^2)^{D/2}} 
exp \left[
-\frac{(\boldsymbol{\xi}-\boldsymbol{u})^2}{2 c_s^2}
\right] 
$$

Where $ \psi $  is the quantity of interest (like fluid density or enthalpy), $c_s^2$ is the lattice speed of sound (aka variance of the distribution) and $ D $ is the number of dimensions.
The continuous definition of the central moments is:

$$
\tilde{\kappa}_{mn} = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty}
(\xi_x - u_x)^m (\xi_y -u_y)^n
\Psi(\psi, \boldsymbol{\xi}, \boldsymbol{u}) 
d \xi_x d \xi_y 
$$


```python
from sympy import Symbol, exp, pi, integrate, oo
from sympy import simplify, Float, preorder_traversal
from sympy.matrices import Matrix, eye, diag
from sympy.interactive.printing import init_printing
from sympy import ccode
import sympy as sp
import numpy as np

# init_printing()
```


```python
ex_D2Q9 = Matrix([0, 1, 0, -1, 0, 1, -1, 1, -1])
ey_D2Q9 = Matrix([0, 0, 1, 0, -1, 1, 1, -1, -1])

# Let us choose the following order of moments
# one can denote the variables as f[0], f[1], f[2], f[3]... 
# or f_00, f_10, f_01, f_20 
# We will use the latter notation.
# observe that f[3]=f_20. It is streamed from direction e[-1,0].

order_of_moments = [ 
    (0, 0), 
    (1, 0),
    (0, 1),
    (2, 0),
    (0, 2),
    (1, 1),
    (2, 1),
    (1, 2),
    (2, 2)]

dzeta_x = Symbol('dzeta_x', real=True)
dzeta_y = Symbol('dzeta_y', real=True)

dzeta2D = Matrix([dzeta_x, dzeta_y])

ux = Symbol('u.x')  # don't set real=True for velocity as it freezes the test suite :/
uy = Symbol('u.y')

u2D = Matrix([ux, uy])

# rho = Symbol(r'\rho', positive=True)
# cs2 = Symbol(r'\sigma', positive=True)

rho = Symbol('rho', positive=True)
cs2 = 1./3.
```


```python
def round_and_simplify(stuff):
    simplified_stuff = simplify(stuff)
    rounded_stuff = simplified_stuff

    for a in preorder_traversal(rounded_stuff):
        if isinstance(a, Float):
            rounded_stuff = rounded_stuff.subs(a, round(a, 10))

    rounded_and_simplified_stuff = simplify(rounded_stuff)
    return rounded_and_simplified_stuff
```

## Task

Fill the body of `get_Maxwellian_DF` function and run the script to calculate (central) moments.



```python
class ContinuousCMTransforms:
    def __init__(self, dzeta, u, rho, cs2):
        """
        :param dzeta: direction (x,y,z)
        :param u: velocity (x,y,z) i.e., mean of the distribution
        :param rho: density (not necessarily m00, for instance in multiphase flows)
        :param cs2: variance of the distribution = (speed of sound)^2,
                    for isothermal LB cs2=1./3;
                    otherwise  cs2 = Symbol('RT', positive=True) 

        """
        self.dzeta = dzeta
        self.u = u
        self.rho = rho
        self.sigma2 = cs2

    def get_Maxwellian_DF(self):
        """
        :return: continuous, local Maxwell-Boltzmann distribution       
        """
        # fill...
        return df

    def get_m(self, mno):
        fun = self.get_Maxwellian_DF()
        for dzeta_i, mno_i in zip(self.dzeta, mno):
            fun *= pow(dzeta_i, mno_i)

        lim = [(dim, -oo, oo) for dim in self.dzeta]
        result = integrate(fun, *lim)
        return round_and_simplify(result)

    def get_cm(self, mno):
        # fill...
        return round_and_simplify(result)

```


```python
# here the zeroth moment is calculated
ccmt = ContinuousCMTransforms(dzeta2D, u2D, rho=rho, cs2=cs2)
row0 = order_of_moments[0]
moment0 = ccmt.get_cm(row0)

moment0
```


```python
# write a line of code to calculate the whole vector of moments

# m_eq = Matrix(# fill...)
m_eq
```


```python
# cm_eq = Matrix(# fill...)
cm_eq
```


```python
# next, print is as 'C' code

def print_code(order_of_moments, lhs,rhs):
    for moment, expr in zip(order_of_moments, rhs):
        mstr = [str(m) for m in moment]
        mstr = ''.join(mstr)
        print(f"double {lhs}_{mstr} = {ccode(expr)};")

print_code(order_of_moments, "cm_eq", cm_eq)
```

## Moments of non-equlibrium Distribution Function

The discrete distribution function are streamed along the lattice links, which are defined by a set of discrete velocities,$\textbf{e}$.
Using the Euleran basis and a D2Q9 space, the discrete velocities read,

$$
\textbf{e} = [\textbf{e}_x, \textbf{e}_y], \\
\textbf{e}_x = [0,1,0,-1,0,1,-1,-1,1]^\top, \\
\textbf{e}_y = [0,0,1,0,-1,1,1,-1,-1]^\top, \\
$$

The discrete, raw and central moments are introduced based on the work of Geier et al. [^5] as,

$$ k_{mn} = \sum_{\alpha}(e_{\alpha x})^m ( e_{\alpha y})^n \Psi_{\alpha} $$

while the central moments are calculated in a moving reference frame i.e., with respect to the fluid velocity:

$$ \tilde{k}_{mn} = \sum_{\alpha} ( e_{\alpha x} - u_x)^m ( e_{\alpha y} - u_y)^n \Psi_{\alpha} $$

where $ \Psi_{\alpha} $ is the distribution function of interest (for example hydrodynamic or enthalpy).

Notice, that the equations can be expressed by matrix transformations [^1][^2][^3][^4].

$$
\boldsymbol{\Upsilon} = \mathbb{M} \boldsymbol{\Psi} \\
\boldsymbol{\tilde{\Upsilon}} = \mathbb{N} \boldsymbol{\Upsilon} = \underbrace{\mathbb{N} \mathbb{M}}_{\mathbb{T}} \boldsymbol{\Psi}
$$


where $\boldsymbol{\Upsilon}$ and $\boldsymbol{\tilde{\Upsilon}}$ denote the raw and central moments, respectively.
From the computational point of view, it is preferred to perform the transformations in two steps as in above (without explicit usage of the $\mathbb{T}$ matrix).

Rows of the transformation matrices are calculated analogously to $k$ and $\tilde{k}$, 
$$
M_{mn} = [ (\textbf{e}_x)^m (\textbf{e}_y)^n ]^\top, \\
T_{mn} = [ (\textbf{e}_x - \mathbb{1} u_x)^m (\textbf{e}_y - \mathbb{1} u_y)^n ]
$$
Then, the matrices are assembled row by row as,

$$
\mathbb{M}  
 = 
 \left[
 M_{00}, 
 M_{10}, 
 M_{01},  
 M_{20},
 M_{02},
 M_{11},
 M_{10},
 M_{01},
 M_{22}
 \right]
  \\
\mathbb{T} = 
 \left[
 T_{00}, 
 T_{10}, 
 T_{01}, 
 T_{20},
 T_{02},
 T_{11},
 T_{10},
 T_{01},
 T_{22}
 \right]
$$

The $\mathbb{N}$ matrix can be found as $\mathbb{N} = \mathbb{T} \mathbb{M}^{-1} $.

Observe that $ \mathbb{M} $ is a fixed matrix while $ \mathbb{N} $ depends on the fluid velocity, $ \textbf{u} $. 

Finally, the set of the central moments can be expressed in vector form as,

$$
\boldsymbol{\tilde{\Upsilon}} = 
[\tilde{k}_{00}, \tilde{k}_{10}, \tilde{k}_{01}, \tilde{k}_{20}, \tilde{k}_{02}, \tilde{k}_{11}, \tilde{k}_{21}, \tilde{k}_{12}, \tilde{k}_{22}]^\top.
$$

The physical interpretation of the raw, zeroth and first order moments of the hydrodynamic  DF corresponds to the values of density, $ \rho $ and momentum $  \rho \textbf{u} $.





```python

class MatrixGenerator:
    def __init__(self, ex, ey, order_of_moments):
        self.ex = ex
        self.ey = ey
        self.order_of_moments = order_of_moments

    def __matrix_maker(self, row_maker_fun):
        M = [row_maker_fun(*row) for row in self.order_of_moments]
        return M

    def get_raw_moments_matrix(self):
        """
        :return: transformation matrix from DF to raw moments
        """
    
        def get_row(m, n):
            row = [pow((self.ex[i]), m) * pow((self.ey[i]), n)  for i in range(0, 9)]
            return row

        m_ = self.__matrix_maker(get_row)
        # M = [get_row(*row) for row in self.order_of_moments] # same as
        return Matrix(m_)

    def get_T_matrix(self):
        """
        :return: transformation matrix from DF to central moments
        """

        def get_row(m, n):
            # fill...
            row = [round_and_simplify(r) for r in row] # simplify the elements in each row
            return row

        m_ = self.__matrix_maker(get_row)
        return Matrix(m_)

```


```python
matrixGenerator = MatrixGenerator(ex_D2Q9, ey_D2Q9, order_of_moments)

Mraw = matrixGenerator.get_raw_moments_matrix()
Mraw
```


```python
Traw = matrixGenerator.get_T_matrix()
Nraw = Traw * Mraw.inv()

Nraw = Matrix([round_and_simplify(Nraw[i,:]) for i in range(9)])
Nraw
```

## Task
We have just generate the matrix of transformation. 
Now, let as create the vector of variables which are going to be transformed.
Implement the `get_symbols` function. It shall return a vector (1-D Matrix, i.e. `Matrix([stuff])` ) having the following form $ [f_{00}, f_{10}, f_{01}, f_{20}, f_{02}, etc...] $


```python

def get_symbols(name, directions):
    print_symbols = []
    # fill...
    return Matrix(print_symbols)

fs = get_symbols("f", order_of_moments)
fs
```


```python
m = Mraw * fs
m
```


```python
print("//raw moments from density-probability functions")
print_code(order_of_moments, "m", m)
```


```python
ms = get_symbols("m", order_of_moments)
cm = Nraw * ms
cm
```


```python
print("//central moments from raw moments")
print_code(order_of_moments, "cm", cm)
```


```python
# RELAXATION MATRIX
omega_v = Symbol('omega_nu', positive=True)
omega_b = Symbol('omega_bulk', positive=True) 

s_plus_D2Q9 = (omega_b + omega_v) / 2
s_minus_D2Q9 = (omega_b - omega_v) / 2

S_relax_hydro_D2Q9 = diag(1, 1, 1, s_plus_D2Q9, s_plus_D2Q9, omega_v, 1, 1, 1)
S_relax_hydro_D2Q9[3, 4] = s_minus_D2Q9
S_relax_hydro_D2Q9[4, 3] = s_minus_D2Q9
```


```python
cm_after_collision = eye(9) * cm + S_relax_hydro_D2Q9 * (cm_eq - cm)
print("//collision in central moments space")
print_code(order_of_moments, "cm_after_collision", cm_after_collision)
```

## Summary

That's the magic - you have learned how perform symbolic computations and generate code from it.
The back-tranformation from central moments to moments, then from moments to distribution function follow the same way.


References:

[^1]: Linlin Fei, Kai Hong Luo, 'Cascaded lattice Boltzmann method for incompressible thermal flows with heat sources and general thermal boundary conditions' Computers and Fluids (2018).

[^2]: Linlin Fei, Kai Hong Luo, Chuandong Lin, Qing Li, 'Modeling incompressible thermal flows using a central-moments-based lattice Boltzmann method' International Journal of Heat and Mass Transfer (2017).

[^3]: Linlin Fei and Kai Hong Luo, 'Consistent forcing scheme in the cascaded lattice Boltzmann method' Physical Review E 96, 053307 (2017).

[^4]: Linlin Fei, Kai H. Luo and Qing Li, 'Three-dimensional cascaded lattice Boltzmann method: Improved implementation and consistent forcing scheme' Physical Review E 97, 053309 (2018)

[^5]: M. Geier, A. Greiner, J. G. Korvink, 'Cascaded digital lattice Boltzmann automata for high Reynolds number flow' Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 73 (2006).

[^6]: Xiaoyi He, Xiaowen Shan, and Gary D. Doolen, 'Discrete Boltzmann equation model for nonideal gases' in Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics (1998).


aec6dd4f42b510b153dbd200e6d067b76b654925
ipynb
Jupyter Notebook
lab4_lbm/lbm_theory_to_fill.ipynb
ggruszczynski/CFDPython
1662ede061fb899d6ed3f89c17877e65f65e521c
CC-BY-3.0
lab4_lbm/lbm_theory_to_fill.ipynb
ggruszczynski/CFDPython
1662ede061fb899d6ed3f89c17877e65f65e521c
CC-BY-3.0
lab4_lbm/lbm_theory_to_fill.ipynb
ggruszczynski/CFDPython
1662ede061fb899d6ed3f89c17877e65f65e521c
CC-BY-3.0
2021-02-05T08:00:02.000Z
2021-02-05T08:00:02.000Z
Qwen/Qwen-72B
1. YES
2. YES