from lbmpy.session import *

from pystencils.sympyextensions import prod
from lbmpy.moments import get_default_moment_set_for_stencil, moments_up_to_order
from lbmpy.creationfunctions import create_lb_method
from lbmpy.equilibrium import discrete_equilibrium_from_matching_moments
from lbmpy.methods import create_from_equilibrium, DensityVelocityComputation
from lbmpy.quadratic_equilibrium_construction import *
from lbmpy.chapman_enskog import ChapmanEnskogAnalysis
from lbmpy.moments import exponent_to_polynomial_representation

Demo: Theoretical Background - LB Equilibrium Construction using quadratic Ansatz

According to book by Wolf-Gladrow “Lattice-Gas Cellular Automata and Lattice Boltzmann Methods” (2005)

Through the Chapman Enskog analysis the following necessary conditions can be found in order for a lattice Boltzmann Method to approximate the Navier Stokes mass and momentum conservation equations. In the Chapman Enskog analysis only the moments of the equilibrium distribution functions are used, thus all conditions are formulated with regard to the moments \(\Pi\) of the equilibrium distribution function \(f^{(eq)}\)

The conditions are: - zeroth moment is the density: \(\Pi_0 = \sum_q f^{(eq)}_q = \rho\) - first moment is the momentum density, or for incompressible models the velocity: - compressible: \(\Pi_\alpha = \sum_q c_{q\alpha} f^{(eq)}_q = \rho u_\alpha\) - incompressible: \(\Pi_\alpha = \sum_q c_{q\alpha} f^{(eq)}_q = u_\alpha\) - second moment is related to the pressure tensor and has to be: \(\Pi_{\alpha\beta} = \sum_q c_{q\alpha} c_{q\beta} f^{(eq)}_q = \rho u_\alpha u_\beta + p \delta_{\alpha\beta}\) - third order moments are also used in the Chapman Enskog expansion. The conditions on these moments are harder to formulate and are investigated later. A commonly used, but overly restrictive choice is \(\Pi_{\alpha\beta\gamma} = p ( \delta_{\alpha\beta} u_\gamma + \delta_{\alpha\gamma} u_\beta + \delta_{\beta\gamma} u_\alpha )\). In Wolf-Gladrows book these conditions on the third order moment are not used but implicitly fulfilled by choosing fixed fractions of the coefficients \(\frac{A_1}{A_2}\) etc.

Now the following generic quadratic ansatz is used for the equilibrium distribution.

\[f^{(eq)}_q = A_{|q|} + B_{|q|} (\mathbf{c}_q \cdot \mathbf{u}^2 ) + D_{|q|} \mathbf{u}^2\]

The free parameters \(A_{|q|}, B_{|q|}, C_{|q|}\) and \(D_{|q|}\) are chosen such that above conditions are fulfilled. The subscript \(|q|\) is an integer and defined as the sum of the absolute values of the corresponding stencil direction. For example: for center \(|q|=0\), for direct neighbors like north, east, top \(|q|=1\) for 2D diagonals like north-west its 2 and for 3D diagnoals like bottom-north-west its 3.

lbmpy can create this quadratic ansatz for use for a given stencil:

# Create Stencil
stencil = LBStencil(Stencil.D2Q9)

# Create quadratic equilibrium ansatz
ansatz = generic_equilibrium_ansatz(stencil)

# Show equilibrium for each stencil direction
plt.figure(figsize=(12,8))
stencil.plot(data=ansatz, textsize=9, slice=True)

Next we define the restrictions obtained through the Chapman Enskog analysis, in the book listed as equations (5.4.2) and following:

moment_restrictions = hydrodynamic_moment_values(dim=len(stencil[0]), compressible=True, up_to_order=2)
moment_restrictions

$\displaystyle \left\{ \left( 0, \ 0\right) : \rho, \ \left( 0, \ 1\right) : \rho u_{1}, \ \left( 0, \ 2\right) : p + \rho u_{1}^{2}, \ \left( 1, \ 0\right) : \rho u_{0}, \ \left( 1, \ 1\right) : \rho u_{0} u_{1}, \ \left( 2, \ 0\right) : p + \rho u_{0}^{2}\right\}$

The parameter up_to_order can be modified to 3. Then the third order restrictions are included as well (see discussion above). Using these moment restrictions, the necessary conditions on the parameter \(A\) to \(D\) can be found.

equations = moment_constraint_equations(stencil, ansatz, moment_restrictions)
sp.Matrix(equations)

$\displaystyle \left[\begin{matrix}8 C_{2} - \rho\\2 B_{1} + 4 B_{2} - \rho\\4 C_{2} + 2 D_{1} + 4 D_{2}\\2 C_{1} + 4 C_{2} + D_{0} + 4 D_{1} + 4 D_{2}\\A_{0} + 4 A_{1} + 4 A_{2} - \rho\\2 A_{1} + 4 A_{2} - p\\2 C_{1} + 4 C_{2} + 2 D_{1} + 4 D_{2} - \rho\end{matrix}\right]$

Since we have still more unknowns than equations, some additional restrictions have to be imposed.

dofs = generic_equilibrium_ansatz_parameters(stencil)
dofs

$\displaystyle \left[ A_{0}, \ A_{1}, \ A_{2}, \ B_{0}, \ B_{1}, \ B_{2}, \ C_{0}, \ C_{1}, \ C_{2}, \ D_{0}, \ D_{1}, \ D_{2}, \ p\right]$

In Wolf-Gladrows book the following arbitrary restrictions are added to the necessary constraints:

\[\frac{A_0}{A_1} = \frac{A_1}{A_2} = \frac{B_1}{B_2} = \frac{D_0}{D_1} =: r\]

additional_restrictions = [
    "A_0 / A_1 - r",
    "A_1 / A_2 - r",
    "B_1 / B_2 - r",
    "D_0 / D_1 - r",
    "D_1 / D_2 - r", # comment out this line to get solution dependent on r
]
additional_restrictions = [sp.sympify(e) for e in additional_restrictions]
additional_restrictions

$\displaystyle \left[ \frac{A_{0}}{A_{1}} - r, \ \frac{A_{1}}{A_{2}} - r, \ \frac{B_{1}}{B_{2}} - r, \ \frac{D_{0}}{D_{1}} - r, \ \frac{D_{1}}{D_{2}} - r\right]$

solveResult = sp.solve(equations + additional_restrictions, dofs + [sp.Symbol("r")], dict=True)
solveResult

$\displaystyle \left[ \left\{ A_{0} : \frac{4 \rho}{9}, \ A_{1} : \frac{\rho}{9}, \ A_{2} : \frac{\rho}{36}, \ B_{1} : \frac{\rho}{3}, \ B_{2} : \frac{\rho}{12}, \ C_{1} : \frac{\rho}{2}, \ C_{2} : \frac{\rho}{8}, \ D_{0} : - \frac{2 \rho}{3}, \ D_{1} : - \frac{\rho}{6}, \ D_{2} : - \frac{\rho}{24}, \ p : \frac{\rho}{3}, \ r : 4\right\}\right]$

The equilibrium we found here is the same as obtained with the usual lbmpy construction technique:

constructed_equilibrium = sp.Matrix(ansatz).subs(solveResult[0]).expand()
lbm_config = LBMConfig(stencil=stencil, compressible=True, zero_centered=False)
normal_lbmpy_equilibrium = create_lb_method(lbm_config=lbm_config).get_equilibrium_terms()
assert constructed_equilibrium == normal_lbmpy_equilibrium

Generalization of above technique

def generic_polynomial(u, coeff_prefix, order=2):
    dim = len(u)
    result = 0
    for idx in moments_up_to_order(order, dim=dim):
        u_prod = prod(u[i] ** exp for i, exp in enumerate(idx))
        coeff = sp.Symbol(("%s_" + ("%d" * dim)) % ((coeff_prefix,) + idx))
        result += coeff * u_prod
    return result


def generic_polynomial_coeffs(dim, coeff_prefix, order=2):
    result = []
    for idx in moments_up_to_order(order, dim=dim):
        result.append(sp.Symbol(("%s_" + ("%d" * dim)) % ((coeff_prefix,) + idx)))
    return result

allMoments = get_default_moment_set_for_stencil(stencil)
allMoments

$\displaystyle \left[ 1, \ x, \ y, \ x^{2}, \ y^{2}, \ x y, \ x^{2} y, \ x y^{2}, \ x^{2} y^{2}\right]$

We use, as before, all constraints for moments up order 2. This time we do not use a quadratic ansatz for the equilibrium distribution or the additional constraints (i.e. the ratios of coefficients being some constant \(r\)). Instead an arbitrary second order polynomial \(u\) is used for the third order moments. The forth order moment does not appear at all in the Chapman Enskog expansion and is thus set to a constant.

So we end up with the following moments

dim = len(stencil[0])
u = sp.symbols("u_:3")[:len(stencil[0])]
moment_restrictions = hydrodynamic_moment_values(dim=len(stencil[0]), compressible=True, up_to_order=2)
moment_restrictions[(2, 1)] = generic_polynomial(u, "a")
moment_restrictions[(1, 2)] = generic_polynomial(u, "b")
moment_restrictions[(2, 2)] = sp.Symbol("M_22")
moment_restrictions = {m: v.subs(sp.Symbol("p"), sp.Symbol("rho") / 3) for m, v in moment_restrictions.items()}
moment_restrictions

$\displaystyle \left\{ \left( 0, \ 0\right) : \rho, \ \left( 0, \ 1\right) : \rho u_{1}, \ \left( 0, \ 2\right) : \rho u_{1}^{2} + \frac{\rho}{3}, \ \left( 1, \ 0\right) : \rho u_{0}, \ \left( 1, \ 1\right) : \rho u_{0} u_{1}, \ \left( 1, \ 2\right) : b_{00} + b_{01} u_{1} + b_{02} u_{1}^{2} + b_{10} u_{0} + b_{11} u_{0} u_{1} + b_{20} u_{0}^{2}, \ \left( 2, \ 0\right) : \rho u_{0}^{2} + \frac{\rho}{3}, \ \left( 2, \ 1\right) : a_{00} + a_{01} u_{1} + a_{02} u_{1}^{2} + a_{10} u_{0} + a_{11} u_{0} u_{1} + a_{20} u_{0}^{2}, \ \left( 2, \ 2\right) : M_{22}\right\}$

Next all parameters are collected..

parameters = generic_polynomial_coeffs(dim, "a") + generic_polynomial_coeffs(dim, "b") + [sp.Symbol("p")]
parameters

$\displaystyle \left[ a_{00}, \ a_{01}, \ a_{10}, \ a_{02}, \ a_{11}, \ a_{20}, \ b_{00}, \ b_{01}, \ b_{10}, \ b_{02}, \ b_{11}, \ b_{20}, \ p\right]$

… and a lbmpy LB method is created. On this method an automatic Chapman Enskog analysis can be conducted to find constraints for the free parameters above.

rr = sp.symbols("omega")
eq_values = {exponent_to_polynomial_representation(m) : v for m, v in moment_restrictions.items()}
r_rates = {exponent_to_polynomial_representation(m) : rr for m in moment_restrictions.keys()}
eq_values, r_rates

$\displaystyle \left( \left\{ 1 : \rho, \ x : \rho u_{0}, \ x^{2} : \rho u_{0}^{2} + \frac{\rho}{3}, \ y : \rho u_{1}, \ y^{2} : \rho u_{1}^{2} + \frac{\rho}{3}, \ x y : \rho u_{0} u_{1}, \ x y^{2} : b_{00} + b_{01} u_{1} + b_{02} u_{1}^{2} + b_{10} u_{0} + b_{11} u_{0} u_{1} + b_{20} u_{0}^{2}, \ x^{2} y : a_{00} + a_{01} u_{1} + a_{02} u_{1}^{2} + a_{10} u_{0} + a_{11} u_{0} u_{1} + a_{20} u_{0}^{2}, \ x^{2} y^{2} : M_{22}\right\}, \ \left\{ 1 : \omega, \ x : \omega, \ x^{2} : \omega, \ y : \omega, \ y^{2} : \omega, \ x y : \omega, \ x y^{2} : \omega, \ x^{2} y : \omega, \ x^{2} y^{2} : \omega\right\}\right)$

cqe = DensityVelocityComputation(stencil, True, False)
equilibrium = discrete_equilibrium_from_matching_moments(stencil=stencil, moment_constraints=eq_values,
                                                         zeroth_order_moment_symbol=cqe.density_symbol,
                                                         first_order_moment_symbols=cqe.velocity_symbols)
method = create_from_equilibrium(stencil, equilibrium, cqe, r_rates)
analysis = ChapmanEnskogAnalysis(method, constants=set(parameters))

analysis.does_approximate_navier_stokes()

$\displaystyle \left\{0, - \frac{b_{10}}{2} + \frac{b_{10}}{\omega}, \frac{a_{01}}{2} - \frac{a_{01}}{\omega} - \frac{b_{10}}{2} + \frac{b_{10}}{\omega}, \frac{a_{01}}{2} - \frac{a_{01}}{\omega} - \frac{a_{10}}{2} + \frac{a_{10}}{\omega} - \frac{\rho}{6} + \frac{\rho}{3 \omega}, \frac{a_{01}}{2} - \frac{a_{01}}{\omega} - \frac{b_{01}}{2} + \frac{b_{01}}{\omega} - \frac{\rho}{6} + \frac{\rho}{3 \omega}, \frac{a_{01}}{2} - \frac{a_{01}}{\omega} + b_{10} - \frac{2 b_{10}}{\omega} - \frac{\rho}{2} + \frac{\rho}{\omega}\right\}$

These constraints can be solved for the free parameters:

solveRes = sp.solve(analysis.does_approximate_navier_stokes(), parameters)
solveRes

$\displaystyle \left[ \right]$

new_moment_restrictions = {a : b.subs(solveRes) for a, b in moment_restrictions.items()}
new_moment_restrictions

$\displaystyle \left\{ \left( 0, \ 0\right) : \rho, \ \left( 0, \ 1\right) : \rho u_{1}, \ \left( 0, \ 2\right) : \rho u_{1}^{2} + \frac{\rho}{3}, \ \left( 1, \ 0\right) : \rho u_{0}, \ \left( 1, \ 1\right) : \rho u_{0} u_{1}, \ \left( 1, \ 2\right) : b_{00} + b_{01} u_{1} + b_{02} u_{1}^{2} + b_{10} u_{0} + b_{11} u_{0} u_{1} + b_{20} u_{0}^{2}, \ \left( 2, \ 0\right) : \rho u_{0}^{2} + \frac{\rho}{3}, \ \left( 2, \ 1\right) : a_{00} + a_{01} u_{1} + a_{02} u_{1}^{2} + a_{10} u_{0} + a_{11} u_{0} u_{1} + a_{20} u_{0}^{2}, \ \left( 2, \ 2\right) : M_{22}\right\}$

All methods constructed with these constraints should theoretically approximate Navier Stokes.