Demo: Create lbmpy Method from Scratch

fa5dc9b658714502856d87b17d2e104e

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]:
from lbmpy.stencils import get_stencil
d2q9 = get_stencil("D2Q9", ordering='walberla')
ps.stencil.plot(d2q9)
../_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_moments_of_continuous_maxwellian_equilibrium

eq_moments = get_moments_of_continuous_maxwellian_equilibrium(moments, order=2, dim=2,
                                                              c_s_sq=sp.Rational(1, 3))
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, cumulant=False)
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:
$$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_{0}}{2} - f_{3} - f_{5} - f_{7} + vel0Term$$
$$u_{1} \leftarrow \frac{F_{1}}{2} - f_{2} + f_{6} - f_{7} - f_{8} + vel1Term$$
$$forceTerm_{0} \leftarrow \left(1 - \frac{\omega}{2}\right) \left(- \frac{4 F_{0} u_{0}}{3} - \frac{4 F_{1} u_{1}}{3}\right)$$
$$forceTerm_{1} \leftarrow \left(1 - \frac{\omega}{2}\right) \left(- \frac{F_{0} u_{0}}{3} + \frac{F_{1} \left(2 u_{1} + 1\right)}{3}\right)$$
$$forceTerm_{2} \leftarrow \left(1 - \frac{\omega}{2}\right) \left(- \frac{F_{0} u_{0}}{3} + \frac{F_{1} \left(2 u_{1} - 1\right)}{3}\right)$$
$$forceTerm_{3} \leftarrow \left(1 - \frac{\omega}{2}\right) \left(\frac{F_{0} \left(2 u_{0} - 1\right)}{3} - \frac{F_{1} u_{1}}{3}\right)$$
$$forceTerm_{4} \leftarrow \left(1 - \frac{\omega}{2}\right) \left(\frac{F_{0} \left(2 u_{0} + 1\right)}{3} - \frac{F_{1} u_{1}}{3}\right)$$
$$forceTerm_{5} \leftarrow \left(1 - \frac{\omega}{2}\right) \left(\frac{F_{0} \left(2 u_{0} - 3 u_{1} - 1\right)}{12} + \frac{F_{1} \left(- 3 u_{0} + 2 u_{1} + 1\right)}{12}\right)$$
$$forceTerm_{6} \leftarrow \left(1 - \frac{\omega}{2}\right) \left(\frac{F_{0} \left(2 u_{0} + 3 u_{1} + 1\right)}{12} + \frac{F_{1} \left(3 u_{0} + 2 u_{1} + 1\right)}{12}\right)$$
$$forceTerm_{7} \leftarrow \left(1 - \frac{\omega}{2}\right) \left(\frac{F_{0} \left(2 u_{0} + 3 u_{1} - 1\right)}{12} + \frac{F_{1} \left(3 u_{0} + 2 u_{1} - 1\right)}{12}\right)$$
$$forceTerm_{8} \leftarrow \left(1 - \frac{\omega}{2}\right) \left(\frac{F_{0} \left(2 u_{0} - 3 u_{1} + 1\right)}{12} + \frac{F_{1} \left(- 3 u_{0} + 2 u_{1} - 1\right)}{12}\right)$$
Main Assignments:
$$d_{0} \leftarrow f_{0} + forceTerm_{0} + \omega \left(- f_{5} - f_{6} - f_{7} - f_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right) - \omega \left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{1}^{2} + \frac{\rho}{3}\right) - \omega \left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{0}^{2} + \frac{\rho}{3}\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} + forceTerm_{1} - \frac{\omega \left(- f_{5} - f_{6} + f_{7} + f_{8} + \frac{\rho u_{1}}{3}\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(- f_{5} - f_{6} - f_{7} - f_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right)}{2} + \frac{\omega \left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{1}^{2} + \frac{\rho}{3}\right)}{2}$$
$$d_{2} \leftarrow f_{2} + forceTerm_{2} + \frac{\omega \left(- f_{5} - f_{6} + f_{7} + f_{8} + \frac{\rho u_{1}}{3}\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(- f_{5} - f_{6} - f_{7} - f_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right)}{2} + \frac{\omega \left(- f_{1} - f_{2} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{1}^{2} + \frac{\rho}{3}\right)}{2}$$
$$d_{3} \leftarrow f_{3} + forceTerm_{3} + \frac{\omega \left(f_{5} - f_{6} + f_{7} - f_{8} + \frac{\rho u_{0}}{3}\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(- f_{5} - f_{6} - f_{7} - f_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right)}{2} + \frac{\omega \left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{0}^{2} + \frac{\rho}{3}\right)}{2}$$
$$d_{4} \leftarrow f_{4} + forceTerm_{4} - \frac{\omega \left(f_{5} - f_{6} + f_{7} - f_{8} + \frac{\rho u_{0}}{3}\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(- f_{5} - f_{6} - f_{7} - f_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right)}{2} + \frac{\omega \left(- f_{3} - f_{4} - f_{5} - f_{6} - f_{7} - f_{8} + \rho u_{0}^{2} + \frac{\rho}{3}\right)}{2}$$
$$d_{5} \leftarrow f_{5} + forceTerm_{5} + \frac{\omega \left(- f_{5} - f_{6} + f_{7} + f_{8} + \frac{\rho u_{1}}{3}\right)}{4} - \frac{\omega \left(f_{5} - f_{6} - f_{7} + f_{8} + \rho u_{0} u_{1}\right)}{4} - \frac{\omega \left(f_{5} - f_{6} + f_{7} - f_{8} + \frac{\rho u_{0}}{3}\right)}{4} + \frac{\omega \left(- f_{5} - f_{6} - f_{7} - f_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right)}{4}$$
$$d_{6} \leftarrow f_{6} + forceTerm_{6} + \frac{\omega \left(- f_{5} - f_{6} + f_{7} + f_{8} + \frac{\rho u_{1}}{3}\right)}{4} + \frac{\omega \left(f_{5} - f_{6} - f_{7} + f_{8} + \rho u_{0} u_{1}\right)}{4} + \frac{\omega \left(f_{5} - f_{6} + f_{7} - f_{8} + \frac{\rho u_{0}}{3}\right)}{4} + \frac{\omega \left(- f_{5} - f_{6} - f_{7} - f_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right)}{4}$$
$$d_{7} \leftarrow f_{7} + forceTerm_{7} - \frac{\omega \left(- f_{5} - f_{6} + f_{7} + f_{8} + \frac{\rho u_{1}}{3}\right)}{4} + \frac{\omega \left(f_{5} - f_{6} - f_{7} + f_{8} + \rho u_{0} u_{1}\right)}{4} - \frac{\omega \left(f_{5} - f_{6} + f_{7} - f_{8} + \frac{\rho u_{0}}{3}\right)}{4} + \frac{\omega \left(- f_{5} - f_{6} - f_{7} - f_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right)}{4}$$
$$d_{8} \leftarrow f_{8} + forceTerm_{8} - \frac{\omega \left(- f_{5} - f_{6} + f_{7} + f_{8} + \frac{\rho u_{1}}{3}\right)}{4} - \frac{\omega \left(f_{5} - f_{6} - f_{7} + f_{8} + \rho u_{0} u_{1}\right)}{4} + \frac{\omega \left(f_{5} - f_{6} + f_{7} - f_{8} + \frac{\rho u_{0}}{3}\right)}{4} + \frac{\omega \left(- f_{5} - f_{6} - f_{7} - f_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\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]:
NameRuntimeAddsMulsDivsTotal
OriginalTerm- 293 261 0 554
sympy_cse38.68 ms 114 67 0 181

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)

Seeing the simplification in action

[11]:
simplification_strategy.show_intermediate_results(collision_rule, symbols=[sp.Symbol("d_7")])
[11]:
Initial Version
$$d_{7} \leftarrow f_{7} + forceTerm_{7} - \frac{\omega \left(- f_{5} - f_{6} + f_{7} + f_{8} + \frac{\rho u_{1}}{3}\right)}{4} + \frac{\omega \left(f_{5} - f_{6} - f_{7} + f_{8} + \rho u_{0} u_{1}\right)}{4} - \frac{\omega \left(f_{5} - f_{6} + f_{7} - f_{8} + \frac{\rho u_{0}}{3}\right)}{4} + \frac{\omega \left(- f_{5} - f_{6} - f_{7} - f_{8} + \frac{\rho u_{0}^{2}}{3} + \frac{\rho u_{1}^{2}}{3} + \frac{\rho}{9}\right)}{4}$$
expand
$$d_{7} \leftarrow - f_{7} \omega + f_{7} + forceTerm_{7} + \frac{\omega \rho u_{0}^{2}}{12} + \frac{\omega \rho u_{0} u_{1}}{4} - \frac{\omega \rho u_{0}}{12} + \frac{\omega \rho u_{1}^{2}}{12} - \frac{\omega \rho u_{1}}{12} + \frac{\omega \rho}{36}$$
replace_second_order_velocity_products
$$d_{7} \leftarrow - f_{7} \omega + f_{7} + forceTerm_{7} + \frac{\omega \rho u_{0}^{2}}{12} - \frac{\omega \rho u_{0}}{12} + \frac{\omega \rho u_{1}^{2}}{12} - \frac{\omega \rho u_{1}}{12} + \frac{\omega \rho \left(u0Pu1^{2} - u_{0}^{2} - u_{1}^{2}\right)}{8} + \frac{\omega \rho}{36}$$
expand
$$d_{7} \leftarrow - f_{7} \omega + f_{7} + forceTerm_{7} + \frac{\omega \rho u0Pu1^{2}}{8} - \frac{\omega \rho u_{0}^{2}}{24} - \frac{\omega \rho u_{0}}{12} - \frac{\omega \rho u_{1}^{2}}{24} - \frac{\omega \rho u_{1}}{12} + \frac{\omega \rho}{36}$$
factor_relaxation_rates
$$d_{7} \leftarrow f_{7} + forceTerm_{7} + \omega \left(- f_{7} + \frac{\rho u0Pu1^{2}}{8} - \frac{\rho u_{0}^{2}}{24} - \frac{\rho u_{0}}{12} - \frac{\rho u_{1}^{2}}{24} - \frac{\rho u_{1}}{12} + \frac{\rho}{36}\right)$$
replace_density_and_velocity
$$d_{7} \leftarrow f_{7} + forceTerm_{7} + \omega \left(- f_{7} + \frac{\rho u0Pu1^{2}}{8} - \frac{\rho u_{0}^{2}}{24} - \frac{\rho u_{0}}{12} - \frac{\rho u_{1}^{2}}{24} - \frac{\rho u_{1}}{12} + \frac{\rho}{36}\right)$$
replace_common_quadratic_and_constant_term
$$d_{7} \leftarrow f_{7} + forceTerm_{7} + \omega \left(- f_{7} + \frac{f_{eq common}}{36} + \frac{\rho u0Pu1^{2}}{8} - \frac{\rho u_{0}}{12} - \frac{\rho u_{1}}{12}\right)$$
factor_density_after_factoring_relaxation_times
$$d_{7} \leftarrow f_{7} + forceTerm_{7} + \omega \left(- f_{7} + \frac{f_{eq common}}{36} + \rho \left(\frac{u0Pu1^{2}}{8} - \frac{u_{0}}{12} - \frac{u_{1}}{12}\right)\right)$$
subexpression_substitution_in_main_assignments
$$d_{7} \leftarrow f_{7} + forceTerm_{7} + \omega \left(- f_{7} + \frac{f_{eq common}}{36} + \rho \left(\frac{u0Pu1^{2}}{8} - \frac{u0Pu1}{12}\right)\right)$$
add_subexpressions_for_divisions
$$d_{7} \leftarrow f_{7} + forceTerm_{7} + \omega \left(- f_{7} + \frac{f_{eq common}}{36} + \rho \left(\frac{u0Pu1^{2}}{8} - \frac{u0Pu1}{12}\right)\right)$$
sympy_cse
$$d_{7} \leftarrow f_{7} + forceTerm_{7} + \omega \left(\rho \left(- \xi_{95} + \xi_{96}\right) + \xi_{64} + \xi_{92}\right)$$