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
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:
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
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)
Since we have still more unknowns than equations, some additional restrictions have to be imposed.
dofs = generic_equilibrium_ansatz_parameters(stencil)
dofs
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
solveResult = sp.solve(equations + additional_restrictions, dofs + [sp.Symbol("r")], dict=True)
solveResult
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
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
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
Next all parameters are collected..
parameters = generic_polynomial_coeffs(dim, "a") + generic_polynomial_coeffs(dim, "b") + [sp.Symbol("p")]
parameters
... 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
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()
These constraints can be solved for the free parameters:
solveRes = sp.solve(analysis.does_approximate_navier_stokes(), parameters)
solveRes
new_moment_restrictions = {a : b.subs(solveRes) for a, b in moment_restrictions.items()}
new_moment_restrictions
All methods constructed with these constraints should theoretically approximate Navier Stokes.