Demo: Create lbmpy Method from Scratch¶
Defining transformation to collision space¶
[2]:
from lbmpy.moments import moment_matrix, moments_up_to_component_order, exponents_to_polynomial_representations
moment_exponents = list(moments_up_to_component_order(2, 2))
moment_exponents
[2]:
$\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]$
[3]:
moments = exponents_to_polynomial_representations(moment_exponents)
moments
[3]:
$\displaystyle \left( 1, \ y, \ y^{2}, \ x, \ x y, \ x y^{2}, \ x^{2}, \ x^{2} y, \ x^{2} y^{2}\right)$
[4]:
d2q9 = LBStencil(Stencil.D2Q9, ordering='walberla')
d2q9.plot()

[5]:
M = moment_matrix(moments, stencil=d2q9)
M
[5]:
$\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]$
Defining the Equilibrium Distribution, and Computation of Conserved Quantities¶
[6]:
from lbmpy.equilibrium import ContinuousHydrodynamicMaxwellian
from lbmpy.methods import DensityVelocityComputation
equilibrium = ContinuousHydrodynamicMaxwellian(dim=2, compressible=True, order=2)
cqc = DensityVelocityComputation(d2q9, True, False)
tuple(zip(moments, equilibrium.moments(moments)))
[6]:
$\displaystyle \left( \left( 1, \ \rho\right), \ \left( y, \ \rho u_{1}\right), \ \left( y^{2}, \ c_{s}^{2} \rho + \rho u_{1}^{2}\right), \ \left( x, \ \rho u_{0}\right), \ \left( x y, \ \rho u_{0} u_{1}\right), \ \left( x y^{2}, \ c_{s}^{2} \rho u_{0}\right), \ \left( x^{2}, \ c_{s}^{2} \rho + \rho u_{0}^{2}\right), \ \left( x^{2} y, \ c_{s}^{2} \rho u_{1}\right), \ \left( x^{2} y^{2}, \ c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2}\right)\right)$
Defining Relaxation Behaviour and Creating the Method¶
[7]:
from lbmpy.methods.creationfunctions import create_from_equilibrium
omega = sp.symbols("omega")
relaxation_rate_dict = {moment : omega for moment in moments}
force_model = forcemodels.Guo(sp.symbols("F_:2"))
method = create_from_equilibrium(d2q9, equilibrium, cqc, relaxation_rate_dict)
method
[7]:
Moment-Based Method | Stencil: D2Q9 | Zero-Centered Storage: ✗ | Force Model: None |
---|
Continuous Hydrodynamic Maxwellian Equilibrium | $f (\rho, \left( u_{0}, \ u_{1}\right), \left( v_{0}, \ v_{1}\right)) = \frac{\rho e^{\frac{- \left(- u_{0} + v_{0}\right)^{2} - \left(- u_{1} + v_{1}\right)^{2}}{2 c_{s}^{2}}}}{2 \pi c_{s}^{2}}$ | ||
---|---|---|---|
Compressible: ✓ | Deviation Only: ✗ | Order: 2 |
Relaxation Info | ||
---|---|---|
Moment | Eq. Value | Relaxation Rate |
$1$ | $\rho$ | $\omega$ |
$y$ | $\rho u_{1}$ | $\omega$ |
$y^{2}$ | $c_{s}^{2} \rho + \rho u_{1}^{2}$ | $\omega$ |
$x$ | $\rho u_{0}$ | $\omega$ |
$x y$ | $\rho u_{0} u_{1}$ | $\omega$ |
$x y^{2}$ | $c_{s}^{2} \rho u_{0}$ | $\omega$ |
$x^{2}$ | $c_{s}^{2} \rho + \rho u_{0}^{2}$ | $\omega$ |
$x^{2} y$ | $c_{s}^{2} \rho u_{1}$ | $\omega$ |
$x^{2} y^{2}$ | $c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2}$ | $\omega$ |
Example of a update equation without simplifications¶
[8]:
collision_rule = method.get_collision_rule()
collision_rule
[8]:
Subexpressions:
$$vel0Term \leftarrow f_{4} + f_{6} + f_{8}$$ |
$$vel1Term \leftarrow f_{1} + f_{5}$$ |
$$\rho \leftarrow f_{0} + f_{2} + f_{3} + f_{7} + vel0Term + vel1Term$$ |
$$u_{0} \leftarrow \frac{- f_{3} - f_{5} - f_{7} + vel0Term}{\rho}$$ |
$$u_{1} \leftarrow \frac{- f_{2} + f_{6} - f_{7} - f_{8} + vel1Term}{\rho}$$ |
$$\delta_{\rho} \leftarrow \rho - 1$$ |
Main Assignments:
$$d_{0} \leftarrow f_{0} + \omega \left(c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2} - f_{5} - f_{6} - f_{7} - f_{8}\right) - \omega \left(c_{s}^{2} \rho - f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{1}^{2}\right) - \omega \left(c_{s}^{2} \rho - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{0}^{2}\right) + \omega \left(- f_{0} - f_{1} - f_{2} - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \rho\right)$$ |
$$d_{1} \leftarrow f_{1} - \frac{\omega \left(c_{s}^{2} \rho u_{1} - f_{5} - f_{6} + f_{7} + f_{8}\right)}{2} + \frac{\omega \left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \rho u_{1}\right)}{2} - \frac{\omega \left(c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2} - f_{5} - f_{6} - f_{7} - f_{8}\right)}{2} + \frac{\omega \left(c_{s}^{2} \rho - f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{1}^{2}\right)}{2}$$ |
$$d_{2} \leftarrow f_{2} + \frac{\omega \left(c_{s}^{2} \rho u_{1} - f_{5} - f_{6} + f_{7} + f_{8}\right)}{2} - \frac{\omega \left(- f_{1} + f_{2} - f_{5} - f_{6} + f_{7} + f_{8} + \rho u_{1}\right)}{2} - \frac{\omega \left(c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2} - f_{5} - f_{6} - f_{7} - f_{8}\right)}{2} + \frac{\omega \left(c_{s}^{2} \rho - f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{1}^{2}\right)}{2}$$ |
$$d_{3} \leftarrow f_{3} + \frac{\omega \left(c_{s}^{2} \rho u_{0} + f_{5} - f_{6} + f_{7} - f_{8}\right)}{2} - \frac{\omega \left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \rho u_{0}\right)}{2} - \frac{\omega \left(c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2} - f_{5} - f_{6} - f_{7} - f_{8}\right)}{2} + \frac{\omega \left(c_{s}^{2} \rho - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{0}^{2}\right)}{2}$$ |
$$d_{4} \leftarrow f_{4} - \frac{\omega \left(c_{s}^{2} \rho u_{0} + f_{5} - f_{6} + f_{7} - f_{8}\right)}{2} + \frac{\omega \left(f_{3} - f_{4} + f_{5} - f_{6} + f_{7} - f_{8} + \rho u_{0}\right)}{2} - \frac{\omega \left(c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2} - f_{5} - f_{6} - f_{7} - f_{8}\right)}{2} + \frac{\omega \left(c_{s}^{2} \rho - f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{0}^{2}\right)}{2}$$ |
$$d_{5} \leftarrow f_{5} - \frac{\omega \left(f_{5} - f_{6} - f_{7} + f_{8} + \rho u_{0} u_{1}\right)}{4} - \frac{\omega \left(c_{s}^{2} \rho u_{0} + f_{5} - f_{6} + f_{7} - f_{8}\right)}{4} + \frac{\omega \left(c_{s}^{2} \rho u_{1} - f_{5} - f_{6} + f_{7} + f_{8}\right)}{4} + \frac{\omega \left(c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2} - f_{5} - f_{6} - f_{7} - f_{8}\right)}{4}$$ |
$$d_{6} \leftarrow f_{6} + \frac{\omega \left(f_{5} - f_{6} - f_{7} + f_{8} + \rho u_{0} u_{1}\right)}{4} + \frac{\omega \left(c_{s}^{2} \rho u_{0} + f_{5} - f_{6} + f_{7} - f_{8}\right)}{4} + \frac{\omega \left(c_{s}^{2} \rho u_{1} - f_{5} - f_{6} + f_{7} + f_{8}\right)}{4} + \frac{\omega \left(c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2} - f_{5} - f_{6} - f_{7} - f_{8}\right)}{4}$$ |
$$d_{7} \leftarrow f_{7} + \frac{\omega \left(f_{5} - f_{6} - f_{7} + f_{8} + \rho u_{0} u_{1}\right)}{4} - \frac{\omega \left(c_{s}^{2} \rho u_{0} + f_{5} - f_{6} + f_{7} - f_{8}\right)}{4} - \frac{\omega \left(c_{s}^{2} \rho u_{1} - f_{5} - f_{6} + f_{7} + f_{8}\right)}{4} + \frac{\omega \left(c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2} - f_{5} - f_{6} - f_{7} - f_{8}\right)}{4}$$ |
$$d_{8} \leftarrow f_{8} - \frac{\omega \left(f_{5} - f_{6} - f_{7} + f_{8} + \rho u_{0} u_{1}\right)}{4} + \frac{\omega \left(c_{s}^{2} \rho u_{0} + f_{5} - f_{6} + f_{7} - f_{8}\right)}{4} - \frac{\omega \left(c_{s}^{2} \rho u_{1} - f_{5} - f_{6} + f_{7} + f_{8}\right)}{4} + \frac{\omega \left(c_{s}^{4} \rho + c_{s}^{2} \rho u_{0}^{2} + c_{s}^{2} \rho u_{1}^{2} - f_{5} - f_{6} - f_{7} - f_{8}\right)}{4}$$ |
Generic simplification strategy - common subexpresssion elimination¶
[9]:
generic_strategy = ps.simp.SimplificationStrategy()
generic_strategy.add(ps.simp.sympy_cse)
generic_strategy.create_simplification_report(collision_rule)
[9]:
Name | Runtime | Adds | Muls | Divs | Total |
---|---|---|---|---|---|
OriginalTerm | - | 245 | 248 | 2 | 495 |
sympy_cse | 30.74 ms | 82 | 36 | 1 | 119 |
A custom simplification strategy for moment-based methods¶
[10]:
simplification_strategy = create_simplification_strategy(method)
simplification_strategy.create_simplification_report(collision_rule)
simplification_strategy.add(ps.simp.sympy_cse)