Symbolic Operations and CLBM

Before start, it is assumed that the reader is already familiar with:

Basic model creation in the TCLB code, see first tutorial

Basics of LB physics, see SRT LBM and DDF Heat Transfer

We will cover following topics:

Symbolic computations (within) TCLB.
Heat transfer on D2Q9 lattice using Central Moments (CM) a.k.a. Cascaded LBM. This is an enhancement to the previous tutorial, where the SRT collision was performed.

Notion of moments

The raw and central moments are introduced based on the work of Geier et al. ⁵ 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 (either hydrodynamic or enthalpy).

Alternatively, formerly equations can be expressed by matrix transformations ¹²³⁴.

Note

Raw moments are obtained from distribution functions through $\boldsymbol{\Upsilon} = \mathbb{M} \boldsymbol{g}$ . Central moments are derived from raw moments by $\boldsymbol{\tilde{\Upsilon}}= \mathbb{N} \boldsymbol{\Upsilon}$ . It is important to notice that $\mathbb{M}$ is a fixed matrix while $\mathbb{N}$ depends on the fluid velocity, $\textbf{u}$ .

The form of the matrices can be found here. The resulting order of the central moments is,

$\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 order moments of the hydrodynamic and enthalpy DF corresponds to the values of density, $\rho$ , and enthalpy $H=\rho c_p T$ . For simplicity we take $cp=1=const$ in the remainder of the tutorial.

$\begin{eqnarray} \rho &=& k_{00} = \sum_\alpha f_\alpha, \\ \rho c_p T &=& k_{00}^H = \sum_\alpha h_\alpha, \end{eqnarray}$

The macroscopic fluid momentum is the first raw moment of the hydrodynamic DF and has a dependency on the forcing term $\mathbf{F}$ .

$\begin{eqnarray} \rho \textbf{u} = \begin{bmatrix} u_{x} \\ u_{y} \\ \end{bmatrix} = \begin{bmatrix} k_{10} \\ k_{01} \\ \end{bmatrix} = \sum_\alpha f_\alpha \textbf{e}_\alpha + \frac{\textbf{F}}{2 \rho} \delta t. \end{eqnarray}$

Collision in the Central Moment Space

The collision process and application of forcing terms is conducted in central moment space to improve computational efficiency ²³,

$\begin{eqnarray} % hydrodynamics \boldsymbol{\tilde{\Upsilon}}^{*}(\textbf{x}, t)& = & (\mathbb{1} - \mathbb{S})\boldsymbol{\tilde{\Upsilon}} + \mathbb{S} \boldsymbol{\tilde{\Upsilon}}^{eq} + (\mathbb{1} - \mathbb{S}/2)\tilde{\textbf{F}}, \\ % heat \boldsymbol{\tilde{\Upsilon}}^{H,*}(\textbf{x}, t)& = & (\mathbb{1} - \mathbb{S}^{H}) \boldsymbol{\tilde{\Upsilon}}^H + \mathbb{S}^{H} \boldsymbol{\tilde{\Upsilon}}^{H, eq} \label{eq:pf_collision_streaming}. \end{eqnarray}$

The relaxation matrices, $\mathbb{S}$ and $\mathbb{S}^{H}$ , have a diagonal form specified as,

$\begin{eqnarray} \mathbb{S}& = &diag\left( [s_0, s_1, s_1], \begin{bmatrix} s_{+}, & s_{-} \\ s_{-}, & s_{+} \\ \end{bmatrix}, [s_\nu , s_3, s_3, s_4] \right), \label{eq:S_relaxation_matrix} \\ \mathbb{S}^{H}& = &diag\left( [s_0^{H}, s_k^{H}, s_k^{H},s_2^{H}, s_2^{H}, s_2^{H} , s_3^{H}, s_3^{H}, s_4^{H}] \right), \label{eq:S_phi_relaxation_matrix} \end{eqnarray}$

where $s_{+} = (s_b + s_v)/2$ and $s_{-} = (s_b - s_v)/2$ . Following the Chapman-Enskog expansion, the kinematic and bulk viscosities are used to calculate corresponding relaxation frequencies $s_v$ and $s_b$ ,

$\begin{eqnarray} s_v = \frac{1}{\tau} = \frac{1}{\frac{\nu}{c_s^2 \delta t } + 1/2}, \\ s_b = \frac{1}{\tau_b} = \frac{1}{\frac{\zeta}{c_s^2 \delta t } + 1/2}. \end{eqnarray}$

The bulk relaxation rate is specified as, $s_b =1$ . Such an approach results in a reduction of non-physical pressure oscillations, giving a bulk equilibrium after each collision and a more stable model. To satisfy conservation laws for density and momentum, the zeroth and first order relaxation constants are set to unity, $s_0 = s_1 = 1$ . For this study, the tunable relaxation frequencies are specified as, $s_3 = s_4 =1$ .

In the case of the Enthalpy DF, the relaxation frequency, $s_v$ , is replaced by,

$\begin{eqnarray} s_k^{H} = \frac{1}{\tau_{k}} = \frac{1}{\frac{k}{c_s^2 \delta t } + 1/2}. \end{eqnarray}$

Here, $s_k^{H}$ is used to relax the first moments only (see ²) while the remainder are conserved by setting the relaxation frequencies to unity, $s_0^{H} = s_2^{H} = s_3^{H} =s_4^{H} = 1$ .

streaming

To perform the streaming step, the relaxed, post-collision central moments need to be transformed back into the velocity distribution space and streamed to neighbours, $\begin{eqnarray} % hydrodynamics f_{\alpha}(\textbf{x} + \textbf{e}\delta t, t + \delta t ) & = & \mathbb{M}^{-1} \mathbb{N}^{-1} \boldsymbol{\tilde{\Upsilon}}^{*}_{\alpha}(\textbf{x}, t), \\ % heat h_{\alpha}(\textbf{x} + \textbf{e}\delta t, t + \delta t ) & = & \mathbb{M}^{-1} \mathbb{N}^{-1} \boldsymbol{\tilde{\Upsilon}}^{H,*}_{\alpha}(\textbf{x}, t). \end{eqnarray}$

Model Creation in TCLB

Dynamics.R

Take the one from the basic D2Q9 Heat Transfer tutorial and add a new node type AddNodeType("CM","COLLISION"), since we are going to implement a new collision kernel.

Dynamics.c

First we have to load R symbols.

<?R
    source("conf.R") 
    source("lib/boundary.R")
    c_header();

# Creating variables for symbolic computations
    f = PV(DensityAll$name[DensityAll$group=="f"])
    h = PV(DensityAll$name[DensityAll$group=="h"])
    rho =  PV("rho")
    u = PV(c("ux","uy"))
    rhoT = PV("rhoT")

# Extracting velocity set
    U = d2q9

?>

Macroscopic Quantities

Now we can replace the tedious and error-prone operarions from the previous tutorial with their symbolic equivalents. Using R language:

CudaDeviceFunction real_t getRho(){
    return <?R C(sum(f)) ?>; // compiles to f[8]+f[7]+f[6]+f[5]+f[4]+f[3]+f[2]+f[1]+f[0];
}

CudaDeviceFunction real_t getT(){
    return (<?R C(sum(h)) ?>)/(<?R C(sum(f)) ?>);
}

CudaDeviceFunction vector_t getRawU(){
    real_t d = getRho();
    vector_t u;
<?R C(PV(c("u.x","u.y", "u.z")), f %*% U) ?>
    u.x /= d;
    u.y /= d;
    u.z = 0;
    return u;
}

CudaDeviceFunction vector_t getU()
{
    real_t localTemperature = getT();
    vector_t u = getRawU();
    real_t m00 = getRho();
    vector_t Force = getForce(localTemperature, m00);
    u.x += Force.x/(2*m00);
    u.y += Force.y/(2*m00);
    u.z = 0;
    return u;
}

Collision

CudaDeviceFunction void CollisionCM()
{
    real_t localTemperature = getT();
    real_t m00 = getRho();

    vector_t Force = getForce(localTemperature, m00);
    vector_t u = getRawU();
    u.x += Force.x/(2*m00);
    u.y += Force.y/(2*m00);

    relax_and_collide_hydro_with_F(f, omega_nu, u, Force);
    relax_and_collide_ADE(h, omega_k, u);
}

We are going to implement following functions in Central Moments Space:

CudaDeviceFunction void SetEquilibrium(real_t x_in[9], real_t rho_Xeq, vector_t u){...}
CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t x_in[9], real_t omega_nu, vector_t u, vector_t F){...}
CudaDeviceFunction void relax_and_collide_ADE(real_t x_in[9], real_t omega_ade, vector_t u) {...}

Symbolic Operations with python

Theory

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

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

$\begin{eqnarray} \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] \label{eq:Maxwellian} \end{eqnarray}$

Where $\psi$ is the quantity of interest (like density or enthalpy) and $D$ is the number of dimensions. The continuous definition of the central moments is:

$\begin{eqnarray} \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 \label{eq:cont_cm_mom_def} \end{eqnarray}$

Now, we will discuss transformation of forcing term.

To compute the continuous transform, a scheme proposed by He et al. ⁶ was expressed in a continous form by Premnath al. ⁷.

$\begin{eqnarray} \Delta \boldsymbol{F}_{He} &=& \dfrac{\boldsymbol{F}}{\rho} \cdot \nabla_u \Psi^{\textit{M-B, eq}}(k_{00}, \boldsymbol{\xi}, \boldsymbol{u}) \nonumber \\ \, &=& \dfrac{\boldsymbol{F}}{\rho} \cdot \frac{{(\boldsymbol{\xi}}-\boldsymbol{u})}{c_s^2} \Psi^{\textit{M-B, eq}}(k_{00}, \boldsymbol{\xi}, \boldsymbol{u}) \label{eq:He_force_cont} \end{eqnarray}$

Practise

Although it is possible to inject python snippets using <?python /> within TCLB code, we will run the symbolic calculations in PyCharm enviroment then copy-paste the output.

Clone the repository containing TCLB tools:

git clone https://github.com/CFD-GO/TCLB_tools.git

Then open the TCLB_tools/Python/symbolic_tools/ directory in PyCharm and create a new pycharm-project. Next, go tutorials folder and open the tutorial_cm_eq_derivation.py script. Run it.

Question

Why the python's timer doesn't work as expected?

It turns out the obtained formulas are simple, thus the numerical effort of Cascaded collision is similar to SRT (~ +10%). As it takes a while to compute the integrals, the result is hardcoded and used afterwards for symbolic code generation.

Collision Kernel

Now, it is time to generate the missing functions:

CudaDeviceFunction void SetEquilibrium(real_t x_in[9], real_t rho_Xeq, vector_t u){...}
CudaDeviceFunction void relax_and_collide_hydro_with_F(real_t x_in[9], real_t omega_nu, vector_t u, vector_t F){...}
CudaDeviceFunction void relax_and_collide_ADE(real_t x_in[9], real_t omega_ade, vector_t u) {...}

To do so, open a prepared tutorial_cm_collision.py script. Fill the missing functionalities #TODO

Hint: After the collision, the central moments have to be back-tranformed to moments, then to density-probability functions - see the streaming step.

Setting up a Simulation

To run the new collision kernel, change the configuration file my_heat_cube.xml from previous tutorial namely,

from <BGK><Box /></BGK> to <CM><Box /></CM>

and uncomment the corresponding switch //case NODE_CM: ... in Dynamics.c.Rt.

References

Tip

A good starting point explaining the 'zoology' of various LBM models is a modern (2017) book 'The Lattice Boltzmann Method: Principles and Practice' written by T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E.M. Viggen .

The implementation of the cascaded collision kernel is greatly insipered by the work of Linlin Fei, Kai H. Luo et al.

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). ↩
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). ↩↩↩
Linlin Fei and Kai Hong Luo, 'Consistent forcing scheme in the cascaded lattice Boltzmann method' Physical Review E 96, 053307 (2017). ↩↩
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) ↩
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). ↩
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). ↩
K. N. Premnath, S. Banerjee, 'Incorporating forcing terms in cascaded lattice Boltzmann approach by method of central moments' in Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 80 (2009). ↩