Demo: Create lbmpy Method from Scratch

eb810a2144b24d35a0b9b594be115d1c

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()
../_images/notebooks_demo_create_method_from_scratch_4_0.png
[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]$
[6]:
from lbmpy.maxwellian_equilibrium import get_equilibrium_values_of_maxwell_boltzmann_function

eq_moments = get_equilibrium_values_of_maxwell_boltzmann_function(moments, order=2, dim=2,
                                                                  c_s_sq=sp.Rational(1, 3),
                                                                  space="moment")
omega = sp.symbols("omega")
relaxation_info = [(moment, eq_value, omega) for moment, eq_value in zip(moments, eq_moments)]
relaxation_info
[6]:
$\displaystyle \left[ \left( 1, \ \rho, \ \omega\right), \ \left( y, \ \rho u_{1}, \ \omega\right), \ \left( y^{2}, \ \rho u_{1}^{2} + \frac{\rho}{3}, \ \omega\right), \ \left( x, \ \rho u_{0}, \ \omega\right), \ \left( x y, \ \rho u_{0} u_{1}, \ \omega\right), \ \left( x y^{2}, \ \frac{\rho u_{0}}{3}, \ \omega\right), \ \left( x^{2}, \ \rho u_{0}^{2} + \frac{\rho}{3}, \ \omega\right), \ \left( x^{2} y, \ \frac{\rho u_{1}}{3}, \ \omega\right), \ \left( x^{2} y^{2}, \ \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}, \ \omega\right)\right]$
[7]:
from lbmpy.methods.creationfunctions import create_generic_mrt

force_model = forcemodels.Guo(sp.symbols("F_:2"))
method = create_generic_mrt(d2q9, relaxation_info, compressible=False, force_model=force_model)
method
[7]:
Moment Eq. Value Relaxation Rate
$1$ $\rho$ $\omega$
$y$ $\rho u_{1}$ $\omega$
$y^{2}$ $\rho u_{1}^{2} + \frac{\rho}{3}$ $\omega$
$x$ $\rho u_{0}$ $\omega$
$x y$ $\rho u_{0} u_{1}$ $\omega$
$x y^{2}$ $\frac{\rho u_{0}}{3}$ $\omega$
$x^{2}$ $\rho u_{0}^{2} + \frac{\rho}{3}$ $\omega$
$x^{2} y$ $\frac{\rho u_{1}}{3}$ $\omega$
$x^{2} y^{2}$ $\frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}$ $\omega$

Example of a update equation without simplifications

[8]:
collision_rule = method.get_collision_rule()
collision_rule
[8]:
Subexpressions:
$$partial_{m m1 e 0} \leftarrow f_{3} + f_{5} + f_{7}$$
$$partial_{m 0 e 0} \leftarrow f_{0} + f_{1} + f_{2}$$
$$partial_{m 1 e 0} \leftarrow f_{4} + f_{6} + f_{8}$$
$$partial_{m m1 e 1} \leftarrow f_{5} - f_{7}$$
$$partial_{m 0 e 1} \leftarrow f_{1} - f_{2}$$
$$partial_{m 1 e 1} \leftarrow f_{6} - f_{8}$$
$$partial_{m m1 e 2} \leftarrow f_{5} + f_{7}$$
$$partial_{m 0 e 2} \leftarrow f_{1} + f_{2}$$
$$partial_{m 1 e 2} \leftarrow f_{6} + f_{8}$$
$$m_{00} \leftarrow partial_{m 0 e 0} + partial_{m 1 e 0} + partial_{m m1 e 0}$$
$$m_{10} \leftarrow partial_{m 1 e 0} - partial_{m m1 e 0}$$
$$m_{01} \leftarrow partial_{m 0 e 1} + partial_{m 1 e 1} + partial_{m m1 e 1}$$
$$m_{20} \leftarrow partial_{m 1 e 0} + partial_{m m1 e 0}$$
$$m_{02} \leftarrow partial_{m 0 e 2} + partial_{m 1 e 2} + partial_{m m1 e 2}$$
$$m_{11} \leftarrow partial_{m 1 e 1} - partial_{m m1 e 1}$$
$$m_{21} \leftarrow partial_{m 1 e 1} + partial_{m m1 e 1}$$
$$m_{12} \leftarrow partial_{m 1 e 2} - partial_{m m1 e 2}$$
$$m_{22} \leftarrow partial_{m 1 e 2} + partial_{m m1 e 2}$$
$$\rho \leftarrow m_{00}$$
$$u_{0} \leftarrow \frac{F_{0}}{2} + m_{10}$$
$$u_{1} \leftarrow \frac{F_{1}}{2} + m_{01}$$
$$M_{0} \leftarrow m_{00}$$
$$M_{1} \leftarrow m_{01}$$
$$M_{2} \leftarrow m_{02}$$
$$M_{3} \leftarrow m_{10}$$
$$M_{4} \leftarrow m_{11}$$
$$M_{5} \leftarrow m_{12}$$
$$M_{6} \leftarrow m_{20}$$
$$M_{7} \leftarrow m_{21}$$
$$M_{8} \leftarrow m_{22}$$
$$M_{post 0} \leftarrow M_{0} + \omega \left(- M_{0} + \rho\right)$$
$$M_{post 1} \leftarrow F_{1} \left(1 - \frac{\omega}{2}\right) + M_{1} + \omega \left(- M_{1} + \rho u_{1}\right)$$
$$M_{post 2} \leftarrow 2 F_{1} u_{1} \left(1 - \frac{\omega}{2}\right) + M_{2} + \omega \left(- M_{2} + \rho u_{1}^{2} + \frac{\rho}{3}\right)$$
$$M_{post 3} \leftarrow F_{0} \left(1 - \frac{\omega}{2}\right) + M_{3} + \omega \left(- M_{3} + \rho u_{0}\right)$$
$$M_{post 4} \leftarrow M_{4} + \omega \left(- M_{4} + \rho u_{0} u_{1}\right) + \left(1 - \frac{\omega}{2}\right) \left(F_{0} u_{1} + F_{1} u_{0}\right)$$
$$M_{post 5} \leftarrow \frac{F_{0} \left(1 - \frac{\omega}{2}\right)}{3} + M_{5} + \omega \left(- M_{5} + \frac{\rho u_{0}}{3}\right)$$
$$M_{post 6} \leftarrow 2 F_{0} u_{0} \left(1 - \frac{\omega}{2}\right) + M_{6} + \omega \left(- M_{6} + \rho u_{0}^{2} + \frac{\rho}{3}\right)$$
$$M_{post 7} \leftarrow \frac{F_{1} \left(1 - \frac{\omega}{2}\right)}{3} + M_{7} + \omega \left(- M_{7} + \frac{\rho u_{1}}{3}\right)$$
$$M_{post 8} \leftarrow M_{8} + \omega \left(- M_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right) + \left(1 - \frac{\omega}{2}\right) \left(\frac{2 F_{0} u_{0}}{3} + \frac{2 F_{1} u_{1}}{3}\right)$$
$$sub_{k to f 8} \leftarrow \frac{1}{2}$$
$$sub_{k to f 9} \leftarrow \frac{1}{4}$$
$$m_{post 00} \leftarrow M_{post 0}$$
$$m_{post 10} \leftarrow M_{post 3}$$
$$m_{post 01} \leftarrow M_{post 1}$$
$$m_{post 20} \leftarrow M_{post 6}$$
$$m_{post 02} \leftarrow M_{post 2}$$
$$m_{post 11} \leftarrow M_{post 4}$$
$$m_{post 21} \leftarrow M_{post 7}$$
$$m_{post 12} \leftarrow M_{post 5}$$
$$m_{post 22} \leftarrow M_{post 8}$$
$$sub_{k to f 0} \leftarrow sub_{k to f 8} \left(m_{post 02} - m_{post 22}\right)$$
$$sub_{k to f 1} \leftarrow sub_{k to f 8} \left(m_{post 01} - m_{post 21}\right)$$
$$sub_{k to f 2} \leftarrow sub_{k to f 8} \left(m_{post 20} - m_{post 22}\right)$$
$$sub_{k to f 3} \leftarrow sub_{k to f 8} \left(- m_{post 10} + m_{post 12}\right)$$
$$sub_{k to f 4} \leftarrow sub_{k to f 9} \left(- m_{post 11} + m_{post 22}\right)$$
$$sub_{k to f 5} \leftarrow sub_{k to f 9} \left(- m_{post 12} + m_{post 21}\right)$$
$$sub_{k to f 6} \leftarrow sub_{k to f 9} \left(m_{post 11} + m_{post 22}\right)$$
$$sub_{k to f 7} \leftarrow sub_{k to f 9} \left(m_{post 12} + m_{post 21}\right)$$
Main Assignments:
$$d_{0} \leftarrow m_{post 00} - m_{post 02} - m_{post 20} + m_{post 22}$$
$$d_{1} \leftarrow sub_{k to f 0} + sub_{k to f 1}$$
$$d_{2} \leftarrow sub_{k to f 0} - sub_{k to f 1}$$
$$d_{3} \leftarrow sub_{k to f 2} + sub_{k to f 3}$$
$$d_{4} \leftarrow sub_{k to f 2} - sub_{k to f 3}$$
$$d_{5} \leftarrow sub_{k to f 4} + sub_{k to f 5}$$
$$d_{6} \leftarrow sub_{k to f 6} + sub_{k to f 7}$$
$$d_{7} \leftarrow sub_{k to f 6} - sub_{k to f 7}$$
$$d_{8} \leftarrow sub_{k to f 4} - sub_{k to f 5}$$

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]:
NameRuntimeAddsMulsDivsTotal
OriginalTerm- 85 68 0 153
sympy_cse31.54 ms 72 49 0 121

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)