Skip to content
Snippets Groups Projects
Commit 02d87427 authored by Sebastian Bindgen's avatar Sebastian Bindgen
Browse files

Remove display of initial state to avoid runtime error

parent 2b542d9c
No related branches found
No related tags found
1 merge request!49Lees Edwards boundary conditions
Show whitespace changes
Inline Side-by-side
doc/notebooks/09_tutorial_lees-edwards.ipynb
+ 7
− 77
View file @ 02d87427
%% Cell type:markdown id: tags:
# Protoype of Lees Edwards 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
%% 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
from lbmpy.methods.creationfunctions import create_from_equilibrium
from lbmpy.relaxationrates import lattice_viscosity_from_relaxation_rate
from pystencils.astnodes import LoopOverCoordinate
from pystencils.slicing import get_slice_before_ghost_layer, get_ghost_region_slice
from pystencils.slicing import get_periodic_boundary_functor
from functools import partial
```
%% Cell type:markdown id: tags:
## Parameters
%% Cell type:code id: tags:
``` python
N = 64 # domain size
omega = 1. # relaxation rate of first component
U_x = 0.1 # shear velocity
shear_dir = 0 # direction of shear flow
shear_dir_normal = 1 # direction normal to the shear plane. Needed for the interpolation.
stencil = get_stencil("D2Q9")
weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))
```
%% 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.
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 = dh.add_array('src', values_per_cell=len(stencil))
dst = dh.add_array_like('dst', 'src')
F = dh.add_array('F', values_per_cell=dh.dim)
ρ = dh.add_array('rho')
u = dh.add_array('u', values_per_cell=dh.dim)
```
%% 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.
%% 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
```
%% 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.
%% 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]),
density_input=ρ,
force_model='luo',
force = F.center_vector,
kernel_type='collide_only',
optimization={'symbolic_field': src})
```
%% Output
/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/__init__.py:676: SymPyDeprecationWarning:
importing sympy.core.numbers with 'from sympy import *' has been
deprecated since SymPy 1.6. Use import sympy.core.numbers instead. See
https://github.com/sympy/sympy/issues/18245 for more info.
deprecated_since_version="1.6").warn()
%% Cell type:markdown id: tags:
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] 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):
if c[shear_dir_normal] == -1:
b = (True, False)
elif c[shear_dir_normal] == 1:
b = (False, True)
else:
b = (False, False)
a = ps.Assignment(a.lhs, a.rhs.replace(points_down, b[0]))
a = ps.Assignment(a.lhs, a.rhs.replace(points_up, b[1]))
ma.append(a)
collision.main_assignments = ma
```
%% Cell type:code id: tags:
``` 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()
```
%% Output
/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/__init__.py:676: SymPyDeprecationWarning:
importing sympy.logic.boolalg with 'from sympy import *' has been
deprecated since SymPy 1.6. Use import sympy.logic.boolalg instead.
See https://github.com/sympy/sympy/issues/18245 for more info.
deprecated_since_version="1.6").warn()
%% Cell type:markdown id: tags:
## Initialization
%% Cell type:code id: tags:
``` python
init = macroscopic_values_setter(collision.method, velocity=(0, 0),
pdfs=src.center_vector, density=ρ.center)
init_kernel = ps.create_kernel(init, ghost_layers=0).compile()
```
%% Output
/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/sympy/__init__.py:676: SymPyDeprecationWarning:
importing sympy.logic.boolalg with 'from sympy import *' has been
deprecated since SymPy 1.6. Use import sympy.logic.boolalg instead.
See https://github.com/sympy/sympy/issues/18245 for more info.
deprecated_since_version="1.6").warn()
%% Cell type:code id: tags:
``` python
def init():
dh.fill(ρ.name, 1.)
dh.run_kernel(init_kernel)
dh.fill(u.name, 0.0)
dh.fill(F.name, 0.0)
```
%% Cell type:markdown id: tags:
## 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.
%% Cell type:code id: tags:
``` python
def get_le_boundary_functor(stencil, offset, ghost_layers=1, thickness=None):
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, :]
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, :]
pdfs[i,ghost_layers-1,:] = res
return functor
```
%% Cell type:markdown id: tags:
## Timeloop
%% Cell type:code id: tags:
``` python
offset = [0.0]
sync_pdfs = dh.synchronization_function([src.name], functor = partial(get_le_boundary_functor, offset = offset))
def time_loop(steps):
dh.all_to_gpu()
for i in range(steps):
dh.run_kernel(collision_kernel)
sync_pdfs()
dh.run_kernel(stream_kernel)
dh.swap(src.name, dst.name)
offset[0] += U_x
dh.all_to_cpu()
```
%% Cell type:code id: tags:
``` python
def plot_v():
plt.subplot(121)
plt.title("$v_A$")
plt.vector_field(dh.gather_array(u.name), step=2)
plt.subplot(122)
plt.vector_field_magnitude(dh.gather_array(u.name))
plt.colorbar()
```
%% Cell type:markdown id: tags:
## Run the simulation
### Initial state
### Initialize all velocities with zero
%% Cell type:code id: tags:
``` python
init()
plot_v()
```
%% Output
/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/matplotlib/quiver.py:686: RuntimeWarning: divide by zero encountered in double_scalars
length = a * (widthu_per_lenu / (self.scale * self.width))
/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/matplotlib/quiver.py:686: RuntimeWarning: invalid value encountered in multiply
length = a * (widthu_per_lenu / (self.scale * self.width))
/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/matplotlib/quiver.py:738: RuntimeWarning: invalid value encountered in less
short = np.repeat(length < minsh, 8, axis=1)
/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/matplotlib/quiver.py:751: RuntimeWarning: invalid value encountered in less
tooshort = length < self.minlength
%% Cell type:markdown id: tags:
### Run the simulation until converged
### Run the simulation to show the flow profile
%% Cell type:code id: tags:
``` python
time = 500
time_loop(time)
plot_v()
```
%% Output
%% Cell type:markdown id: tags:
## Checking the results
This transient shear flow (planar Couette flow) can be modeled with the Navier Stokes equation. We use this as a testcase here.
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.
%% 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
```
%% Cell type:code id: tags:
``` python
nu = lattice_viscosity_from_relaxation_rate(omega)
h = N
k_max=100
X = np.linspace(0.5, h-0.5, h)
U = w(X, time, nu, U_x, h, k_max)
plt.plot(U, X, label='analytical')
X_lb = np.linspace(0.5, h-0.5, h)
plt.plot(dh.gather_array(u.name)[0,:,0], X_lb, 'o-', label='lbmpy')
plt.legend()
np.testing.assert_array_almost_equal(U, dh.gather_array(u.name)[0,:,0], decimal = 5)
```
%% Output
%% 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.
%% Cell type:code id: tags:
``` python
def time_loop_2(steps):
dh.all_to_gpu()
for i in range(steps):
dh.run_kernel(collision_kernel)
sync_pdfs()
dh.run_kernel(stream_kernel)
dh.swap(src.name, dst.name)
dh.all_to_cpu()
def init_2():
dh.fill(ρ.name, 1.)
dh.run_kernel(init_kernel)
dh.fill(u.name, 0.0)
dh.fill(F.name, 0.0)
dh.cpu_arrays[F.name][N//3,1,:] = [1e-2,-1e-1]
```
%% Cell type:markdown id: tags:
## Run the simulation without any offset and a constant offset
%% Cell type:code id: tags:
``` python
offset[0] = 0
init_2()
time = 20
time_loop_2(time)
plot_v()
vel_unshifted = np.array(dh.gather_array(u.name)[:,-3:-1,:])
```
%% Output
%% Cell type:code id: tags:
``` python
offset[0] = 10
init_2()
time = 20
time_loop_2(time)
plot_v()
vel_shifted = np.array(dh.gather_array(u.name)[:,-3:-1,:])
```
%% Output
%% Cell type:markdown id: tags:
### We role the array back by the offset and compare it to the results of the unshifted measurement
%% Cell type:code id: tags:
``` python
vel_rolled = np.roll(vel_shifted, -offset[0], axis=0)
np.testing.assert_array_almost_equal(vel_unshifted, vel_rolled)
```
%% Cell type:markdown id: tags:
### Finally let's run the simualtion with fixed offset a while longer
%% Cell type:code id: tags:
``` python
init_2()
time = 500
time_loop_2(time)
plot_v()
```
%% Output
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment