Skip to content
GitLab
Commit 107c35da authored by Michael Kuron's avatar Michael Kuron
Browse files

Merge branch 'warnings' into 'master' 

Fix warnings

See merge request pycodegen/lbmpy!59
parents dd98681c 38a17b0b
Hide whitespace changes
Inline Side-by-side
doc/notebooks/demo_thermalized_lbm.ipynb
View file @ 107c35da
%% Cell type:code id: tags: %% Cell type:code id: tags:
   
``` python ``` python
import pystencils as ps import pystencils as ps
from lbmpy.session import * from lbmpy.session import *
from lbmpy.moments import is_shear_moment, get_order from lbmpy.moments import is_shear_moment, get_order
``` ```
   
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
   
# Demo: Thermalized (Randomized) LBM # 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. 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.
   
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
   
## 1) Creating the method definition ## 1) Creating the method definition
   
We begin with a standard SRT method... We begin with a standard SRT method...
   
%% Cell type:code id: tags: %% Cell type:code id: tags:
   
``` python ``` python
random_number_symbols = sp.symbols("rand_:3") random_number_symbols = sp.symbols("rand_:3")
method = create_lb_method(method='srt', relaxation_rate=1.8) method = create_lb_method(method='srt', relaxation_rate=1.8)
method method
``` ```
   
%% Output %% Output
   
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fee5d180fd0> <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f53ba131fa0>
   
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
   
...and modify its collision table. The `create_lb_method_from_existing` function provides a convenient way to do this. ...and modify its collision table. The `create_lb_method_from_existing` function provides a convenient way to do this.
We pass a custom function that receives a row of the collision table and returns a modified version of it. We pass a custom function that receives a row of the collision table and returns a modified version of it.
Our modification rule adds symbols that stand for random numbers to some of the moments. Our modification rule adds symbols that stand for random numbers to some of the moments.
   
%% Cell type:code id: tags: %% Cell type:code id: tags:
   
``` python ``` python
def modification_func(moment, equilibrium, relaxation_rate): def modification_func(moment, equilibrium, relaxation_rate):
if is_shear_moment(moment, dim=method.dim): if is_shear_moment(moment, dim=method.dim):
equilibrium += random_number_symbols[0] * 0.001 equilibrium += random_number_symbols[0] * 0.001
elif get_order(moment) > 2: elif get_order(moment) > 2:
equilibrium += random_number_symbols[1] * 0.001 equilibrium += random_number_symbols[1] * 0.001
return moment, equilibrium, relaxation_rate return moment, equilibrium, relaxation_rate
   
   
thermalized_method = create_lb_method_from_existing(method, modification_func) thermalized_method = create_lb_method_from_existing(method, modification_func)
thermalized_method.override_weights(method.weights) thermalized_method.override_weights(method.weights)
thermalized_method thermalized_method
``` ```
   
%% Output %% Output
   
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7fee54d8ee80> <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f53b1d12a90>
   
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
   
## 2) Creating the kernel equations ## 2) Creating the kernel equations
   
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
   
Next we have to define rules how to compute the quantities $rand_0$ and $rand_1$. 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. 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.
   
%% Cell type:code id: tags: %% Cell type:code id: tags:
   
``` python ``` python
dh = ps.create_data_handling(domain_size=(80, 80)) dh = ps.create_data_handling(domain_size=(80, 80))
   
seed_type = np.uint32 seed_type = np.uint32
max_seed_type = np.iinfo(seed_type).max max_seed_type = np.iinfo(seed_type).max
   
# Initialize the seed array # Initialize the seed array
seedField = dh.add_array('seed', dtype=seed_type, values_per_cell=len(random_number_symbols)) seedField = dh.add_array('seed', dtype=seed_type, values_per_cell=len(random_number_symbols))
for b in dh.iterate(): for b in dh.iterate():
random_field = np.random.randint(0, high=max_seed_type, dtype=seed_type, size=b['seed'].shape) random_field = np.random.randint(0, high=max_seed_type, dtype=seed_type, size=b['seed'].shape)
np.copyto(b['seed'], random_field) np.copyto(b['seed'], random_field)
   
debug_output = dh.add_array('dbg') debug_output = dh.add_array('dbg')
   
linear_congruent_rng_eqs = [ps.Assignment(seedField(i), seedField(i) * 1664525 + 1013904223) linear_congruent_rng_eqs = [ps.Assignment(seedField(i), seedField(i) * 1664525 + 1013904223)
for i, _ in enumerate(random_number_symbols)] for i, _ in enumerate(random_number_symbols)]
floatEqs = [ps.Assignment(ps.TypedSymbol(s.name, np.float), seedField(i) / max_seed_type) floatEqs = [ps.Assignment(ps.TypedSymbol(s.name, np.float64), seedField(i) / max_seed_type)
for i, s in enumerate(random_number_symbols)] for i, s in enumerate(random_number_symbols)]
   
rng_eqs = linear_congruent_rng_eqs + floatEqs + [ps.Assignment(debug_output.center, seedField(0) / max_seed_type)] rng_eqs = linear_congruent_rng_eqs + floatEqs + [ps.Assignment(debug_output.center, seedField(0) / max_seed_type)]
rng_eqs rng_eqs
``` ```
   
%% Output %% Output
   
$\displaystyle \left[ {{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 \frac{{{seed}_{(0,0)}^{0}}}{4294967295}, \ rand_{1} \leftarrow \frac{{{seed}_{(0,0)}^{1}}}{4294967295}, \ rand_{2} \leftarrow \frac{{{seed}_{(0,0)}^{2}}}{4294967295}, \ {{dbg}_{(0,0)}} \leftarrow \frac{{{seed}_{(0,0)}^{0}}}{4294967295}\right]$ $\displaystyle \left[ {{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 \frac{{{seed}_{(0,0)}^{0}}}{4294967295}, \ rand_{1} \leftarrow \frac{{{seed}_{(0,0)}^{1}}}{4294967295}, \ rand_{2} \leftarrow \frac{{{seed}_{(0,0)}^{2}}}{4294967295}, \ {{dbg}_{(0,0)}} \leftarrow \frac{{{seed}_{(0,0)}^{0}}}{4294967295}\right]$
⎡ ⎡
⎢seed_C__0 := 1664525⋅seed_C__0 + 1013904223, seed_C__1 := 1664525⋅seed_C__1 + ⎢seed_C__0 := 1664525⋅seed_C__0 + 1013904223, seed_C__1 := 1664525⋅seed_C__1 +
⎣ ⎣
seed_C__0 seed_C__0
1013904223, seed_C__2 := 1664525⋅seed_C__2 + 1013904223, rand₀ := ──────────, 1013904223, seed_C__2 := 1664525⋅seed_C__2 + 1013904223, rand₀ := ──────────,
4294967295 4294967295
seed_C__1 seed_C__2 seed_C__0 ⎤ seed_C__1 seed_C__2 seed_C__0 ⎤
rand₁ := ──────────, rand₂ := ──────────, dbg_C := ──────────⎥ rand₁ := ──────────, rand₂ := ──────────, dbg_C := ──────────⎥
4294967295 4294967295 4294967295⎦ 4294967295 4294967295 4294967295⎦
   
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
   
These assignments are then added to the LB collision rule. These assignments are then added to the LB collision rule.
   
%% Cell type:code id: tags: %% Cell type:code id: tags:
   
``` python ``` python
collision_rule = create_lb_collision_rule(lb_method=thermalized_method) collision_rule = create_lb_collision_rule(lb_method=thermalized_method)
collision_rule.subexpressions += rng_eqs collision_rule.subexpressions += rng_eqs
``` ```
   
%% Cell type:markdown id: tags: %% Cell type:markdown id: tags:
   
Finally, lets test our method by running a lid-driven-cavity with it. Finally, lets test our method by running a lid-driven-cavity with it.
   
%% Cell type:code id: tags: %% Cell type:code id: tags:
   
``` python ``` python
ldc = create_lid_driven_cavity(data_handling=dh, collision_rule=collision_rule, lid_velocity=0.05, ldc = create_lid_driven_cavity(data_handling=dh, collision_rule=collision_rule, lid_velocity=0.05,
kernel_params={'rand_0': 0, 'rand_1': 0}) kernel_params={'rand_0': 0, 'rand_1': 0})
``` ```
   
%% Cell type:code id: tags: %% Cell type:code id: tags:
   
``` python ``` python
#show_code(ldc.ast) #show_code(ldc.ast)
``` ```
   
%% Cell type:code id: tags: %% Cell type:code id: tags:
   
``` python ``` python
ldc.run(100) ldc.run(100)
``` ```
   
%% Cell type:code id: tags: %% Cell type:code id: tags:
   
``` python ``` python
plt.figure(dpi=200) plt.figure(dpi=200)
plt.vector_field(ldc.velocity[:, :]); plt.vector_field(ldc.velocity[:, :]);
``` ```
   
%% Output %% Output