Skip to content

GitLab

Commit e404ad34 authored by Sebastian Bindgen's avatar Sebastian Bindgen
Browse files

Corrections of Markus Holzer.

parent c42b6883
Hide whitespace changes
Inline Side-by-side
doc/notebooks/09_tutorial_lees-edwards.ipynb
View file @ e404ad34
%% Cell type:markdown id: tags:
# Protoype of Lees Edwards with lbmpy
# Lees Edwards boundary conditions with lbmpy
This examples shows how to implement Lees Edwards boundary conditions folowing the pricniples discussed in: Wagner, A.J., Pagonabarraga, I. Lees–Edwards Boundary Conditions for Lattice Boltzmann. Journal of Statistical Physics 107, 521–537 (2002). https://doi.org/10.1023/A:1014595628808
This example shows how to implement Lees Edwards boundary conditions following the principles discussed in Wagner, A.J., Pagonabarraga, I. Lees–Edwards Boundary Conditions for Lattice Boltzmann. Journal of Statistical Physics 107, 521–537 (2002). https://doi.org/10.1023/A:1014595628808
%% Cell type:code id: tags:
``` python
from lbmpy.session import *
......@@ -40,13 +40,12 @@
%% Cell type:markdown id: tags:
## Data structures
We allocate a set of PDFs.
For later testing we need a force field. This will be allcated as well.
We allocate a set of PDFs `src` and `dst`, the density field `ρ` and the velocity field `u`.
For later testing, we also need a force field `F`. This will be allocated as well.
To run the simulation on GPU, change the `default_target` to gpu
%% Cell type:code id: tags:
......@@ -67,60 +66,51 @@
%% Cell type:markdown id: tags:
## Kernels
Following Wagner et al., we need to find all the populations that will cross the boundary in the direction normal to the shearing direction and alter their equilibrium distribution.
Hence, we construct a piec wise function that fulfills this.
Hence, we construct a piecewise function that fulfils this.
%% Cell type:code id: tags:
``` python
counters = [LoopOverCoordinate.get_loop_counter_symbol(i) for i in
range(dim)]
(points_up, points_down) = (sp.Symbol('points_up'),
sp.Symbol('points_down'))
U_p = sp.Piecewise((1,
sp.And(ps.data_types.type_all_numbers(
counters[1] <= 1, 'int'),
points_down)),
(-1,
sp.And(ps.data_types.type_all_numbers(
counters[1] >= src.shape[1] - 2, 'int'),
points_up)),
(0, True)) \
* U_x
counters = [LoopOverCoordinate.get_loop_counter_symbol(i) for i in range(dim)]
points_up = sp.Symbol('points_up')
points_down = sp.Symbol('points_down')
U_p = sp.Piecewise((1, sp.And(ps.data_types.type_all_numbers(counters[1] <= 1, 'int'), points_down)),
(-1, sp.And(ps.data_types.type_all_numbers(counters[1] >= src.shape[1] - 2, 'int'),
points_up)), (0, True)) * U_x
```
%% Cell type:markdown id: tags:
For the LB update we will use a velocity input in the shearing direction with the magnitude U_x that is defined further up.
For the LB update, we will use a velocity input in the shearing direction with the magnitude `U_x` that is defined further up.
%% Cell type:code id: tags:
``` python
collision = create_lb_update_rule(stencil=stencil,
relaxation_rate=omega,
compressible=True,
velocity_input=u.center_vector+sp.Matrix(
[U_p, 0]),
velocity_input=u.center_vector+sp.Matrix([U_p, 0]),
density_input=ρ,
force_model='luo',
force=F.center_vector,
kernel_type='collide_only',
optimization={'symbolic_field': src})
```
%% Cell type:markdown id: tags:
We need to get the populations that cross the upper and lower boundary respectively.
We need to get the populations that cross the upper and lower boundary, respectively.
%% Cell type:code id: tags:
``` python
to_insert = [s.lhs for s in collision.subexpressions
if collision.method.first_order_equilibrium_moment_symbols[
shear_dir]
if collision.method.first_order_equilibrium_moment_symbols[shear_dir]
in s.free_symbols]
for s in to_insert:
collision = collision.new_with_inserted_subexpression(s)
ma = []
for a, c in zip(collision.main_assignments, collision.method.stencil):
......@@ -140,14 +130,12 @@
``` python
stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
{'density': ρ, 'velocity': u})
opts = {
'target': dh.default_target}
stream_kernel = ps.create_kernel(stream, **opts).compile()
collision_kernel = ps.create_kernel(collision, **opts).compile()
stream_kernel = ps.create_kernel(stream, target=dh.default_target).compile()
collision_kernel = ps.create_kernel(collision, target=dh.default_target).compile()
```
%% Cell type:markdown id: tags:
## Initialization
......@@ -174,11 +162,11 @@
## Interpolation
%% Cell type:markdown id: tags:
After applying the regular periodic boundary conditions we interpolate back in the original cells by using a linear interpolation scheme. In this step the corners are not special anymore so that we can just use the entire upper and lower slab.
After applying the normal periodic boundary conditions, we interpolate back in the original cells by using a linear interpolation scheme. In this step, the corners are not special anymore so that we can use the entire upper and lower slab.
%% Cell type:code id: tags:
``` python
def get_le_boundary_functor(stencil, offset, ghost_layers=1, thickness=None):
......@@ -186,31 +174,30 @@
functor_2 = get_periodic_boundary_functor(stencil, ghost_layers, thickness)
def functor(pdfs, **_):
functor_2(pdfs)
weight = np.fmod(offset[0] + N, 1.)
# First, we interpolate the upper pdfs
for i in range(len(pdfs[:, ghost_layers, :])):
ind1 = int(np.floor(i - offset[0]) % N)
ind2 = int(np.ceil(i - offset[0]) % N)
res = (1-weight) * pdfs[ind1, ghost_layers, :] + \
weight * pdfs[ind2, ghost_layers, :]
weight * pdfs[ind2, ghost_layers, :]
pdfs[i, -ghost_layers, :] = res
# Second, we interpolate the lower pdfs
for i in range(len(pdfs[:, -ghost_layers, :])):
ind1 = int(np.floor(i + offset[0]) % N)
ind2 = int(np.ceil(i + offset[0]) % N)
res = (1-weight) * pdfs[ind1, -ghost_layers-1, :] + \
weight * pdfs[ind2, -ghost_layers-1, :]
weight * pdfs[ind2, -ghost_layers-1, :]
pdfs[i, ghost_layers-1, :] = res
return functor
```
......@@ -222,13 +209,11 @@
``` python
offset = [0.0]
sync_pdfs = dh.synchronization_function([src.name],
functor=partial(
get_le_boundary_functor,
offset=offset))
functor=partial(get_le_boundary_functor, offset=offset))
def time_loop(steps):
dh.all_to_gpu()
for i in range(steps):
......@@ -290,22 +275,22 @@
The solution of the Navier Stokes equation in this particular case yields:
$ w(x,t) = v \left[ \frac{x}{h} - \frac{1}{2} + \sum_{k=1} ^\infty \frac{1}{\pi k} e^{-\frac{4 \pi^2 \nu k^2}{h^2} t} \sin \left( \frac{2 \pi}{h} kx \right) \right] $
with $u$ as the resulting velocity, $v$ as shear velocity, $h$ as the height, $\nu$ as kinematic viscosity and $t$ as time.
with $w$ as the resulting velocity, $v$ as shear velocity, $h$ as the height, $\nu$ as kinematic viscosity and $t$ as time.
%% Cell type:code id: tags:
``` python
def w(x, t, nu, v=1.0, h=1.0, k_max=10):
w = x/h-0.5
for k in np.arange(1, k_max+1):
w += 1.0/(np.pi*k) * \
np.exp(-4*np.pi**2*nu*k**2/h**2*t) * \
np.sin(2*np.pi/h*k*x)
return v*w
w = x / h - 0.5
for k in np.arange(1, k_max + 1):
w += 1.0 / (np.pi * k) * \
np.exp(-4 * np.pi**2 * nu * k**2 / h**2 * t) * \
np.sin(2 * np.pi / h * k * x)
return v * w
```
%% Cell type:code id: tags:
``` python
......@@ -333,11 +318,11 @@
%% Cell type:markdown id: tags:
## Testing of the interpolation scheme
We redefine our init function here. In order to test the interpolation scheme we send an impuls across the boundary in the shear gradient direction. If interpolated correctly it should appear shifted to the right by the preselected offset. Therefor we also need to fix the offset and redefine the time loop.
We redefine our init function here. In order to test the interpolation scheme, we send an impulse across the boundary in the shear gradient direction. If interpolated correctly, it should appear shifted to the right by the preselected offset. Therefore we also need to fix the offset and redefine the time loop.
%% Cell type:code id: tags:
``` python
def time_loop_2(steps):
......
Markdown is supported
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