Skip to content
Snippets Groups Projects
Commit 78beeb86 authored by Martin Bauer's avatar Martin Bauer
Browse files

Merge branch 'shanchen' into 'master' 

Shan-Chen tutorial

See merge request pycodegen/lbmpy!8
parents 835f8e4d 55512e7d
Branches
No related merge requests found
Hide whitespace changes
Inline Side-by-side
doc/notebooks/07_tutorial_shanchen_twophase.ipynb 0 → 100644
+ 334
− 0
View file @ 78beeb86
%% Cell type:markdown id: tags:
# Shan-Chen Two-Phase Single-Component Lattice Boltzmann
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.updatekernels import create_stream_pull_with_output_kernel
from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
```
%% Cell type:markdown id: tags:
This is based on section 9.3.2 of Krüger et al.'s "The Lattice Boltzmann Method", Springer 2017 (http://www.lbmbook.com).
Sample code is available at [https://github.com/lbm-principles-practice/code/](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp).
%% Cell type:markdown id: tags:
## Parameters
%% Cell type:code id: tags:
``` python
N = 64
omega_a = 1.
g_aa = -4.7
rho0 = 1.
```
%% Cell type:markdown id: tags:
## Data structures
%% Cell type:code id: tags:
``` python
dh = ps.create_data_handling((N, N), periodicity=True, default_target='cpu')
method_a = create_lb_method(relaxation_rate=omega_a, compressible=True)
src_a = dh.add_array('src_a', values_per_cell=len(method_a.stencil))
dst_a = dh.add_array_like('dst_a', 'src_a')
ρ_a = dh.add_array('rho_a')
```
%% Cell type:markdown id: tags:
## Force & combined velocity
%% Cell type:markdown id: tags:
The force on the fluid is
$\vec{F}_A(\vec{x})=-\psi(\rho_A(\vec{x}))g_{AA}\sum\limits_{i=1}^{19}w_i\psi(\rho_A(\vec{x}+\vec{c}_i))\vec{c}_i$
with
$\psi(\rho)=\rho_0\left[1-\exp(-\rho/\rho_0)\right]$.
%% Cell type:code id: tags:
``` python
def psi(dens):
return rho0 * (1. - sp.exp(-dens / rho0));
```
%% Cell type:code id: tags:
``` python
zero_vec = sp.Matrix([0] * dh.dim)
stencil, weights = method_a.stencil, method_a.weights
force_a = sum((psi(ρ_a[d]) * w_d * sp.Matrix(d)
for d, w_d in zip(stencil, weights)),
zero_vec) * psi(ρ_a.center) * -1 * g_aa
```
%% Cell type:markdown id: tags:
## Kernels
%% Cell type:code id: tags:
``` python
stream_a = create_stream_pull_with_output_kernel(method_a, src_a, dst_a, {'density': ρ_a})
# TODO use method above
collision_a = create_lb_update_rule(relaxation_rate=omega_a,
compressible=True,
force_model='guo',
force=force_a,
kernel_type='collide_only',
optimization={'symbolic_field': src_a})
opts = {'cpu_openmp': False}
stream_a_kernel = ps.create_kernel(stream_a, **opts).compile()
collision_a_kernel = ps.create_kernel(collision_a, **opts).compile()
```
%% Cell type:markdown id: tags:
## Initialization
%% Cell type:code id: tags:
``` python
init_a = macroscopic_values_setter(method_a, velocity=(0, 0),
pdfs=src_a.center_vector, density=ρ_a.center)
init_a_kernel = ps.create_kernel(init_a, ghost_layers=0).compile()
```
%% Cell type:code id: tags:
``` python
def init():
for x in range(N):
for y in range(N):
if (x-N/2)**2 + (y-N/2)**2 <= 15**2:
dh.fill(ρ_a.name, 2.1, slice_obj=[x,y])
else:
dh.fill(ρ_a.name, 0.15, slice_obj=[x,y])
dh.run_kernel(init_a_kernel)
```
%% Cell type:markdown id: tags:
## Timeloop
%% Cell type:code id: tags:
``` python
sync_pdfs = dh.synchronization_function([src_a.name])
sync_ρs = dh.synchronization_function([ρ_a.name])
def time_loop(steps):
dh.all_to_gpu()
for i in range(steps):
sync_ρs()
dh.run_kernel(collision_a_kernel)
sync_pdfs()
dh.run_kernel(stream_a_kernel)
dh.swap(src_a.name, dst_a.name)
dh.all_to_cpu()
```
%% Cell type:code id: tags:
``` python
def plot_ρs():
plt.title("$\\rho_A$")
plt.scalar_field(dh.gather_array(ρ_a.name), vmin=0, vmax=2.5)
plt.colorbar()
```
%% Cell type:markdown id: tags:
## Run the simulation
### Initial state
%% Cell type:code id: tags:
``` python
init()
plot_ρs()
```
%% Output
%% Cell type:markdown id: tags:
### Check the first time step against reference data
The reference data was obtained with the [sample code](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp) after making the following changes:
```c++
const int nsteps = 1000;
const int noutput = 1;
```
Remove the next cell if you changed the parameters at the beginning of this notebook.
%% Cell type:code id: tags:
``` python
init()
time_loop(1)
ref_a = np.array([0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.136756, 0.220324, 1.2382, 2.26247, 2.26183, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.26183, 2.26247, 1.2382, 0.220324, 0.136756, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15])
assert np.allclose(dh.gather_array(ρ_a.name)[N//2], ref_a)
```
%% Cell type:markdown id: tags:
### Run the simulation until converged
%% Cell type:code id: tags:
``` python
init()
time_loop(1000)
plot_ρs()
```
%% Output
%% Cell type:code id: tags:
``` python
```
doc/notebooks/08_tutorial_shanchen_twocomponent.ipynb 0 → 100644
+ 392
− 0
View file @ 78beeb86
%% Cell type:markdown id: tags:
# Shan-Chen Two-Component Lattice Boltzmann
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from lbmpy.updatekernels import create_stream_pull_with_output_kernel
from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
```
%% Cell type:markdown id: tags:
This is based on section 9.3.3 of Krüger et al.'s "The Lattice Boltzmann Method", Springer 2017 (http://www.lbmbook.com).
Sample code is available at [https://github.com/lbm-principles-practice/code/](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp).
%% Cell type:markdown id: tags:
## Parameters
%% Cell type:code id: tags:
``` python
N = 64
omega_a = 1.
omega_b = 1.
g_aa = 0.
g_ab = g_ba = 6.
g_bb = 0.
rho0 = 1.
```
%% Cell type:markdown id: tags:
## Data structures
%% Cell type:code id: tags:
``` python
dh = ps.create_data_handling((N, N), periodicity=True, default_target='cpu')
method_a = create_lb_method(relaxation_rate=omega_a, compressible=True)
method_b = create_lb_method(relaxation_rate=omega_b, compressible=True)
src_a = dh.add_array('src_a', values_per_cell=len(method_a.stencil))
dst_a = dh.add_array_like('dst_a', 'src_a')
src_b = dh.add_array('src_b', values_per_cell=len(method_b.stencil))
dst_b = dh.add_array_like('dst_b', 'src_b')
ρ_a = dh.add_array('rho_a')
ρ_b = dh.add_array('rho_b')
u_a = dh.add_array('u_a', values_per_cell=dh.dim)
u_b = dh.add_array('u_b', values_per_cell=dh.dim)
```
%% Cell type:markdown id: tags:
## Force & combined velocity
%% Cell type:markdown id: tags:
The force between the two components is
$\vec{F}_k(\vec{x})=-\psi(\rho_k(\vec{x}))\sum\limits_{k^\prime\in\{A,B\}}g_{kk^\prime}\sum\limits_{i=1}^{19}w_i\psi(\rho_{k^\prime}(\vec{x}+\vec{c}_i))\vec{c}_i$
for $k\in\{A,B\}$
and with
$\psi(\rho)=\rho_0\left[1-\exp(-\rho/\rho_0)\right]$.
%% Cell type:code id: tags:
``` python
def psi(dens):
return rho0 * (1. - sp.exp(-dens / rho0));
```
%% Cell type:code id: tags:
``` python
zero_vec = sp.Matrix([0] * dh.dim)
force_a = zero_vec
for factor, ρ in zip([g_aa, g_ab], [ρ_a, ρ_b]):
stencil, weights = method_a.stencil, method_a.weights
force_a += sum((psi(ρ[d]) * w_d * sp.Matrix(d)
for d, w_d in zip(stencil, weights)),
zero_vec) * psi(ρ_a.center) * -1 * factor
force_b = zero_vec
for factor, ρ in zip([g_ba, g_bb], [ρ_a, ρ_b]):
stencil, weights = method_b.stencil, method_b.weights
force_b += sum((psi(ρ[d]) * w_d * sp.Matrix(d)
for d, w_d in zip(stencil, weights)),
zero_vec) * psi(ρ_b.center) * -1 * factor
```
%% Cell type:markdown id: tags:
The barycentric velocity, which is used in place of the individual components' velocities in the equilibrium distribution and Guo force term, is
$\vec{u}=\frac{1}{\rho_a+\rho_b}\left(\rho_a\vec{u}_a+\frac{1}{2}\vec{F}_a+\rho_b\vec{u}_b+\frac{1}{2}\vec{F}_b\right)$.
%% Cell type:code id: tags:
``` python
u_full = [(ρ_a.center * u_a(i) + force_a[i]/2 +
ρ_b.center * u_b(i) + force_b[i]/2) / (ρ_a.center + ρ_b.center)
for i in range(dh.dim)]
```
%% Cell type:markdown id: tags:
## Kernels
%% Cell type:code id: tags:
``` python
stream_a = create_stream_pull_with_output_kernel(method_a, src_a, dst_a, {'density': ρ_a, 'velocity': u_a})
stream_b = create_stream_pull_with_output_kernel(method_b, src_b, dst_b, {'density': ρ_b, 'velocity': u_b})
# TODO use method above
collision_a = create_lb_update_rule(relaxation_rate=omega_a,
compressible=True,
velocity_input=u_full, density_input=ρ_a,
force_model='guo',
force=force_a,
kernel_type='collide_only',
optimization={'symbolic_field': src_a})
collision_b = create_lb_update_rule(relaxation_rate=omega_b,
compressible=True,
velocity_input=u_full, density_input=ρ_b,
force_model='guo',
force=force_b,
kernel_type='collide_only',
optimization={'symbolic_field': src_b})
opts = {'cpu_openmp': False}
stream_a_kernel = ps.create_kernel(stream_a, **opts).compile()
stream_b_kernel = ps.create_kernel(stream_b, **opts).compile()
collision_a_kernel = ps.create_kernel(collision_a, **opts).compile()
collision_b_kernel = ps.create_kernel(collision_b, **opts).compile()
```
%% Cell type:markdown id: tags:
## Initialization
%% Cell type:code id: tags:
``` python
init_a = macroscopic_values_setter(method_a, velocity=(0, 0),
pdfs=src_a.center_vector, density=ρ_a.center)
init_b = macroscopic_values_setter(method_b, velocity=(0, 0),
pdfs=src_b.center_vector, density=ρ_b.center)
init_a_kernel = ps.create_kernel(init_a, ghost_layers=0).compile()
init_b_kernel = ps.create_kernel(init_b, ghost_layers=0).compile()
```
%% Cell type:code id: tags:
``` python
def init():
dh.fill(ρ_a.name, 0.1, slice_obj=ps.make_slice[:, :0.5])
dh.fill(ρ_a.name, 0.9, slice_obj=ps.make_slice[:, 0.5:])
dh.fill(ρ_b.name, 0.9, slice_obj=ps.make_slice[:, :0.5])
dh.fill(ρ_b.name, 0.1, slice_obj=ps.make_slice[:, 0.5:])
dh.run_kernel(init_a_kernel)
dh.run_kernel(init_b_kernel)
dh.fill(u_a.name, 0.0)
dh.fill(u_b.name, 0.0)
```
%% Cell type:markdown id: tags:
## Timeloop
%% Cell type:code id: tags:
``` python
sync_pdfs = dh.synchronization_function([src_a.name, src_b.name])
sync_ρs = dh.synchronization_function([ρ_a.name, ρ_b.name])
def time_loop(steps):
dh.all_to_gpu()
for i in range(steps):
sync_ρs()
dh.run_kernel(collision_a_kernel)
dh.run_kernel(collision_b_kernel)
sync_pdfs()
dh.run_kernel(stream_a_kernel)
dh.run_kernel(stream_b_kernel)
dh.swap(src_a.name, dst_a.name)
dh.swap(src_b.name, dst_b.name)
dh.all_to_cpu()
```
%% Cell type:code id: tags:
``` python
def plot_ρs():
plt.subplot(1,2,1)
plt.title("$\\rho_A$")
plt.scalar_field(dh.gather_array(ρ_a.name), vmin=0, vmax=2)
plt.colorbar()
plt.subplot(1,2,2)
plt.title("$\\rho_B$")
plt.scalar_field(dh.gather_array(ρ_b.name), vmin=0, vmax=2)
plt.colorbar()
```
%% Cell type:markdown id: tags:
## Run the simulation
### Initial state
%% Cell type:code id: tags:
``` python
init()
plot_ρs()
```
%% Output
%% Cell type:markdown id: tags:
### Check the first time step against reference data
The reference data was obtained with the [sample code](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp) after making the following changes:
```c++
const int nsteps = 1000;
const int noutput = 1;
const int nfluids = 2;
const double gA = 0;
```
Remove the next cell if you changed the parameters at the beginning of this notebook.
%% Cell type:code id: tags:
``` python
init()
time_loop(1)
ref_a = np.array([0.133183, 0.0921801, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.0921801, 0.133183, 0.719568, 1.05507, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 1.05507, 0.719568])
ref_b = np.array([0.719568, 1.05507, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 1.05507, 0.719568, 0.133183, 0.0921801, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.0921801, 0.133183])
assert np.allclose(dh.gather_array(ρ_a.name)[0], ref_a)
assert np.allclose(dh.gather_array(ρ_b.name)[0], ref_b)
```
%% Cell type:markdown id: tags:
### Run the simulation until converged
%% Cell type:code id: tags:
``` python
init()
time_loop(1000)
plot_ρs()
```
%% Output
doc/sphinx/tutorials.rst
+ 2
− 0
View file @ 78beeb86
......@@ -14,6 +14,8 @@ You can open the notebooks directly to play around with the code examples.
/notebooks/04_tutorial_nondimensionalization_and_scaling.ipynb
/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb
/notebooks/06_tutorial_thermal_lbm.ipynb
/notebooks/07_tutorial_shanchen_twophase.ipynb
/notebooks/08_tutorial_shanchen_twocomponent.ipynb
/notebooks/demo_stencils.ipynb
/notebooks/demo_create_method_from_scratch.ipynb
/notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment