Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
pycodegen
lbmpy
Commits
78ab6911
Commit
78ab6911
authored
Dec 20, 2020
by
Sebastian Bindgen
Browse files
PEP8 notebook
parent
02d87427
Changes
1
Hide whitespace changes
Inline
Side-by-side
doc/notebooks/09_tutorial_lees-edwards.ipynb
View file @
78ab6911
...
...
@@ -7,17 +7,19 @@
%% 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.macroscopic_value_kernels
import
macroscopic_values_getter
from
lbmpy.macroscopic_value_kernels
import
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
pystencils.slicing
import
get_slice_before_ghost_layer
from
pystencils.slicing
import
get_ghost_region_slice
from
pystencils.slicing
import
get_periodic_boundary_functor
from
functools
import
partial
```
%% Cell type:markdown id: tags:
...
...
@@ -25,17 +27,17 @@
%% Cell type:code id: tags:
```
python
N
=
64
# domain size
omega
=
1.
# relaxation rate of first component
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
.
shear_dir_normal
=
1
# direction normal to shear plane
, for
interpolation
stencil
=
get_stencil
(
"D2Q9"
)
weights
=
get_weights
(
stencil
,
c_s_sq
=
sp
.
Rational
(
1
,
3
))
weights
=
get_weights
(
stencil
,
c_s_sq
=
sp
.
Rational
(
1
,
3
))
```
%% Cell type:markdown id: tags:
## Data structures
...
...
@@ -48,11 +50,13 @@
%% Cell type:code id: tags:
```
python
dim
=
len
(
stencil
[
0
])
dh
=
ps
.
create_data_handling
((
N
,
)
*
dim
,
periodicity
=
True
,
default_target
=
'cpu'
)
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
)
...
...
@@ -68,44 +72,56 @@
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
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
})
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
})
```
%% 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
]
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
:
...
...
@@ -121,12 +137,12 @@
```
%% Cell type:code id: tags:
```
python
stream
=
create_stream_pull_with_output_kernel
(
collision
.
method
,
src
,
dst
,
{
'density'
:
ρ
,
'velocity'
:
u
})
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
()
...
...
@@ -137,12 +153,12 @@
## Initialization
%% Cell type:code id: tags:
```
python
init
=
macroscopic_values_setter
(
collision
.
method
,
velocity
=
(
0
,
0
),
pdfs
=
src
.
center_vector
,
density
=
ρ
.
center
)
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
()
```
%% Cell type:code id: tags:
...
...
@@ -163,40 +179,40 @@
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
):
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
,:])):
for
i
in
range
(
len
(
pdfs
[:,
ghost_layers
,
:])):
ind1
=
int
(
np
.
floor
(
i
-
offset
[
0
])
%
N
)
ind2
=
int
(
np
.
ceil
(
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
,:])):
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
)
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
weight
*
pdfs
[
ind2
,
-
ghost_layers
-
1
,
:]
pdfs
[
i
,
ghost_layers
-
1
,
:]
=
res
return
functor
```
%% Cell type:markdown id: tags:
...
...
@@ -205,20 +221,24 @@
%% Cell type:code id: tags:
```
python
offset
=
[
0.0
]
sync_pdfs
=
dh
.
synchronization_function
([
src
.
name
],
functor
=
partial
(
get_le_boundary_functor
,
offset
=
offset
))
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
()
```
...
...
@@ -276,33 +296,37 @@
%% 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
nu
=
lattice_viscosity_from_relaxation_rate
(
omega
)
h
=
N
k_max
=
100
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
.
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
)
np
.
testing
.
assert_array_almost_equal
(
U
,
dh
.
gather_array
(
u
.
name
)[
0
,
:,
0
],
decimal
=
5
)
```
%%%% Output: display_data
...
...
@@ -318,24 +342,24 @@
```
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
]
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
...
...
@@ -348,11 +372,11 @@
init_2
()
time
=
20
time_loop_2
(
time
)
plot_v
()
vel_unshifted
=
np
.
array
(
dh
.
gather_array
(
u
.
name
)[:,
-
3
:
-
1
,:])
vel_unshifted
=
np
.
array
(
dh
.
gather_array
(
u
.
name
)[:,
-
3
:
-
1
,
:])
```
%%%% Output: display_data
...
...
@@ -365,11 +389,11 @@
init_2
()
time
=
20
time_loop_2
(
time
)
plot_v
()
vel_shifted
=
np
.
array
(
dh
.
gather_array
(
u
.
name
)[:,
-
3
:
-
1
,:])
vel_shifted
=
np
.
array
(
dh
.
gather_array
(
u
.
name
)[:,
-
3
:
-
1
,
:])
```
%%%% Output: display_data
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment
