Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 6653 additions and 2112 deletions
doc/index.rst
View file @ 32380aab
...@@ -6,6 +6,9 @@ lbmpy ...@@ -6,6 +6,9 @@ lbmpy
:maxdepth: 2 :maxdepth: 2
sphinx/tutorials.rst sphinx/tutorials.rst
sphinx/methods.rst
sphinx/boundary_conditions.rst
sphinx/forcemodels.rst
sphinx/api.rst sphinx/api.rst
......
doc/notebooks/00_tutorial_lbmpy_walberla_overview.ipynb
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/01_tutorial_predefinedScenarios.ipynb
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/02_tutorial_boundary_setup.ipynb
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/03_tutorial_lbm_formulation.ipynb
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/04_tutorial_cumulant_LBM.ipynb 0 → 100644
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/04_tutorial_nondimensionalization_and_scaling.ipynb deleted 100644 → 0
View file @ 9254f55f
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/05_tutorial_modifying_method_smagorinsky.ipynb deleted 100644 → 0
View file @ 9254f55f
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/05_tutorial_nondimensionalization_and_scaling.ipynb 0 → 100644
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb 0 → 100644
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/06_tutorial_thermal_lbm.ipynb deleted 100644 → 0
View file @ 9254f55f
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
from pystencils.timeloop import TimeLoop
```
%% Cell type:markdown id: tags:
# Tutorial 06: Coupling two LBM simulations for thermal simulations
In this notebook we demonstrate how to run a thermal lattice Boltzmann simulation.
We use a separate set of distribution functions to solve the advection-diffusion equation for the temperature.
The zeroth moment of these additional pdfs corresponds to temperature.
The thermal LB step is coupled to the normal hydrodynamic LB scheme using the Boussinesq approximation. The force on the liquid is proportional to the relative temperature. The hydrodynamic LB method computes the fluid velocity which in turn enters the thermal scheme, completing the two-way coupling.
To set this up in *lbmpy* we create first a `data handling` object and create an array to store the temperature.
%% Cell type:code id: tags:
``` python
domain_size = (100, 50)
gpu = False
dh = ps.create_data_handling(domain_size)
temperature_field = dh.add_array("T", gpu=gpu)
dh.fill('T', val=1.0)
```
%% Cell type:markdown id: tags:
Next, we define how to compute the local force from the temperature field:
%% Cell type:code id: tags:
``` python
gravity = sp.Matrix([0, -1e-2])
force = -gravity * (temperature_field(0) - 1.0)
```
%% Cell type:markdown id: tags:
Now, we can create both LB steps.
The coupling is created by passing the following parameters to the hydrodynamic step,
- `compute_velocity_in_every_step`: usually the velocity is not computed/stored in every time step, only for output and plotting reasons. In the coupled algorithm we have to make sure that after every time step the current velocity is stored in an array, since it enters the thermal LB scheme.
- `force`: we can simply pass our sympy expression for the force here, as long as the `LatticeBoltzmannStep` operates on a data handling that stores the fields that are referenced in the expression. This is why we have to create the data handling first and pass it to both Step objects
and to the thermal step
- `compute_density_in_every_step`: density corresponds to the temperature here, which we need to be computed in every time step, since it enters the force expression for the hydrodynamic scheme
- `equilibrium_order`: for the thermal LB method a first order accurate equilibrium is sufficient. This is slightly faster to compute than the normal equilibrium of order 2
- `velocity_input_array_name`: the velocity entering the thermal equilibrium equation is not computed as first moment of the thermal pdfs. Instead, the hydrodynamic velocity is used here.
%% Cell type:code id: tags:
``` python
optimization = {'target': 'cpu' if gpu else 'cpu', 'openmp': 2}
hydro_step = LatticeBoltzmannStep(data_handling=dh, name='hydro', optimization=optimization,
relaxation_rate=1.8,
compute_velocity_in_every_step=True,
force=force)
thermal_step = LatticeBoltzmannStep(data_handling=dh, name='thermal', optimization=optimization,
relaxation_rate=1.8, density_data_name="T",
compute_density_in_every_step=True,
equilibrium_order=1,
velocity_input_array_name=hydro_step.velocity_data_name)
```
%% Cell type:markdown id: tags:
We add `NoSlip` boundary conditions on all four walls for both schemes with the exception of the left wall of the thermal scheme. Here we set a heated wall, where the temperature is fixed to `1.01`. This kind of Dirichlet boundary condition can be set using a `FixedDensity` LB boundary.
%% Cell type:code id: tags:
``` python
add_box_boundary(hydro_step.boundary_handling)
add_box_boundary(thermal_step.boundary_handling)
thermal_step.boundary_handling.set_boundary(FixedDensity(1.01), slice_from_direction('W', dh.dim))
plt.figure(figsize=(20, 5))
plt.subplot(1, 2, 1)
plt.title("Hydrodynamic")
plt.boundary_handling(hydro_step.boundary_handling)
plt.subplot(1, 2, 2)
plt.title("Thermal")
plt.boundary_handling(thermal_step.boundary_handling)
```
%% Output
%% Cell type:code id: tags:
``` python
def run(time_steps):
hydro_step.pre_run()
thermal_step.pre_run()
for t in range(time_steps):
hydro_step.time_step()
thermal_step.time_step()
hydro_step.post_run()
thermal_step.post_run()
hydro_step.time_steps_run += time_steps
thermal_step.time_steps_run += time_steps
```
%% Cell type:code id: tags:
``` python
run(5000)
assert dh.max('T') < 2
plt.figure(figsize=(16, 5))
plt.subplot(1, 2, 1)
plt.scalar_field(thermal_step.density[:, :])
plt.subplot(1, 2, 2)
plt.vector_field(hydro_step.velocity[:, :]);
```
%% Output
%% Cell type:code id: tags:
``` python
assert np.isfinite(hydro_step.velocity[:, :].max())
assert np.isfinite(thermal_step.density[:, :].max())
```
doc/notebooks/07_tutorial_shanchen_twophase.ipynb deleted 100644 → 0
View file @ 9254f55f
%% 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
from lbmpy.maxwellian_equilibrium import get_weights
```
%% 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.
stencil = get_stencil("D2Q9")
weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))
```
%% Cell type:markdown id: tags:
## Data structures
%% Cell type:code id: tags:
``` python
dim = len(stencil[0])
dh = ps.create_data_handling((N,)*dim, periodicity=True, default_target='cpu')
src = dh.add_array('src', values_per_cell=len(stencil))
dst = dh.add_array_like('dst', 'src')
ρ = dh.add_array('rho')
```
%% 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)
force = sum((psi(ρ[d]) * w_d * sp.Matrix(d)
for d, w_d in zip(stencil, weights)), zero_vec) * psi(ρ.center) * -1 * g_aa
```
%% Cell type:markdown id: tags:
## Kernels
%% Cell type:code id: tags:
``` python
collision = create_lb_update_rule(stencil=stencil,
relaxation_rate=omega_a,
compressible=True,
force_model='guo',
force=force,
kernel_type='collide_only',
optimization={'symbolic_field': src})
stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'density': ρ})
opts = {'cpu_openmp': False,
'target': dh.default_target}
stream_kernel = ps.create_kernel(stream, **opts).compile()
collision_kernel = ps.create_kernel(collision, **opts).compile()
```
%% Cell type:markdown id: tags:
## Initialization
%% Cell type:code id: tags:
``` python
method_without_force = create_lb_method(stencil=stencil, relaxation_rate=omega_a, compressible=True)
init_assignments = macroscopic_values_setter(method_without_force, velocity=(0, 0),
pdfs=src.center_vector, density=ρ.center)
init_kernel = ps.create_kernel(init_assignments, 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(ρ.name, 2.1, slice_obj=[x,y])
else:
dh.fill(ρ.name, 0.15, slice_obj=[x,y])
dh.run_kernel(init_kernel)
```
%% Cell type:markdown id: tags:
## Timeloop
%% Cell type:code id: tags:
``` python
sync_pdfs = dh.synchronization_function([src.name])
sync_ρs = dh.synchronization_function([ρ.name])
def time_loop(steps):
dh.all_to_gpu()
for i in range(steps):
sync_ρs()
dh.run_kernel(collision_kernel)
sync_pdfs()
dh.run_kernel(stream_kernel)
dh.swap(src.name, dst.name)
dh.all_to_cpu()
```
%% Cell type:code id: tags:
``` python
def plot_ρs():
plt.title("$\\rho$")
plt.scalar_field(dh.gather_array(ρ.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 = 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(ρ.name)[N//2], ref)
```
%% Cell type:markdown id: tags:
### Run the simulation until converged
%% Cell type:code id: tags:
``` python
init()
time_loop(1000)
plot_ρs()
```
%% Output
doc/notebooks/07_tutorial_thermal_lbm.ipynb 0 → 100644
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/08_tutorial_shanchen_twocomponent.ipynb deleted 100644 → 0
View file @ 9254f55f
%% 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
from lbmpy.maxwellian_equilibrium import get_weights
```
%% 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 # domain size
omega_a = 1. # relaxation rate of first component
omega_b = 1. # relaxation rate of second component
# interaction strength
g_aa = 0.
g_ab = g_ba = 6.
g_bb = 0.
rho0 = 1.
stencil = get_stencil("D2Q9")
weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))
```
%% Cell type:markdown id: tags:
## Data structures
We allocate two sets of PDF's, one for each phase. Additionally, for each phase there is one field to store its density and velocity.
To run the simulation on GPU, change the `default_target` to gpu
%% Cell type:code id: tags:
``` python
dim = len(stencil[0])
dh = ps.create_data_handling((N, ) * dim, periodicity=True, default_target='cpu')
src_a = dh.add_array('src_a', values_per_cell=len(stencil))
dst_a = dh.add_array_like('dst_a', 'src_a')
src_b = dh.add_array('src_b', values_per_cell=len(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
The two LB methods are coupled using a force term. Its symbolic representation is created in the next cells.
The force value is not written to a field, but directly evaluated inside the collision kernel.
%% 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]):
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]):
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
collision_a = create_lb_update_rule(stencil=stencil,
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(stencil=stencil,
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})
```
%% Cell type:code id: tags:
``` python
stream_a = create_stream_pull_with_output_kernel(collision_a.method, src_a, dst_a,
{'density': ρ_a, 'velocity': u_a})
stream_b = create_stream_pull_with_output_kernel(collision_b.method, src_b, dst_b,
{'density': ρ_b, 'velocity': u_b})
opts = {'cpu_openmp': 1, # number of threads when running on CPU
'target': dh.default_target}
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(collision_a.method, velocity=(0, 0),
pdfs=src_a.center_vector, density=ρ_a.center)
init_b = macroscopic_values_setter(collision_b.method, 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() # collision kernel evaluates force values, that depend on neighboring ρ'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/notebooks/08_tutorial_shanchen_twophase.ipynb 0 → 100644
View file @ 32380aab
%% 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
from lbmpy.maxwellian_equilibrium import get_weights
```
%% 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.
stencil = LBStencil(Stencil.D2Q9)
weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))
```
%% Cell type:markdown id: tags:
## Data structures
%% Cell type:code id: tags:
``` python
dh = ps.create_data_handling((N,) * stencil.D, periodicity=True, default_target=ps.Target.CPU)
src = dh.add_array('src', values_per_cell=stencil.Q)
dst = dh.add_array_like('dst', 'src')
ρ = dh.add_array('rho')
```
%% 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}^{q}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] * stencil.D)
force = sum((psi(ρ[d]) * w_d * sp.Matrix(d)
for d, w_d in zip(stencil, weights)), zero_vec) * psi(ρ.center) * -1 * g_aa
```
%% Cell type:markdown id: tags:
## Kernels
%% Cell type:code id: tags:
``` python
lbm_config = LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True,
force_model=ForceModel.GUO, force=force, kernel_type='collide_only')
collision = create_lb_update_rule(lbm_config=lbm_config,
optimization={'symbolic_field': src})
stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'density': ρ})
config = ps.CreateKernelConfig(target=dh.default_target, cpu_openmp=False)
stream_kernel = ps.create_kernel(stream, config=config).compile()
collision_kernel = ps.create_kernel(collision, config=config).compile()
```
%% Cell type:markdown id: tags:
## Initialization
%% Cell type:code id: tags:
``` python
method_without_force = create_lb_method(LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True))
init_assignments = macroscopic_values_setter(method_without_force, velocity=(0, 0),
pdfs=src.center_vector, density=ρ.center)
init_kernel = ps.create_kernel(init_assignments, ghost_layers=0, config=config).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(ρ.name, 2.1, slice_obj=[x,y])
else:
dh.fill(ρ.name, 0.15, slice_obj=[x,y])
dh.run_kernel(init_kernel)
```
%% Cell type:markdown id: tags:
## Timeloop
%% Cell type:code id: tags:
``` python
sync_pdfs = dh.synchronization_function([src.name])
sync_ρs = dh.synchronization_function([ρ.name])
def time_loop(steps):
dh.all_to_gpu()
for i in range(steps):
sync_ρs()
dh.run_kernel(collision_kernel)
sync_pdfs()
dh.run_kernel(stream_kernel)
dh.swap(src.name, dst.name)
dh.all_to_cpu()
```
%% Cell type:code id: tags:
``` python
def plot_ρs():
plt.figure(dpi=200)
plt.title("$\\rho$")
plt.scalar_field(dh.gather_array(ρ.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:
### Run the simulation until converged
%% Cell type:code id: tags:
``` python
init()
time_loop(1000)
plot_ρs()
```
%% Output
%% Cell type:code id: tags:
``` python
assert np.isfinite(dh.gather_array(ρ.name)).all()
```
doc/notebooks/09_tutorial_shanchen_twocomponent.ipynb 0 → 100644
View file @ 32380aab
%% 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
from lbmpy.maxwellian_equilibrium import get_weights
```
%% 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 # domain size
omega_a = 1. # relaxation rate of first component
omega_b = 1. # relaxation rate of second component
# interaction strength
g_aa = 0.
g_ab = g_ba = 6.
g_bb = 0.
rho0 = 1.
stencil = LBStencil(Stencil.D2Q9)
weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
```
%% Cell type:markdown id: tags:
## Data structures
We allocate two sets of PDF's, one for each phase. Additionally, for each phase there is one field to store its density and velocity.
To run the simulation on GPU, change the `default_target` to gpu
%% Cell type:code id: tags:
``` python
dim = stencil.D
dh = ps.create_data_handling((N, ) * dim, periodicity=True, default_target=ps.Target.CPU)
src_a = dh.add_array('src_a', values_per_cell=stencil.Q)
dst_a = dh.add_array_like('dst_a', 'src_a')
src_b = dh.add_array('src_b', values_per_cell=stencil.Q)
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=stencil.D)
u_b = dh.add_array('u_b', values_per_cell=stencil.D)
u_bary = dh.add_array_like('u_bary', u_a.name)
f_a = dh.add_array('f_a', values_per_cell=stencil.D)
f_b = dh.add_array_like('f_b', f_a.name)
```
%% Cell type:markdown id: tags:
## Force & combined velocity
The two LB methods are coupled using a force term. Its symbolic representation is created in the next cells.
The force value is not written to a field, but directly evaluated inside the collision kernel.
%% Cell type:markdown id: tags:
The force between the two components is
$\mathbf{F}_k(\mathbf{x})=-\psi(\rho_k(\mathbf{x}))\sum\limits_{k^\prime\in\{A,B\}}g_{kk^\prime}\sum\limits_{i=1}^{q}w_i\psi(\rho_{k^\prime}(\mathbf{x}+\mathbf{c}_i))\mathbf{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]):
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]):
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
f_expressions = ps.AssignmentCollection([
ps.Assignment(f_a.center_vector, force_a),
ps.Assignment(f_b.center_vector, force_b)
])
```
%% 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 = list(ps.Assignment(u_bary.center_vector,
(ρ_a.center * u_a.center_vector + ρ_b.center * u_b.center_vector) / (ρ_a.center + ρ_b.center)))
```
%% Cell type:markdown id: tags:
## Kernels
%% Cell type:code id: tags:
``` python
lbm_config_a = LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True,
velocity_input=u_bary, density_input=ρ_a, force_model=ForceModel.GUO,
force=f_a, kernel_type='collide_only')
lbm_config_b = LBMConfig(stencil=stencil, relaxation_rate=omega_b, compressible=True,
velocity_input=u_bary, density_input=ρ_b, force_model=ForceModel.GUO,
force=f_b, kernel_type='collide_only')
collision_a = create_lb_update_rule(lbm_config=lbm_config_a,
optimization={'symbolic_field': src_a})
collision_b = create_lb_update_rule(lbm_config=lbm_config_b,
optimization={'symbolic_field': src_b})
```
%% Cell type:code id: tags:
``` python
stream_a = create_stream_pull_with_output_kernel(collision_a.method, src_a, dst_a,
{'density': ρ_a, 'velocity': u_a})
stream_b = create_stream_pull_with_output_kernel(collision_b.method, src_b, dst_b,
{'density': ρ_b, 'velocity': u_b})
config = ps.CreateKernelConfig(target=dh.default_target)
stream_a_kernel = ps.create_kernel(stream_a, config=config).compile()
stream_b_kernel = ps.create_kernel(stream_b, config=config).compile()
collision_a_kernel = ps.create_kernel(collision_a, config=config).compile()
collision_b_kernel = ps.create_kernel(collision_b, config=config).compile()
force_kernel = ps.create_kernel(f_expressions, config=config).compile()
u_kernel = ps.create_kernel(u_full, config=config).compile()
```
%% Cell type:markdown id: tags:
## Initialization
%% Cell type:code id: tags:
``` python
init_a = macroscopic_values_setter(collision_a.method, velocity=(0, 0),
pdfs=src_a.center_vector, density=ρ_a.center)
init_b = macroscopic_values_setter(collision_b.method, 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()
dh.run_kernel(init_a_kernel)
dh.run_kernel(init_b_kernel)
```
%% 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.fill(f_a.name, 0.0)
dh.fill(f_b.name, 0.0)
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() # force values depend on neighboring ρ's
dh.run_kernel(force_kernel)
dh.run_kernel(u_kernel)
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.figure(dpi=200)
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:
### Run the simulation until converged
%% Cell type:code id: tags:
``` python
init()
time_loop(10000)
plot_ρs()
```
%% Output
%% Cell type:code id: tags:
``` python
assert np.isfinite(dh.gather_array(ρ_a.name)).all()
assert np.isfinite(dh.gather_array(ρ_b.name)).all()
```
doc/notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb 0 → 100644
View file @ 32380aab
%% Cell type:markdown id: tags:
# The conservative Allen-Cahn model for high Reynolds number, two phase flow with large-density and viscosity constrast
%% Cell type:code id: tags:
``` python
from pystencils.session import *
from lbmpy.session import *
from pystencils.boundaries import BoundaryHandling
from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
from lbmpy.phasefield_allen_cahn.kernel_equations import *
from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti, AllenCahnParameters
from lbmpy.advanced_streaming import LBMPeriodicityHandling
from lbmpy.boundaries import NoSlip, LatticeBoltzmannBoundaryHandling
```
%% Cell type:markdown id: tags:
If `cupy` is installed the simulation automatically runs on GPU
%% Cell type:code id: tags:
``` python
try:
import cupy
except ImportError:
cupy = None
gpu = False
target = ps.Target.CPU
print('No cupy installed')
if cupy:
gpu = True
target = ps.Target.GPU
```
%% Cell type:markdown id: tags:
The conservative Allen-Cahn model (CACM) for two-phase flow is based on the work of Fakhari et al. (2017) [Improved locality of the phase-field lattice-Boltzmann model for immiscible fluids at high density ratios](http://dx.doi.org/10.1103/PhysRevE.96.053301). The model can be created for two-dimensional problems as well as three-dimensional problems, which have been described by Mitchell et al. (2018) [Development of a three-dimensional
phase-field lattice Boltzmann method for the study of immiscible fluids at high density ratios](http://dx.doi.org/10.1103/PhysRevE.96.053301). Furthermore, cascaded lattice Boltzmann methods can be combined with the model which was described in [A cascaded phase-field lattice Boltzmann model for the simulation of incompressible, immiscible fluids with high density contrast](http://dx.doi.org/10.1016/j.camwa.2019.08.018)
The CACM is suitable for simulating highly complex two phase flow problems with high-density ratios and high Reynolds numbers. In this tutorial, an overview is provided on how to derive the model with lbmpy. For this, the model is defined with two LBM populations. One for the interface tracking, which we call the phase-field LB step and one for recovering the hydrodynamic properties. The latter is called the hydrodynamic LB step.
%% Cell type:markdown id: tags:
## Geometry Setup
First of all, the stencils for the phase-field LB step as well as the stencil for the hydrodynamic LB step are defined. According to the stencils, the simulation can be performed in either 2D- or 3D-space. For 2D simulations, only the D2Q9 stencil is supported. For 3D simulations, the D3Q15, D3Q19 and the D3Q27 stencil are supported. Note here that the cascaded LBM can not be derived for D3Q15 stencils.
%% Cell type:code id: tags:
``` python
stencil_phase = LBStencil(Stencil.D2Q9)
stencil_hydro = LBStencil(Stencil.D2Q9)
assert(len(stencil_phase[0]) == len(stencil_hydro[0]))
dimensions = len(stencil_phase[0])
```
%% Cell type:markdown id: tags:
Definition of the domain size
%% Cell type:code id: tags:
``` python
# domain
L0 = 256
domain_size = (L0, 4 * L0)
```
%% Cell type:markdown id: tags:
## Parameter definition
The next step is to calculate all parameters which are needed for the simulation. In this example, a Rayleigh-Taylor instability test case is set up. The parameter calculation for this setup is already implemented in lbmpy and can be used with the dimensionless parameters which describe the problem.
%% Cell type:code id: tags:
``` python
# time step
timesteps = 8000
# reference time
reference_time = 4000
# calculate the parameters for the RTI
parameters = calculate_parameters_rti(reference_length=L0,
reference_time=reference_time,
density_heavy=1.0,
capillary_number=0.44,
reynolds_number=3000,
atwood_number=0.998,
peclet_number=1000,
density_ratio=1000,
viscosity_ratio=100)
```
%% Cell type:markdown id: tags:
This function returns a `AllenCahnParameters` class. It is struct like class holding all parameters for the conservative Allen Cahn model:
%% Cell type:code id: tags:
``` python
parameters
```
%% Output
<lbmpy.phasefield_allen_cahn.parameter_calculation.AllenCahnParameters at 0x1329b88b0>
%% Cell type:markdown id: tags:
## Fields
As a next step all fields which are needed get defined. To do so, we create a `datahandling` object. More details about it can be found in the third tutorial of the [pystencils framework]( http://pycodegen.pages.walberla.net/pystencils/). This object holds all fields and manages the kernel runs.
%% Cell type:code id: tags:
``` python
# create a datahandling object
dh = ps.create_data_handling((domain_size), periodicity=(True, False), parallel=False, default_target=target)
# pdf fields. g is used for hydrodynamics and h for the interface tracking
g = dh.add_array("g", values_per_cell=len(stencil_hydro))
dh.fill("g", 0.0, ghost_layers=True)
h = dh.add_array("h",values_per_cell=len(stencil_phase))
dh.fill("h", 0.0, ghost_layers=True)
g_tmp = dh.add_array("g_tmp", values_per_cell=len(stencil_hydro))
dh.fill("g_tmp", 0.0, ghost_layers=True)
h_tmp = dh.add_array("h_tmp",values_per_cell=len(stencil_phase))
dh.fill("h_tmp", 0.0, ghost_layers=True)
# velocity field
u = dh.add_array("u", values_per_cell=dh.dim)
dh.fill("u", 0.0, ghost_layers=True)
# phase-field
C = dh.add_array("C")
dh.fill("C", 0.0, ghost_layers=True)
C_tmp = dh.add_array("C_tmp")
dh.fill("C_tmp", 0.0, ghost_layers=True)
```
%% Cell type:markdown id: tags:
As a next step the relaxation time is stated in a symbolic form. It is calculated via interpolation.
%% Cell type:code id: tags:
``` python
rho_L = parameters.symbolic_density_light
rho_H = parameters.symbolic_density_heavy
density = rho_L + C.center * (rho_H - rho_L)
body_force = [0, 0, 0]
body_force[1] = parameters.symbolic_gravitational_acceleration * density
```
%% Cell type:markdown id: tags:
## Definition of the lattice Boltzmann methods
%% Cell type:markdown id: tags:
For both LB steps, a weighted orthogonal MRT (WMRT) method is used. It is also possible to change the method to a simpler SRT scheme or a more complicated CLBM scheme. The CLBM scheme can be obtained with `Method.CENTRAL_MOMENT`. Note here that the hydrodynamic LB step is formulated as an incompressible velocity-based LBM. Thus, the velocity terms can not be removed from the equilibrium in the central moment space.
%% Cell type:code id: tags:
``` python
w_c = parameters.symbolic_omega_phi
config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
delta_equilibrium=False,
force=sp.symbols("F_:2"), velocity_input=u,
weighted=True, relaxation_rates=[0, w_c, w_c, 1, 1, 1, 1, 1, 1],
output={'density': C_tmp})
method_phase = create_lb_method(lbm_config=config_phase)
method_phase
```
%% Output
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x1346b22e0>
%% Cell type:code id: tags:
``` python
omega = parameters.omega(C)
config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
weighted=True, relaxation_rates=[omega, 1, 1, 1],
force=sp.symbols("F_:2"), output={'velocity': u})
method_hydro = create_lb_method(lbm_config=config_hydro)
method_hydro
```
%% Output
<lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x137772a30>
%% Cell type:markdown id: tags:
## Initialization
%% Cell type:markdown id: tags:
The probability distribution functions (pdfs) are initialised with the equilibrium distribution for the LB methods.
%% Cell type:code id: tags:
``` python
h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)
h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()
g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()
```
%% Cell type:markdown id: tags:
Following this, the phase field is initialised directly in python.
%% Cell type:code id: tags:
``` python
# initialize the domain
def Initialize_distributions():
Nx = domain_size[0]
Ny = domain_size[1]
for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
x = np.zeros_like(block.midpoint_arrays[0])
x[:, :] = block.midpoint_arrays[0]
y = np.zeros_like(block.midpoint_arrays[1])
y[:, :] = block.midpoint_arrays[1]
y -= 2 * L0
tmp = 0.1 * Nx * np.cos((2 * np.pi * x) / Nx)
init_values = 0.5 + 0.5 * np.tanh((y - tmp) / (parameters.interface_thickness / 2))
block["C"][:, :] = init_values
block["C_tmp"][:, :] = init_values
if gpu:
dh.all_to_gpu()
dh.run_kernel(h_init, **parameters.symbolic_to_numeric_map)
dh.run_kernel(g_init)
```
%% Cell type:code id: tags:
``` python
Initialize_distributions()
plt.scalar_field(dh.gather_array(C.name))
```
%% Output
<matplotlib.image.AxesImage at 0x137ad6bb0>
%% Cell type:markdown id: tags:
## Source Terms
%% Cell type:markdown id: tags:
For the Allen-Cahn LB step, the Allen-Cahn equation needs to be applied as a source term. Here, a simple forcing model is used which is directly applied in the moment space:
$$
F_i^\phi (\boldsymbol{x}, t) = \Delta t \frac{\left[1 - 4 \left(\phi - \phi_0\right)^2\right]}{\xi} w_i \boldsymbol{c}_i \cdot \frac{\nabla \phi}{|{\nabla \phi}|},
$$
where $\phi$ is the phase-field, $\phi_0$ is the interface location, $\Delta t$ it the timestep size $\xi$ is the interface width, $\boldsymbol{c}_i$ is the discrete direction from stencil_phase and $w_i$ are the weights. Furthermore, the equilibrium needs to be shifted:
$$
\bar{h}^{eq}_\alpha = h^{eq}_\alpha - \frac{1}{2} F^\phi_\alpha
$$
The hydrodynamic force is given by:
$$
F_i (\boldsymbol{x}, t) = \Delta t w_i \frac{\boldsymbol{c}_i \boldsymbol{F}}{\rho c_s^2},
$$
where $\rho$ is the interpolated density and $\boldsymbol{F}$ is the source term which consists of the pressure force
$$
\boldsymbol{F}_p = -p^* c_s^2 \nabla \rho,
$$
the surface tension force:
$$
\boldsymbol{F}_s = \mu_\phi \nabla \phi
$$
and the viscous force term:
$$
F_{\mu, i}^{\mathrm{MRT}} = - \frac{\nu}{c_s^2 \Delta t} \left[\sum_{\beta} c_{\beta i} c_{\beta j} \times \sum_{\alpha} \Omega_{\beta \alpha}(g_\alpha - g_\alpha^{\mathrm{eq}})\right] \frac{\partial \rho}{\partial x_j}.
$$
In the above equations $p^*$ is the normalised pressure which can be obtained from the zeroth order moment of the hydrodynamic distribution function $g$. The lattice speed of sound is given with $c_s$ and the chemical potential is $\mu_\phi$. Furthermore, the viscosity is $\nu$ and $\Omega$ is the moment-based collision operator. Note here that the hydrodynamic equilibrium is also adjusted as shown above for the phase-field distribution functions.
For CLBM methods the forcing is applied directly in the central moment space. This is done with the `CentralMomentMultiphaseForceModel`. Furthermore, the GUO force model is applied here to be consistent with [A cascaded phase-field lattice Boltzmann model for the simulation of incompressible, immiscible fluids with high density contrast](http://dx.doi.org/10.1016/j.camwa.2019.08.018). Here we refer to equation D.7 which can be derived for 3D stencils automatically with lbmpy.
%% Cell type:code id: tags:
``` python
force_h = interface_tracking_force(C, stencil_phase, parameters)
hydro_force = hydrodynamic_force(C, method_hydro, parameters, body_force)
```
%% Cell type:markdown id: tags:
## Definition of the LB update rules
%% Cell type:markdown id: tags:
The update rule for the phase-field LB step is defined as:
$$
h_i (\boldsymbol{x} + \boldsymbol{c}_i \Delta t, t + \Delta t) = h_i(\boldsymbol{x}, t) + \Omega_{ij}^h(\bar{h_j}^{eq} - h_j)|_{(\boldsymbol{x}, t)} + F_i^\phi(\boldsymbol{x}, t).
$$
In our framework the pull scheme is applied as streaming step. Furthermore, the update of the phase-field is directly integrated into the kernel. As a result of this, a second temporary phase-field is needed.
%% Cell type:code id: tags:
``` python
lbm_optimisation = LBMOptimisation(symbolic_field=h, symbolic_temporary_field=h_tmp)
allen_cahn_update_rule = create_lb_update_rule(lbm_config=config_phase,
lbm_optimisation=lbm_optimisation)
allen_cahn_update_rule = add_interface_tracking_force(allen_cahn_update_rule, force_h)
ast_kernel = ps.create_kernel(allen_cahn_update_rule, target=dh.default_target, cpu_openmp=True)
kernel_allen_cahn_lb = ast_kernel.compile()
```
%% Cell type:markdown id: tags:
The update rule for the hydrodynmaic LB step is defined as:
$$
g_i (\boldsymbol{x} + \boldsymbol{c}_i \Delta t, t + \Delta t) = g_i(\boldsymbol{x}, t) + \Omega_{ij}^g(\bar{g_j}^{eq} - g_j)|_{(\boldsymbol{x}, t)} + F_i(\boldsymbol{x}, t).
$$
Here, the push scheme is applied which is easier due to the data access required for the viscous force term. Furthermore, the velocity update is directly done in the kernel.
%% Cell type:code id: tags:
``` python
force_Assignments = hydrodynamic_force_assignments(u, C, method_hydro, parameters, body_force)
lbm_optimisation = LBMOptimisation(symbolic_field=g, symbolic_temporary_field=g_tmp)
hydro_lb_update_rule = create_lb_update_rule(lbm_config=config_hydro,
lbm_optimisation=lbm_optimisation)
hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters,
config_hydro)
ast_kernel = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
kernel_hydro_lb = ast_kernel.compile()
```
%% Cell type:markdown id: tags:
## Boundary Conditions
%% Cell type:markdown id: tags:
As a last step suitable boundary conditions are applied
%% Cell type:code id: tags:
``` python
# periodic Boundarys for g, h and C
periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {"openmp": True})
periodic_BC_g = LBMPeriodicityHandling(stencil=stencil_hydro, data_handling=dh, pdf_field_name=g.name,
streaming_pattern=config_hydro.streaming_pattern)
periodic_BC_h = LBMPeriodicityHandling(stencil=stencil_phase, data_handling=dh, pdf_field_name=h.name,
streaming_pattern=config_phase.streaming_pattern)
# No slip boundary for the phasefield lbm
bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',
target=dh.default_target, name='boundary_handling_h',
streaming_pattern=config_phase.streaming_pattern)
# No slip boundary for the velocityfield lbm
bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,
target=dh.default_target, name='boundary_handling_g',
streaming_pattern=config_hydro.streaming_pattern)
contact_angle = BoundaryHandling(dh, C.name, stencil_hydro, target=dh.default_target)
contact = ContactAngle(90, parameters.interface_thickness)
wall = NoSlip()
if dimensions == 2:
bh_allen_cahn.set_boundary(wall, make_slice[:, 0])
bh_allen_cahn.set_boundary(wall, make_slice[:, -1])
bh_hydro.set_boundary(wall, make_slice[:, 0])
bh_hydro.set_boundary(wall, make_slice[:, -1])
contact_angle.set_boundary(contact, make_slice[:, 0])
contact_angle.set_boundary(contact, make_slice[:, -1])
else:
bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])
bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])
bh_hydro.set_boundary(wall, make_slice[:, 0, :])
bh_hydro.set_boundary(wall, make_slice[:, -1, :])
contact_angle.set_boundary(contact, make_slice[:, 0, :])
contact_angle.set_boundary(contact, make_slice[:, -1, :])
bh_allen_cahn.prepare()
bh_hydro.prepare()
contact_angle.prepare()
```
%% Cell type:markdown id: tags:
## Full timestep
%% Cell type:code id: tags:
``` python
# definition of the timestep for the immiscible fluids model
def timeloop():
# Solve the interface tracking LB step with boundary conditions
periodic_BC_h()
bh_allen_cahn()
dh.run_kernel(kernel_allen_cahn_lb, **parameters.symbolic_to_numeric_map)
# Solve the hydro LB step with boundary conditions
periodic_BC_g()
bh_hydro()
dh.run_kernel(kernel_hydro_lb, **parameters.symbolic_to_numeric_map)
dh.swap("C", "C_tmp")
# apply the three phase-phase contact angle
contact_angle()
# periodic BC of the phase-field
periodic_BC_C()
# field swaps
dh.swap("h", "h_tmp")
dh.swap("g", "g_tmp")
```
%% Cell type:code id: tags:
``` python
Initialize_distributions()
frames = 300
steps_per_frame = (timesteps//frames) + 1
if 'is_test_run' not in globals():
def run():
for i in range(steps_per_frame):
timeloop()
if gpu:
dh.to_cpu("C")
return dh.gather_array(C.name)
animation = plt.scalar_field_animation(run, frames=frames, rescale=True)
set_display_mode('video')
res = display_animation(animation)
else:
timeloop()
res = None
res
```
%% Output
<IPython.core.display.HTML object>
%% Cell type:markdown id: tags:
Note that the video is played for 10 seconds while the simulation time is only 2 seconds!
doc/notebooks/11_tutorial_Non_Newtonian_Flow.ipynb 0 → 100644
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/12_Thermocapillary_flows_heated_channel.ipynb 0 → 100644
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.
doc/notebooks/13_Thermocapillary_flows_droplet_motion.ipynb 0 → 100644
View file @ 32380aab
This source diff could not be displayed because it is too large. You can view the blob instead.