import pystencils as ps
import sympy as sp
from lbmpy.session import *
from lbmpy.moments import is_shear_moment, get_order
Demo: Thermalized (Randomized) LBM#
This demo notebook shows how to modify the LB collision operator to account for microscopic fluctuations. This technique is also called thermalized or randomized LBM. In these methods a random fluctuation is added to the equilibrium moments. In this simple example we implement a thermalized model by writing our own simple linear congruent random number generator, the draws indepedent random numbers on each cell. A seed is stored for each cell since all cells are processed in parallel.
1) Creating the method definition#
In the thermalized LBM, the equilibrium values of moments are altered by small, random fluctuations. We achieve this by creating a new class to represent this randomized equilibrium, inheriting from ContinuousHydrodynamicMaxwellian, but overriding the computation of raw moments to insert our random numbers. The overridden equilibrium looks like this:
class ThermalizedEquilibrium(ContinuousHydrodynamicMaxwellian):
def __init__(self, random_number_symbols, **kwargs):
super().__init__(**kwargs)
self.random_number_symbols = random_number_symbols
def moment(self, exponent_tuple_or_polynomial):
value = super().moment(exponent_tuple_or_polynomial)
if is_shear_moment(exponent_tuple_or_polynomial, dim=self.dim):
value += self.random_number_symbols[0] * 0.001
elif get_order(exponent_tuple_or_polynomial) > 2:
value += self.random_number_symbols[1] * 0.001
return value
We use the low-level function create_from_equilibrium to set up a method using this altered equilibrium. This requires also constructing a DensityVelocityComputation instance, and the full collision info dictionary from scratch.
from lbmpy.fluctuatinglb import ThermalizedEquilibrium
from lbmpy.methods import DensityVelocityComputation, CollisionSpaceInfo, create_from_equilibrium
from lbmpy.moments import get_default_moment_set_for_stencil
random_number_symbols = sp.symbols("rand_:3")
stencil = LBStencil(Stencil.D2Q9)
equilibrium = ThermalizedEquilibrium(random_number_symbols, dim=stencil.D, compressible=False, c_s_sq=sp.Rational(1,3))
cqc = DensityVelocityComputation(stencil, False, False)
# SRT-Type Collision Info
moments = get_default_moment_set_for_stencil(stencil)
r_rate = 1.8
r_rate_dict = {m : r_rate for m in moments}
c_space = CollisionSpaceInfo(CollisionSpace.RAW_MOMENTS)
thermalized_method = create_from_equilibrium(stencil, equilibrium, cqc, r_rate_dict, collision_space_info=c_space)
thermalized_method
| 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{3 \delta_{\rho} e^{- \frac{3 v_{0}^{2}}{2} - \frac{3 v_{1}^{2}}{2}}}{2 \pi} + \frac{3 e^{- \frac{3 \left(- u_{0} + v_{0}\right)^{2}}{2} - \frac{3 \left(- u_{1} + v_{1}\right)^{2}}{2}}}{2 \pi}$ | ||
|---|---|---|---|
| Compressible: ✗ | Deviation Only: ✗ | Order: ∞ | |
| Relaxation Info | ||
|---|---|---|
| Moment | Eq. Value | Relaxation Rate |
| $1$ | $\rho$ | $1.8$ |
| $x$ | $u_{0}$ | $1.8$ |
| $y$ | $u_{1}$ | $1.8$ |
| $x^{2}$ | $0.001 rand_{0} + \frac{\rho}{3} + u_{0}^{2}$ | $1.8$ |
| $y^{2}$ | $0.001 rand_{0} + \frac{\rho}{3} + u_{1}^{2}$ | $1.8$ |
| $x y$ | $0.001 rand_{0} + u_{0} u_{1}$ | $1.8$ |
| $x^{2} y$ | $0.001 rand_{1} + u_{0}^{2} u_{1} + \frac{u_{1}}{3}$ | $1.8$ |
| $x y^{2}$ | $0.001 rand_{1} + u_{0} u_{1}^{2} + \frac{u_{0}}{3}$ | $1.8$ |
| $x^{2} y^{2}$ | $0.001 rand_{1} + \frac{\rho}{9} + u_{0}^{2} u_{1}^{2} + \frac{u_{0}^{2}}{3} + \frac{u_{1}^{2}}{3}$ | $1.8$ |
2) Creating the kernel equations#
Next we have to define rules how to compute the quantities \(rand_0\) and \(rand_1\). We do this using a linear congruent RNG on each cell. We pass in a seed array, and in each time step this seed array is updated with the new random numbers.
try:
from pystencils import tcast
except ImportError:
from pystencils.typing import CastFunc as tcast
dh = ps.create_data_handling(domain_size=(80, 80))
seed_type = np.uint32
max_seed_type = np.iinfo(seed_type).max
# Initialize the seed array
seedField = dh.add_array('seed', dtype=seed_type, values_per_cell=len(random_number_symbols))
for b in dh.iterate():
random_field = np.random.randint(0, high=max_seed_type, dtype=seed_type, size=b['seed'].shape)
np.copyto(b['seed'], random_field)
debug_output = dh.add_array('dbg')
linear_congruent_rng_eqs = [ps.Assignment(seedField(i), seedField(i) * 1664525 + 1013904223)
for i, _ in enumerate(random_number_symbols)]
floatEqs = [ps.Assignment(ps.TypedSymbol(s.name, np.float64), tcast(seedField(i) / max_seed_type, np.float64))
for i, s in enumerate(random_number_symbols)]
rng_eqs = linear_congruent_rng_eqs + floatEqs + [ps.Assignment(debug_output.center, tcast(seedField(0) / max_seed_type, np.float64))]
rng_eqs
These assignments are then added to the LB collision rule.
collision_rule = create_lb_collision_rule(lb_method=thermalized_method)
collision_rule.subexpressions = rng_eqs + collision_rule.subexpressions
collision_rule
Subexpressions:
| $${seed}_{(0,0)}^{0} \leftarrow_{} 1664525 {seed}_{(0,0)}^{0} + 1013904223$$ |
| $${seed}_{(0,0)}^{1} \leftarrow_{} 1664525 {seed}_{(0,0)}^{1} + 1013904223$$ |
| $${seed}_{(0,0)}^{2} \leftarrow_{} 1664525 {seed}_{(0,0)}^{2} + 1013904223$$ |
| $$rand_{0} \leftarrow_{} \operatorname{CastFunc}{\left(\frac{{seed}_{(0,0)}^{0}}{4294967295},\text{BasicType} \right)}$$ |
| $$rand_{1} \leftarrow_{} \operatorname{CastFunc}{\left(\frac{{seed}_{(0,0)}^{1}}{4294967295},\text{BasicType} \right)}$$ |
| $$rand_{2} \leftarrow_{} \operatorname{CastFunc}{\left(\frac{{seed}_{(0,0)}^{2}}{4294967295},\text{BasicType} \right)}$$ |
| $${dbg}_{(0,0)} \leftarrow_{} \operatorname{CastFunc}{\left(\frac{{seed}_{(0,0)}^{0}}{4294967295},\text{BasicType} \right)}$$ |
| $$\xi_{1} \leftarrow_{} 0.6$$ |
| $$\xi_{2} \leftarrow_{} 1.8$$ |
| $$\xi_{3} \leftarrow_{} 0.0018$$ |
| $$\xi_{4} \leftarrow_{} 0.8$$ |
| $$\xi_{5} \leftarrow_{} 0.2$$ |
| $$\xi_{6} \leftarrow_{} \frac{1}{2}$$ |
| $$\xi_{7} \leftarrow_{} \frac{1}{4}$$ |
| $$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}$$ |
| $$M_{post 0} \leftarrow_{} 1.0 m_{00}$$ |
| $$M_{post 1} \leftarrow_{} 1.0 m_{10}$$ |
| $$M_{post 2} \leftarrow_{} 1.0 m_{01}$$ |
| $$M_{post 3} \leftarrow_{} m_{00} \xi_{1} + m_{10}^{2} \xi_{2} - m_{20} \xi_{4} + rand_{0} \xi_{3}$$ |
| $$M_{post 4} \leftarrow_{} m_{00} \xi_{1} + m_{01}^{2} \xi_{2} - m_{02} \xi_{4} + rand_{0} \xi_{3}$$ |
| $$M_{post 5} \leftarrow_{} m_{01} m_{10} \xi_{2} - m_{11} \xi_{4} + rand_{0} \xi_{3}$$ |
| $$M_{post 6} \leftarrow_{} m_{01} m_{10}^{2} \xi_{2} + m_{01} \xi_{1} - m_{21} \xi_{4} + rand_{1} \xi_{3}$$ |
| $$M_{post 7} \leftarrow_{} m_{01}^{2} m_{10} \xi_{2} + m_{10} \xi_{1} - m_{12} \xi_{4} + rand_{1} \xi_{3}$$ |
| $$M_{post 8} \leftarrow_{} m_{00} \xi_{5} + m_{01}^{2} m_{10}^{2} \xi_{2} - m_{22} \xi_{4} + rand_{1} \xi_{3} + \xi_{1} \left(m_{01}^{2} + m_{10}^{2}\right)$$ |
| $$sub_{k to f 10} \leftarrow_{} \xi_{6} \left(M_{post 4} - M_{post 8}\right)$$ |
| $$sub_{k to f 11} \leftarrow_{} \xi_{6} \left(M_{post 2} - M_{post 6}\right)$$ |
| $$sub_{k to f 12} \leftarrow_{} \xi_{6} \left(M_{post 3} - M_{post 8}\right)$$ |
| $$sub_{k to f 13} \leftarrow_{} \xi_{6} \left(- M_{post 1} + M_{post 7}\right)$$ |
| $$sub_{k to f 14} \leftarrow_{} \xi_{7} \left(- M_{post 5} + M_{post 8}\right)$$ |
| $$sub_{k to f 15} \leftarrow_{} \xi_{7} \left(M_{post 6} - M_{post 7}\right)$$ |
| $$sub_{k to f 16} \leftarrow_{} \xi_{7} \left(M_{post 5} + M_{post 8}\right)$$ |
| $$sub_{k to f 17} \leftarrow_{} \xi_{7} \left(M_{post 6} + M_{post 7}\right)$$ |
Main Assignments:
| $$d_{0} \leftarrow_{} M_{post 0} - M_{post 3} - M_{post 4} + M_{post 8}$$ |
| $$d_{1} \leftarrow_{} sub_{k to f 10} + sub_{k to f 11}$$ |
| $$d_{2} \leftarrow_{} sub_{k to f 10} - sub_{k to f 11}$$ |
| $$d_{3} \leftarrow_{} sub_{k to f 12} + sub_{k to f 13}$$ |
| $$d_{4} \leftarrow_{} sub_{k to f 12} - sub_{k to f 13}$$ |
| $$d_{5} \leftarrow_{} sub_{k to f 14} + sub_{k to f 15}$$ |
| $$d_{6} \leftarrow_{} sub_{k to f 16} + sub_{k to f 17}$$ |
| $$d_{7} \leftarrow_{} sub_{k to f 16} - sub_{k to f 17}$$ |
| $$d_{8} \leftarrow_{} sub_{k to f 14} - sub_{k to f 15}$$ |
Finally, lets test our method by running a lid-driven-cavity with it.
ldc = create_lid_driven_cavity(data_handling=dh, collision_rule=collision_rule, lid_velocity=0.05,
kernel_params={'rand_0': 0, 'rand_1': 0})
ldc.run(100)
plt.figure(dpi=200)
plt.vector_field(ldc.velocity[:, :]);
assert np.isfinite(dh.max('ldc_velocity'))