Skip to content
GitLab
Explore
Sign in
Primary navigation
Search or go to…
Project
lbmpy
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Model registry
Operate
Environments
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Terms and privacy
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
pycodegen
lbmpy
Commits
ef7c7239
Commit
ef7c7239
authored
5 years ago
by
Martin Bauer
Browse files
Options
Downloads
Patches
Plain Diff
Re-added compiled-in boundaries with tests
are much faster in certain setups
parent
d4ac8672
No related branches found
No related tags found
No related merge requests found
Pipeline
#20454
passed
5 years ago
Stage: test
Stage: deploy
Changes
2
Pipelines
2
Hide whitespace changes
Inline
Side-by-side
Showing
2 changed files
lbmpy/boundaries/boundaries_in_kernel.py
+134
-0
134 additions, 0 deletions
lbmpy/boundaries/boundaries_in_kernel.py
lbmpy_tests/test_compiled_in_boundaries.ipynb
+175
-0
175 additions, 0 deletions
lbmpy_tests/test_compiled_in_boundaries.ipynb
with
309 additions
and
0 deletions
lbmpy/boundaries/boundaries_in_kernel.py
0 → 100644
+
134
−
0
View file @
ef7c7239
import
sympy
as
sp
from
lbmpy.boundaries.boundaryhandling
import
BoundaryOffsetInfo
,
LbmWeightInfo
from
pystencils.assignment
import
Assignment
from
pystencils.astnodes
import
Block
,
Conditional
,
LoopOverCoordinate
,
SympyAssignment
from
pystencils.data_types
import
type_all_numbers
from
pystencils.field
import
Field
from
pystencils.simp.assignment_collection
import
AssignmentCollection
from
pystencils.simp.simplifications
import
sympy_cse_on_assignment_list
from
pystencils.stencil
import
inverse_direction
from
pystencils.sympyextensions
import
fast_subs
def
direction_indices_in_direction
(
direction
,
stencil
):
for
i
,
offset
in
enumerate
(
stencil
):
for
d_i
,
o_i
in
zip
(
direction
,
offset
):
if
(
d_i
==
1
and
o_i
==
1
)
or
(
d_i
==
-
1
and
o_i
==
-
1
):
yield
i
break
def
boundary_substitutions
(
lb_method
):
stencil
=
lb_method
.
stencil
w
=
lb_method
.
weights
replacements
=
{}
for
idx
,
offset
in
enumerate
(
stencil
):
symbolic_offset
=
BoundaryOffsetInfo
.
offset_from_dir
(
idx
,
dim
=
lb_method
.
dim
)
for
sym
,
value
in
zip
(
symbolic_offset
,
offset
):
replacements
[
sym
]
=
value
replacements
[
BoundaryOffsetInfo
.
inv_dir
(
idx
)]
=
stencil
.
index
(
inverse_direction
(
offset
))
replacements
[
LbmWeightInfo
.
weight_of_direction
(
idx
)]
=
w
[
idx
]
return
replacements
def
border_conditions
(
direction
,
field
,
ghost_layers
=
1
):
abs_direction
=
tuple
(
-
e
if
e
<
0
else
e
for
e
in
direction
)
assert
sum
(
abs_direction
)
==
1
idx
=
abs_direction
.
index
(
1
)
val
=
direction
[
idx
]
loop_ctrs
=
[
LoopOverCoordinate
.
get_loop_counter_symbol
(
i
)
for
i
in
range
(
len
(
direction
))]
loop_ctr
=
loop_ctrs
[
idx
]
gl
=
ghost_layers
border_condition
=
sp
.
Eq
(
loop_ctr
,
gl
if
val
<
0
else
field
.
shape
[
idx
]
-
gl
-
1
)
if
ghost_layers
==
0
:
return
type_all_numbers
(
border_condition
,
loop_ctr
.
dtype
)
else
:
other_min
=
[
sp
.
Ge
(
c
,
gl
)
for
c
in
loop_ctrs
if
c
!=
loop_ctr
]
other_max
=
[
sp
.
Lt
(
c
,
field
.
shape
[
i
]
-
gl
)
for
i
,
c
in
enumerate
(
loop_ctrs
)
if
c
!=
loop_ctr
]
result
=
sp
.
And
(
border_condition
,
*
other_min
,
*
other_max
)
return
type_all_numbers
(
result
,
loop_ctr
.
dtype
)
def
transformed_boundary_rule
(
boundary
,
accessor_func
,
field
,
direction_symbol
,
lb_method
,
**
kwargs
):
tmp_field
=
field
.
new_field_with_different_name
(
"
t
"
)
rule
=
boundary
(
tmp_field
,
direction_symbol
,
lb_method
,
**
kwargs
)
bsubs
=
boundary_substitutions
(
lb_method
)
rule
=
[
a
.
subs
(
bsubs
)
for
a
in
rule
]
accessor_writes
=
accessor_func
(
tmp_field
,
lb_method
.
stencil
)
to_replace
=
set
()
for
assignment
in
rule
:
to_replace
.
update
({
fa
for
fa
in
assignment
.
rhs
.
atoms
(
Field
.
Access
)
if
fa
.
field
==
tmp_field
})
def
compute_replacement
(
fa
):
f
=
fa
.
index
[
0
]
shift
=
accessor_writes
[
f
].
offsets
new_index
=
tuple
(
a
+
b
for
a
,
b
in
zip
(
fa
.
offsets
,
shift
))
return
field
[
new_index
](
accessor_writes
[
f
].
index
[
0
])
substitutions
=
{
fa
:
compute_replacement
(
fa
)
for
fa
in
to_replace
}
all_assignments
=
[
assignment
.
subs
(
substitutions
)
for
assignment
in
rule
]
main_assignments
=
[
a
for
a
in
all_assignments
if
isinstance
(
a
.
lhs
,
Field
.
Access
)]
sub_expressions
=
[
a
for
a
in
all_assignments
if
not
isinstance
(
a
.
lhs
,
Field
.
Access
)]
assert
len
(
main_assignments
)
==
1
ac
=
AssignmentCollection
(
main_assignments
,
sub_expressions
).
new_without_subexpressions
()
return
ac
.
main_assignments
[
0
].
rhs
def
boundary_conditional
(
boundary
,
direction
,
read_of_next_accessor
,
lb_method
,
output_field
,
cse
=
False
):
stencil
=
lb_method
.
stencil
tmp_field
=
output_field
.
new_field_with_different_name
(
"
t
"
)
dir_indices
=
direction_indices_in_direction
(
direction
,
stencil
)
assignments
=
[]
for
direction_idx
in
dir_indices
:
rule
=
boundary
(
tmp_field
,
direction_idx
,
lb_method
,
index_field
=
None
)
boundary_subs
=
boundary_substitutions
(
lb_method
)
rule
=
[
a
.
subs
(
boundary_subs
)
for
a
in
rule
]
rhs_substitutions
=
{
tmp_field
(
i
):
sym
for
i
,
sym
in
enumerate
(
lb_method
.
post_collision_pdf_symbols
)}
offset
=
stencil
[
direction_idx
]
inv_offset
=
inverse_direction
(
offset
)
inv_idx
=
stencil
.
index
(
inv_offset
)
lhs_substitutions
=
{
tmp_field
[
offset
](
inv_idx
):
read_of_next_accessor
(
output_field
,
stencil
)[
inv_idx
]}
rule
=
[
Assignment
(
a
.
lhs
.
subs
(
lhs_substitutions
),
a
.
rhs
.
subs
(
rhs_substitutions
))
for
a
in
rule
]
ac
=
AssignmentCollection
([
rule
[
-
1
]],
rule
[:
-
1
]).
new_without_subexpressions
()
assignments
+=
ac
.
main_assignments
border_cond
=
border_conditions
(
direction
,
output_field
,
ghost_layers
=
1
)
if
cse
:
assignments
=
sympy_cse_on_assignment_list
(
assignments
)
assignments
=
[
SympyAssignment
(
a
.
lhs
,
a
.
rhs
)
for
a
in
assignments
]
return
Conditional
(
border_cond
,
Block
(
assignments
))
def
update_rule_with_push_boundaries
(
collision_rule
,
field
,
boundary_spec
,
accessor
,
read_of_next_accessor
):
method
=
collision_rule
.
method
loads
=
[
Assignment
(
a
,
b
)
for
a
,
b
in
zip
(
method
.
pre_collision_pdf_symbols
,
accessor
.
read
(
field
,
method
.
stencil
))]
stores
=
[
Assignment
(
a
,
b
)
for
a
,
b
in
zip
(
accessor
.
write
(
field
,
method
.
stencil
),
method
.
post_collision_pdf_symbols
)]
result
=
collision_rule
.
copy
()
result
.
subexpressions
=
loads
+
result
.
subexpressions
result
.
main_assignments
+=
stores
for
direction
,
boundary
in
boundary_spec
.
items
():
cond
=
boundary_conditional
(
boundary
,
direction
,
read_of_next_accessor
,
method
,
field
)
result
.
main_assignments
.
append
(
cond
)
if
'
split_groups
'
in
result
.
simplification_hints
:
substitutions
=
{
b
:
a
for
a
,
b
in
zip
(
accessor
.
write
(
field
,
method
.
stencil
),
method
.
post_collision_pdf_symbols
)}
new_split_groups
=
[]
for
split_group
in
result
.
simplification_hints
[
'
split_groups
'
]:
new_split_groups
.
append
([
fast_subs
(
e
,
substitutions
)
for
e
in
split_group
])
result
.
simplification_hints
[
'
split_groups
'
]
=
new_split_groups
return
result
This diff is collapsed.
Click to expand it.
lbmpy_tests/test_compiled_in_boundaries.ipynb
0 → 100644
+
175
−
0
View file @
ef7c7239
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from lbmpy.session import *\n",
"from lbmpy.boundaries.boundaries_in_kernel import update_rule_with_push_boundaries\n",
"from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter\n",
"from collections import OrderedDict\n",
"from time import perf_counter"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Version 1: compile-in boundaries"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"domain_size = (32, 32, 32)\n",
"relaxation_rate = 1.8\n",
"time_steps = 100\n",
"lid_velocity = 0.05"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"dh = create_data_handling(domain_size, default_target='cpu')\n",
"pdfs = dh.add_array('pdfs', values_per_cell=19)\n",
"u = dh.add_array('u', values_per_cell=len(domain_size))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"boundaries = OrderedDict((\n",
" ((0, 1, 0), UBB([lid_velocity, 0, 0])), \n",
" ((1, 0, 0), NoSlip()),\n",
" ((-1, 0, 0), NoSlip()),\n",
" ((0, -1, 0), NoSlip()),\n",
" ((0, 0, 1), NoSlip()),\n",
" ((0, 0, -1), NoSlip()),\n",
"))\n",
"opt = {'symbolic_field': pdfs, 'cse_global': False, 'cse_pdfs': True}\n",
"cr_even = create_lb_collision_rule(stencil=\"D3Q19\", relaxation_rate=relaxation_rate, compressible=False, optimization=opt)\n",
"cr_odd = create_lb_collision_rule(stencil=\"D3Q19\", relaxation_rate=relaxation_rate, compressible=False, optimization=opt)\n",
"update_rule_aa_even = update_rule_with_push_boundaries(cr_even, pdfs, boundaries, AAEvenTimeStepAccessor, AAOddTimeStepAccessor.read)\n",
"update_rule_aa_odd = update_rule_with_push_boundaries(cr_odd, pdfs, boundaries, AAOddTimeStepAccessor, AAEvenTimeStepAccessor.read)\n",
"\n",
"getter_assignments = macroscopic_values_getter(update_rule_aa_even.method, velocity=u.center_vector,\n",
" pdfs=pdfs.center_vector, density=None)\n",
"\n",
"getter_kernel = ps.create_kernel(getter_assignments, target=dh.default_target).compile()\n",
"even_kernel = ps.create_kernel(update_rule_aa_even, target=dh.default_target, ghost_layers=1).compile()\n",
"odd_kernel = ps.create_kernel(update_rule_aa_odd, target=dh.default_target, ghost_layers=1).compile()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def init():\n",
" dh.fill(pdfs.name, 0, ghost_layers=True)\n",
"\n",
"def aa_time_loop(steps=100):\n",
" assert steps % 2 == 0, \"Works only for an even number of time steps\"\n",
" dh.all_to_gpu()\n",
" for i in range(steps // 2):\n",
" dh.run_kernel(odd_kernel)\n",
" dh.run_kernel(even_kernel)\n",
" dh.run_kernel(getter_kernel) \n",
" dh.all_to_cpu()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0wAAAFlCAYAAADYoWhgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3df6xc5X3n8c9n7vUPfkOwoQRDTIo3CqEbolguWlopCYU43bQmW1BMsw3SIrm/kJJuuyvSCtqiRCpSG6pu2Uik0DhsKUSkbG5bJ4QEIjZV1vGlkIBxUG4oLRd7sR3zm/jHnfnuH3OcjC8z937PuTN3Zu68X9LRnTnzned5Zs6Zc+c7z3Oe44gQAAAAAOCNav1uAAAAAAAMKhImAAAAAOiAhAkAAAAAOiBhAgAAAIAOSJgAAAAAoAMSJgAAAADoYHwxK1u1alWsXbt2MasEBs4Tu59PxS076Ugq7szlL6frrid/IxlTIxW3wjPJ8vKXL8hG1uxk3bm4XuhXze5BzZHcMt2uO1tvL2RrbqTfm+47lGxkPXKf/WXOffYlac+Rk1Nxr710XCruwjefma4bWIoeeeSR/RGxut/tKOP97z0hfnigXvn5j3z30P0RsbGLTeqJRU2Y1q5dq8nJycWsEhg4b/vjW1Jxb/756VTcx97ytXTdrzRyX1xOqB1Kxb1t2d5kefkvYfXkF8ATarmvnyd6WbrujLFkoiZJtWSCWuvyV+kxd3/wQD1y27DbdWfrlfKJS/b9nlHuS8ArjcOpuGUl3ptsoj91JBf3YvKz/+bxV1JxkvTJ3R9Ixe34yoWpuMk//J103cBSZPtf+92GsvYfqGv7/WsqP3/ZWT9Y1cXm9AxD8gAAAACgg0XtYQIAAACwVESpkQDDioQJAAAAQGmh/HDoYUbCBAAAAKCSRnKiqGHGOUwAAAAA0AE9TAAAAABKC4XqwZA8AAAAAGiLc5gAdF2M5Q4s0//37FTcX3xqc7ruxljumi0HT89du+hHp+dG9c4cnwqTJDWSl01Kx63Ivd8xlisvef3PIjZZd7bMPg6iji5fddX5q8L2oMzci6nlrsssZ6/ZWOI7xfjruTYuezVX3rLXcpUfV+IClHv+Q+5DEycv/fMbgFEVkuokTAAAAADQ3ij0MM37e6Xtlba/bfs7tnfa/uNi/Xm2t9v+vu17bC/vfXMBAAAAYPFkBngckvS+iHinpIskbbR9saSbJd0SEeskvSDp2t41EwAAAMAgCUn1iMrLsJg3YYqmo6OklxVLSHqfpHuL9VslXdGTFgIAAAAYSI0FLMMidQ6T7TFJj0g6X9Ktkn4g6cWIOHpK7LSk3BnqAAAAAIZeKJj04aiIqEu6yPapku6T9PZ2Ye2ea3uLpC2SdO6551ZsJgAAAICBElJ96edL5SapjYgXJX1D0sWSTrV9NOFaI2l3h+fcFhHrI2L96tWrF9JWAAAAAFhUmVnyVhc9S7J9nKRfkLRL0kOSrizCrpH0pV41EgAAAMBgCXEO01FnSdpanMdUk/SFiPgH209Kutv2JyU9Kun2HrYTAAAAwECx6urylc0H0LwJU0R8V9K72qx/WtKGXjQKWMre+r/2puJ2/bc3peL+5Zfyl0Bbd9er8wdJWvnsS6m4UxvJ34deztUrST5uZSouxsfSZaasyL2PjeU9uN73WO6fTdS6/E+p2+VJUiM3mN3JOLlEG2dy+2Pt0JHu1p2dGveHL+biJHk8t5/FySfkCqzlRuB/7zdPzZUn6bizXk7Fnf3pZbkCfyddNYABEUof9odaqXOYAAAAAGCU9OCnUgAAAACjgCF5AAAAANBGiIQJAAAAADpqBAkTAAAAALzBqPQwMekDAAAAAHRADxMAAACA0kJWfQT6X0iYAAAAAFTCOUwAAAAA0MaonMNEwgQAAACgAqseDMkD0G0RqbBV38p9PI+cmP9l58CFJ6XiTvh/x6Xilh84mIrzKcen4iSp9tqhXJmv5+rW4SOpsHjx5VRcLVmeJEW9ngzM7RPKltdIltcLtdz+aCf321r+H7GXL8sFjuc+W07G6biVqbBYdVquPEn15bm6D/5U7rP10ltz783yH6bCJEl+9pRU3PiBfflCAWAAkTABAAAAKC0kNUZg0oel/woBAAAA9ERdrrxk2N5o+ynbU7avb/P4Ctv3FI9vt722WL/B9mPF8h3bH2p5zjO2Hy8em5yvDfQwAQAAACgtorfnMNkek3SrpMskTUvaYXsiIp5sCbtW0gsRcb7tzZJulvRhSU9IWh8RM7bPkvQd238fETPF894bEfsz7aCHCQAAAMAg2iBpKiKejojDku6WtGlWzCZJW4vb90q61LYj4vWW5GilmiMIKyFhAgAAAFBJQ668JJwt6dmW+9PFurYxRYL0kqTTJcn2z9reKelxSb/RkkCFpK/afsT2lvkawZA8AAAAAKU1r8O0oP6XVbPOIbotIm5rud8uq5rdU9QxJiK2S3qH7bdL2mr7yxFxUNIlEbHb9hmSHrD9vYh4uFMjSZgAAAAAVLDgc5j2R8T6OR6flnROy/01knZ3iJm2PS7pFEkHWgMiYpft1yRdKGkyInYX6/favk/NoX8dEyaG5AEAAAAo7ei04lWXhB2S1tk+z/ZySZslTcyKmZB0TXH7SkkPRkQUzxmXJNtvkfQ2Sc/YPsH2ScX6EyRdruYEER3RwwQAAABg4BQz3F0n6X5JY5LuiIidtm9Ss6doQtLtku60PaVmz9Lm4uk/J+l620ckNST9VkTst/1WSfcVF1Afl3RXRHxlrnaQMAGLrHHycam4k/7tcCouxnLXMZCk135qWSru8Em5zufxV8dSca7nJ6ZpvOmEVFxt5fJc3OuH0nWnNBrp0NqRmfmDJCkqT9zTXqPL5UlSLb+fpThXXozn9jFJUq27gyYax69IxdWPz32uStW9PPe6f7Q692+8kfxvf9r38vv3+KHcflY/aWW6TADDpx5d/v8wS0Rsk7Rt1robW24flHRVm+fdKenONuuflvTOMm0gYQIAAABQWsgLnfRhKJAwAQAAAKik0cML1w4KEiYAAAAApXVhWvGhsPRfIQAAAABURA8TAAAAgNJC7vmkD4OAhAkAAABAJcnrKQ01EiYAAAAApUVI9RGY9GHpv0IAAAAAqIgeJgAAAAAVWA1xDhOALhvb+1IqzqeuTMWF8weqk6YPp2Mz6ivHUnHjr8+ky6wdPFK1OW3FymWpOB+pJwuMfN3LkofYZJluJOvObZaeiFpyf8zutyX2b43nXngsS8aN5equHWmk4urH5//lHjkxF7vi5dx+u+LldNVptcO5/XHswKvdrxzAQAiNxpA8EiYAAAAAlYzCdZhImAAAAACUFrIaIzCt+NJPCQEAAACgInqYAAAAAFTCkDwAAAAAaCMkNZj0AQAAAADasepMKw4AAAAAbzQqPUxL/xUCAAAAQEX0MAEAAACohCF5ALquceCFVNyyl09LxR0+ZXm67rHXZ1JxjkjFNcZzndSeaaTiJMlH6unYlGzV9WRg8r0pw41kmdm63YN/XunX3eW6y7zf2W2Y3W/rubqjkavXM/nXMv6j3Ocg28asIyfmvxYsf+lILnDvDyu2BsCgi/BIDMkjYQIAAABQSX0EEqZ5X6Htc2w/ZHuX7Z22P1as/yPbz9l+rFh+sffNBQAAAIDFk+lhmpH0uxHxz7ZPkvSI7QeKx26JiD/tXfMAAAAADKKQ1OAcJiki9kjaU9x+xfYuSWf3umEAAAAABpkZkjeb7bWS3iVpe7HqOtvftX2H7dwZ6gAAAACGXvM6TK68DIt0wmT7RElflPTxiHhZ0mck/bSki9TsgfqzDs/bYnvS9uS+ffu60GQAAAAAg6CuWuVlWKRaanuZmsnS30TE30lSRDwfEfWIaEj6rKQN7Z4bEbdFxPqIWL969eputRsAAAAAem7ec5hsW9LtknZFxKdb1p9VnN8kSR+S9ERvmggAAABg0ISGa2hdVZlZ8i6R9GuSHrf9WLHu9yVdbfsiNYcvPiPp13vSQgAAAAADqTFEQ+uqysyS9021v3T7tu43BwAAAMAwiJDq9DAB6DavXJGLe/G1VNyKQ0fSdTeOX56re6aRiht/9XCu4jI/Pjl34PXhmVx5M/USlSdEpENdz72P6TKT740ayXrLyG6X7GtOivGxEsG599EHk5+ZZbl/kbEyF1fL1iup9noutpGsO8Zy2++453+UipOk2qsHc4HJ9xHAcBqFIXlLvw8NAAAAACriZx8AAAAApTUnfVj6/S9L/xUCAAAA6Im6XHnJsL3R9lO2p2xf3+bxFbbvKR7fbnttsX6D7ceK5Tu2P5QtczZ6mAAAAACUFurtOUy2xyTdKukySdOSdtieiIgnW8KulfRCRJxve7OkmyV9WM1LHq2PiBnbZ0n6ju2/L5o9X5nHoIcJAAAAwCDaIGkqIp6OiMOS7pa0aVbMJklbi9v3SrrUtiPi9Yg4OkPUSjUTpWyZxyBhAgAAAFBB8xymqoukVbYnW5Ytsyo4W9KzLfeni3VtY4oE6SVJp0uS7Z+1vVPS45J+o3g8U+YxGJIHAAAAoJJG8lykDvZHxPo5Hm9X+OxrSHSMiYjtkt5h++2Sttr+crLMY5AwAQAAAChtES5cOy3pnJb7ayTt7hAzbXtc0imSDhzbzthl+zVJFybLPAZD8gAAAABUssAhefPZIWmd7fNsL5e0WdLErJgJSdcUt6+U9GBERPGccUmy/RZJb5P0TLLMY9DDBCy2xpy9vj/mw0dy5UWuPKn7v5C4Xs8FJsNKmckV6nqju/WWeL/V6GPd3danuj0zf8xRMdblPTy5j9UOdrdaSQrnfrFNv+JkYO3VQ9kS08eoPu61AIZcMcPddZLulzQm6Y6I2Gn7JkmTETEh6XZJd9qeUrNnaXPx9J+TdL3tI5Iakn4rIvZLUrsy52oHCRMAAACA0poXru3pkDxFxDZJ22atu7Hl9kFJV7V53p2S7syWORcSJgAAAACVLHDSh6FAwgQAAACgtF5fuHZQMOkDAAAAAHRADxMAAACASpKz3Q01EiYAAAAA5UXvJ30YBCRMAAAAAEoLMekDAAAAAHQ0Cj1MS3/QIQAAAABURA8TsMgar76Wihs7/rhcgcuXpeuuvXowFRe15K9FY2O5uHo9FyfJ9UYyMNnGiFxcI1lvI1leL2RfS/a96Wfd2X1spsS+k92Gyf02lKvbR2Zy5Y0nPy+SnIz167m60/t3mX0nuW3itdfzZQIYKqMyrTgJEwAAAIBKSJgAAAAAoI0Qs+QBAAAAQEejMEsekz4AAAAAQAf0MAEAAAAoLziHCQAAAADaYpY8AAAAAJjDKCRMnMMEAAAAAB3QwwQAAACgNKYVB9BXcehQKs61Eh3FtdxBzfVkeTPZwBIiulteI1let+sto591Z2XbmI7rwT/YRjYwt9+6kWyjk5+rMp+XbGy3t0sJcTB3jIp6esMAGEJBwgQAAAAA7Y3CdZhImAAAAACUFiMyrTiTPgAAAABAB/QwAQAAAKiEc5gAAAAAoC1myQMAAACAjuhhAgAAAIA2Qkz6AAAAAAAjjR4mAAAAAOXFcFx7faFImIDFVst1Xcehw12v2itWJAOT3ev1ei5ubCwXV6buRiNfZr90u4215KCA7HbpRd39es1lzCTbmN0Xx7LvTYlvFdlvINk2Zqs9eCgffCgZmzzmARhOXLgWAAAAANoIjcakD/P+LGb7HNsP2d5le6ftjxXr32T7AdvfL/6e1vvmAgAAAMDiyYwjmJH0uxHxdkkXS/pt2xdIul7S1yNinaSvF/cBAAAAjITmdZiqLsNi3oQpIvZExD8Xt1+RtEvS2ZI2SdpahG2VdEWvGgkAAABg8ERUXzJsb7T9lO0p22/ooLG9wvY9xePbba8t1l9m+xHbjxd/39fynG8UZT5WLGfM1YZS5zAVDXiXpO2SzoyIPVIzqZqvIgAAAABLSy/PYbI9JulWSZdJmpa0w/ZERDzZEnatpBci4nzbmyXdLOnDkvZL+qWI2G37Qkn3q9npc9RHImIy04709EO2T5T0RUkfj4iXSzxvi+1J25P79u3LPg0AAADAAGv2FLnykrBB0lREPB0RhyXdreYot1ato97ulXSpbUfEoxGxu1i/U9JK28npgo+VSphsL1MzWfqbiPi7YvXzts8qHj9L0t52z42I2yJifUSsX716dZU2AgAAABg9Z0t6tuX+tI7tJTomJiJmJL0k6fRZMb8i6dGIaL0ewl8Xw/FusOe+RkNmljxLul3Sroj4dMtDE5KuKW5fI+lL85UFAAAAYOlY4KQPq46ORCuWLbOKb5fIzD77ac4Y2+9Qc5jer7c8/pGI+BlJP18svzbXa8ycw3RJUcjjth8r1v2+pD+R9AXb10r6N0lXJcoCAAAAsERkJ2/oYH9ErJ/j8WlJ57TcXyNpd4eYadvjkk6RdECSbK+RdJ+kj0bED37S5niu+PuK7bvUHPr3+U6NmDdhiohvqn3mJkmXzvd8AAAAAEtTjy9cu0PSOtvnSXpO0mZJvzor5uiot29JulLSgxERtk+V9I+SPhER/3Q0uEiqTo2I/cVpRx+U9LW5GlFqljwAXVCv5+IajWR5yThJcfhILrCWO/jNM+T3J7KvuZ8W+BPZoujn+9ivusvUm90fa8n5jrL7xEyXP9M9ENm6G/n3O7LbZhg+/wAqCaUnb6hWfsSM7evUnOFuTNIdEbHT9k2SJiNiQs1Th+60PaVmz9Lm4unXSTpf0g22byjWXS7pNUn3F8nSmJrJ0mfnagcJEwAAAICBFBHbJG2bte7GltsH1ebUoIj4pKRPdij23WXaQMIEAAAAoJIhGJ+xYCRMAAAAAMqLnp/DNBBImAAAAABUMwJdTMkzXwEAAABg9NDDBAAAAKAShuQBAAAAQAfDcFWOhSJhAgAAAFBaiB4mAAAAAGgvJJEwAei2SF713sk4RSNf+cxMLm4sNx9MjI3lyutBf73d3QN0jMKYgiUuvUc0SnxmMvq470T2M51tY6PEa0mWmT3mAcCgImECAAAAUMko/N5IwgQAAACgGhImAAAAAGjHTPoAAAAAAB2NQA9T7sxuAAAAABhB9DABAAAAKC+4DhMAAAAAdDYCQ/JImAAAAABUtPR7mDiHCQAAAAA6oIcJGFBRb+QCs3GSPJ78jSR7Fbp6PVlxid9marlfqqJfV8prlKg3ktumzPuDjiL7Ng7DVRazn63sa0nutzEzkytPUmTbCGBpG4JD6kKRMAEAAACohoQJAAAAANoIScySBwAAAADtDcMo54Vi4DwAAAAAdEAPEwAAAIBqRqCHiYQJAAAAQDWcwwQAAAAA7ZkeJgAAAABoIzQSQ/KY9AEAAAAAOqCHCVhs7vLvFPV6PnYsWXcjGZd+KY1sYL7ubosSbUyXmfzZLZLb0Mlx4r2Y4zVbdz/1a9/JKrOPdXsbZusuczxpJNvY7WMegAFizmECAAAAgI5GYEgeCRMAAACAakYgYaKfHAAAAAA6oIcJAAAAQDX0MAEAAABAG6HmpA9VlwTbG20/ZXvK9vVtHl9h+57i8e221xbrL7P9iO3Hi7/va3nOu4v1U7b/wp57ZiMSJgAAAACVOKov85Ztj0m6VdIHJF0g6WrbF8wKu1bSCxFxvqRbJN1crN8v6Zci4mckXSPpzpbnfEbSFknrimXjXO0gYQIAAABQTSxgmd8GSVMR8XREHJZ0t6RNs2I2Sdpa3L5X0qW2HRGPRsTuYv1OSSuL3qizJJ0cEd+KiJD0eUlXzNUIEiYAAAAA/bDK9mTLsmXW42dLerbl/nSxrm1MRMxIeknS6bNifkXSoxFxqIifnqfMYzDpAwAAAIB+2B8R6+d4vN25RbP7puaMsf0ONYfpXV6izGPQwwQAAACgkl6ew6Rm7885LffXSNrdKcb2uKRTJB0o7q+RdJ+kj0bED1ri18xT5jHoYQIGVS03e0wZUW+k4jyWLLCR/M2lzE8zUc/FzT2hTXnRx3lRG9m6+zl3a7LuHuy3ebn9u296sY/Vk5+X9D5WQnZbJ5sIYEglZ7uraIekdbbPk/ScpM2SfnVWzISakzp8S9KVkh6MiLB9qqR/lPSJiPinHzc3Yo/tV2xfLGm7pI9K+h9zNWLerzG277C91/YTLev+yPZzth8rll/MvGIAAAAAS8RCJnxI/I5TnJN0naT7Je2S9IWI2Gn7Jtu/XITdLul021OS/quko1OPXyfpfEk3tOQsZxSP/aakv5I0JekHkr48VzsyPUyfk/SXas4g0eqWiPjTxPMBAAAAoLSI2CZp26x1N7bcPijpqjbP+6SkT3Yoc1LShdk2zNvDFBEPqxgHCAAAAAA/1ttpxQfCQiZ9uM72d4she6d1rUUAAAAAhkKPJ30YCFUTps9I+mlJF0naI+nPOgXa3nJ0bvV9+/ZVrA4AAADAwKGHqb2IeD4i6hHRkPRZNa/C2yn2tohYHxHrV69eXbWdAAAAALDoKiVMts9qufshSU90igUAAACwRI1AD9O8s+TZ/ltJ75G0yva0pD+U9B7bF6n5Up+R9Os9bCMAAACAATNs5yJVNW/CFBFXt1l9ew/aAgAAAGCY9PbCtQMhcx0mAKOmkfy5qNZIlldi9G8teeCNPv2klX1veiGS73c/ldnW3dbHqlN6se8ky4x+fV4ALH0jcHgZ9H8vAAAAANA39DABAAAAqIRzmAAAAACgExImAAAAAGhjRGbJ4xwmAAAAAOiAHiYAAAAA1YxADxMJEwAAAIBqSJgAAAAAoD3OYQIAAACAEUYPEzDsGo187NhYl+tO/qxUK9HGxhL6HSdKvO5uym6XMmrOxXX7NbvE/tCL191N/dofyihzPAGAEUHCBAAAAKCaAf+tqhtImAAAAACUNyLXYSJhAgAAAFANCRMAAAAAdDACCdMSOrsaAAAAALqLHiYAAAAApVmcwwQAAAAAnZEwAQAAAEAbIzJLHucwAQAAAEAH9DABo6TRSIVFrbu/pThXbVOtTPAS0RiCn+e63caac3FRYn9wl38DLFN3Rg+2c0SyzORnHwBKG4J/YQtFwgQAAACgGhImAAAAAGhvFM5hImECAAAAUM0IJExM+gAAAAAAHZAwAQAAACgvFrgk2N5o+ynbU7avb/P4Ctv3FI9vt722WH+67Ydsv2r7L2c95xtFmY8VyxlztYEheQAAAAAq6eU5TLbHJN0q6TJJ05J22J6IiCdbwq6V9EJEnG97s6SbJX1Y0kFJN0i6sFhm+0hETGbaQQ8TAAAAgGp628O0QdJURDwdEYcl3S1p06yYTZK2FrfvlXSpbUfEaxHxTTUTpwUhYQIAAABQiaP6ImmV7cmWZcus4s+W9GzL/eliXduYiJiR9JKk0xNN/+tiON4Ntue8OCBD8gAAAAD0w/6IWD/H4+0Smdl9U5mY2T4SEc/ZPknSFyX9mqTPdwomYQJGSSQHGjcaubharpM6svVKcrLqtNqcPxoNlTLvY7/M8yNdb0W3d56kRm679GT7ZT+rQ7DvABhSvT28TEs6p+X+Gkm7O8RM2x6XdIqkA3MVGhHPFX9fsX2XmkP/OiZMDMkDAAAAUF7vZ8nbIWmd7fNsL5e0WdLErJgJSdcUt6+U9GDM8QuV7XHbq4rbyyR9UNITczWCHiYAAAAApVntx8N1S0TM2L5O0v2SxiTdERE7bd8kaTIiJiTdLulO21Nq9ixt/nH77GcknSxpue0rJF0u6V8l3V8kS2OSvibps3O1g4QJAAAAwECKiG2Sts1ad2PL7YOSrurw3LUdin13mTaQMAEAAACoZgROkSRhAgAAAFBJLy9cOyhImAAAAABUQ8IEAAAAAB2MQMLEtOIAAAAA0AE9TAAAAADKC85hAgAAAIDOSJgA9E0jdwSKEgNr3WjkAmvJQuv1ZMX5y9pFsm5ny0y+j6r14NJ72bpHUT+3S1b2M9j5gvILqDv5Wc3WnSyv1Gth/wag0ehhmvebie07bO+1/UTLujfZfsD294u/p/W2mQAAAAAGTixgGRKZn3I/J2njrHXXS/p6RKyT9PXiPgAAAAAsKfMmTBHxsKQDs1ZvkrS1uL1V0hVdbhcAAACAAeeovgyLqucwnRkReyQpIvbYPqOLbQIAAAAw6IZsaF1VPb8Ok+0ttidtT+7bt6/X1QEAAABYLJzD1NHzts+SpOLv3k6BEXFbRKyPiPWrV6+uWB0AAAAALL6qCdOEpGuK29dI+lJ3mgMAAABgGFijcQ5TZlrxv5X0LUlvsz1t+1pJfyLpMtvfl3RZcR8AAADAKBmBIXnzTvoQEVd3eOjSLrcFAAAAwBBxLy7ePWCqzpIHYFA08geqSA7CdaORC6z1YN6YZN2RrNt2st7+HfCj2/9sstuvjOT7nX0t6e3SC8lt3dftkq07+3lJl7f0v/gA6KIh6ymqquez5AEAAADAsKKHCQAAAEAlwzR5Q1UkTAAAAACqIWECAAAAgPboYQIAAACATkYgYWLSBwAAAADogB4mAAAAAOUFQ/IAAAAAoDMSJgAAAAB4I4seJgC9EI1knHNxLnEqYiN3VItkkW4kX0utRBudfN0YDcl9tq+yn4Mo8VqSZUa2zF68j+ljWTIOAAYUCRMAAACAasr8GDSkSJgAAAAAVMKQPAAAAABoJ8SkDwAAAADQiUfgNEUuXAsAAABgINneaPsp21O2r2/z+Arb9xSPb7e9tlh/uu2HbL9q+y9nPefdth8vnvMX9twzTpEwAQAAAKgmFrDMw/aYpFslfUDSBZKutn3BrLBrJb0QEedLukXSzcX6g5JukPR7bYr+jKQtktYVy8a52kHCBAAAAKASR/UlYYOkqYh4OiIOS7pb0qZZMZskbS1u3yvpUtuOiNci4ptqJk4/aa99lqSTI+Jb0bw2w+clXTFXI0iYAAAAAJQXak4rXnWZ39mSnm25P12saxsTETOSXpJ0+jxlTs9T5jGY9AEAAABAJQucVnyV7cmW+7dFxG2txbd5zuwaMzELiSdhAoZelJiext3tVI7kxercKNHGWrKNyTIjWd4853tiCcnut9l9LH3RxhKfg3Qbu63M8QQAFm5/RKyf4/FpSee03F8jaXeHmGnb45JOkXRgnjLXzFPmMRiSB/u7CFoAAAx0SURBVAAAAKCaHk76IGmHpHW2z7O9XNJmSROzYiYkXVPcvlLSgzHHr04RsUfSK7YvLmbH+6ikL83VCHqYAAAAAJRmLXhI3pwiYsb2dZLulzQm6Y6I2Gn7JkmTETEh6XZJd9qeUrNnafOP22c/I+lkScttXyHp8oh4UtJvSvqcpOMkfblYOiJhAgAAAFBefvKGBVQR2yRtm7XuxpbbByVd1eG5azusn5R0YbYNDMkDAAAAgA7oYQIAAABQSS+H5A0KEiYAAAAA1ZAwAQAAAEB79DABAAAAQDshqbH0MyYmfQAAAACADuhhAgZU9OAXG9caucBG8reUmqs3ZsDMcY27YzSvcTfgavwWNjKyx4lIfvZL6MUxCsAQGoFDAQkTAAAAgEo4hwkAAAAAOunxhWsHAQkTAAAAgEpGoYeJge4AAAAA0AE9TAAAAADKCzHpAwAAAAC0Y0nmHCYAAAAA6KD7Vy0YOJzDBAAAAAAd0MMEAAAAoBKG5AEAAABAO0z6AADzaOSOklFi8K8bycHQtWShXS4vSvySZjsdi/bKvN9p2X0iW3eyvFKvJfnZ6rboU70AhlVw4dr52H5G0iuS6pJmImJ9NxoFAAAAYPCNwoVru9HD9N6I2N+FcgAAAABgoDAkDwAAAEA1IzAkb6HTioekr9p+xPaWbjQIAAAAwBAIyY3qy7BYaA/TJRGx2/YZkh6w/b2IeLg1oEiktkjSueeeu8DqAAAAAAwMepjmFhG7i797Jd0naUObmNsiYn1ErF+9evVCqgMAAACARVU5YbJ9gu2Tjt6WdLmkJ7rVMAAAAAADLhawDImFDMk7U9J9xXVGxiXdFRFf6UqrAAAAAAw8j8CQvMoJU0Q8LemdXWwLAAAAgGFCwgSg69IHluT0Mc6PrI1Grm7Xul/3wGskX3OtxPud3NZFTz1QXuT22+xnvxd1j8KXKWBkhdJfV4bZEvq2AwAAAADdRQ8TAAAAgNKs4BwmAAAAAOiIhAkAAAAAOiBhAgAAAIA2mPQBAAAAAPrH9kbbT9mesn19m8dX2L6neHy77bUtj32iWP+U7fe3rH/G9uO2H7M9OV8b6GECAAAAUEkvJ32wPSbpVkmXSZqWtMP2REQ82RJ2raQXIuJ825sl3Szpw7YvkLRZ0jskvVnS12z/u4ioF897b0Tsz7SDHiYAAAAA1URUX+a3QdJURDwdEYcl3S1p06yYTZK2FrfvlXSpmxc43CTp7og4FBH/ImmqKK80EiYAAAAAFSwgWWomTKtsT7YsW2ZVcLakZ1vuTxfr2sZExIyklySdPs9zQ9JXbT/Sps43YEgeMOyixNmWXkK/kTSSr7uWfM129bagtOjFEI7sPjECMzpVVuZ4AgALtz8i1s/xeLt/zrMP4p1i5nruJRGx2/YZkh6w/b2IeLhTI5bQtycAAAAAiybU6yF505LOabm/RtLuTjG2xyWdIunAXM+NiKN/90q6T/MM1SNhAgAAAFBNYwHL/HZIWmf7PNvL1ZzEYWJWzISka4rbV0p6MJrDGCYkbS5m0TtP0jpJ37Z9gu2TJMn2CZIul/TEXI1gSB4AAACASno5S15EzNi+TtL9ksYk3RERO23fJGkyIiYk3S7pTttTavYsbS6eu9P2FyQ9KWlG0m9HRN32mZLua84LoXFJd0XEV+ZqBwkTAAAAgGp6fF5oRGyTtG3Wuhtbbh+UdFWH535K0qdmrXta0jvLtIEheQAAAADQAT1MAAAAAMoLSY2lP/MoCRMAAACACtKz3Q01EiYAAAAA1ZAwAQAAAEAHJEwA+iZ7AHK7C1kvsOrkeGTXchdRUCM/v0x0eSoaN5JtrCUrzpZXpsxuq3V/n+jbGPUy73f2M5MsM7r9JaDMexjJNvbz3IER+JIEABIJEwAAAIAqmPQBAAAAADqJdI/4MCNhAgAAAFDNCAzP5cK1AAAAANABPUwAAAAAyuMcJgAAAACYwwgMySNhAgAAAFANCRMAAAAAtBMjkTAx6QMAAAAAdEAPEwAAAIDyQlKD6zABGHSlusK7e1CLRraTup4vNBnqmlNx6XfHydeSrFeS7C63sZZsY7LeUro95CL5DzbK1JudqanLF1mMfs4Q1e0LRo7A0BoAXTYCxw0SJgAAAADVkDABAAAAQDsxEtdhYtIHAAAAAOiAHiYAAAAA5YUU3T6XcgCRMAEAAACoZgSG5JEwAQAAAKhmBCZ94BwmAAAAAOiAHiYAAAAA5UVw4VoAAAAA6GgEhuSRMAGjpNsHtah3t7wyVffrBy07Hdr1fyEeglHU/ZwtaQT+aQPAoIkR6GFa0H9f2xttP2V7yvb13WoUAAAAgEEXzR+rqi5DonLCZHtM0q2SPiDpAklX276gWw0DAAAAgH5byJC8DZKmIuJpSbJ9t6RNkp7sRsMAAAAADLAQ12Gax9mSnm25Py3pZxfWHAAAAABDo5/nri6ShSRM7c58fkOKaXuLpC2SdO655y6gOgAAAACDIiTFCPQwLWTSh2lJ57TcXyNp9+ygiLgtItZHxPrVq1cvoDoAAAAAAyOi2cNUdRkSC0mYdkhaZ/s828slbZY00Z1mAQAAAED/VR6SFxEztq+TdL+kMUl3RMTOrrUMAAAAwEAbhSF5C7pwbURsk7StS20BAAAAMEyGaGhdVY5FvGiU7X2S/nXRKkQ7qyTt73cj8AZsl8HEdhlMbJfBxHYZTGyXwdRuu7wlIobqhH/bX1HztVS1PyI2dqs9vbKoCRP6z/ZkRKzvdztwLLbLYGK7DCa2y2BiuwwmtstgYrsMl4VM+gAAAAAASxoJEwAAAAB0QMI0em7rdwPQFttlMLFdBhPbZTCxXQYT22UwsV2GCOcwAQAAAEAH9DABAAAAQAckTCPC9lW2d9pu2F4/67FP2J6y/ZTt9/erjaPK9sbivZ+yfX2/2zOqbN9he6/tJ1rWvcn2A7a/X/w9rZ9tHDW2z7H9kO1dxfHrY8V6tksf2V5p+9u2v1Nslz8u1p9ne3uxXe6xvbzfbR1FtsdsP2r7H4r7bJc+s/2M7cdtP2Z7sljHcWyIkDCNjick/SdJD7eutH2BpM2S3iFpo6T/aXts8Zs3mor3+lZJH5B0gaSri22Cxfc5NT8Dra6X9PWIWCfp68V9LJ4ZSb8bEW+XdLGk3y4+H2yX/jok6X0R8U5JF0naaPtiSTdLuqXYLi9IuraPbRxlH5O0q+U+22UwvDciLmqZSpzj2BAhYRoREbErIp5q89AmSXdHxKGI+BdJU5I2LG7rRtoGSVMR8XREHJZ0t5rbBIssIh6WdGDW6k2Stha3t0q6YlEbNeIiYk9E/HNx+xU1vwSeLbZLX0XTq8XdZcUSkt4n6d5iPdulD2yvkfQfJf1Vcd9iuwwqjmNDhIQJZ0t6tuX+dLEOi4P3f7CdGRF7pOaXd0ln9Lk9I8v2WknvkrRdbJe+K4Z9PSZpr6QHJP1A0osRMVOEcCzrjz+X9N8lNYr7p4vtMghC0ldtP2J7S7GO49gQGe93A9A9tr8m6afaPPQHEfGlTk9rs46pExcP7z8wD9snSvqipI9HxMvNH83RTxFRl3SR7VMl3Sfp7e3CFrdVo832ByXtjYhHbL/n6Oo2oWyXxXdJROy2fYakB2x/r98NQjkkTEtIRPxChadNSzqn5f4aSbu70yIk8P4PtudtnxURe2yfpeav6VhEtpepmSz9TUT8XbGa7TIgIuJF299Q8xyzU22PF70ZHMsW3yWSftn2L0paKelkNXuc2C59FhG7i797bd+n5nB8jmNDhCF5mJC02fYK2+dJWifp231u0yjZIWldMYvRcjUn4Jjoc5vwExOSriluXyOpU08teqA4/+J2Sbsi4tMtD7Fd+sj26qJnSbaPk/QLap5f9pCkK4swtssii4hPRMSaiFir5v+SByPiI2K79JXtE2yfdPS2pMvVnIiL49gQ4cK1I8L2hyT9D0mrJb0o6bGIeH/x2B9I+i9qzkj18Yj4ct8aOoKKXwP/XNKYpDsi4lN9btJIsv23kt4jaZWk5yX9oaT/LekLks6V9G+SroqI2RNDoEds/5yk/yPpcf3knIzfV/M8JrZLn9j+92qepD6m5g+vX4iIm2y/Vc2Ja94k6VFJ/zkiDvWvpaOrGJL3exHxQbZLfxXv/33F3XFJd0XEp2yfLo5jQ4OECQAAAAA6YEgeAAAAAHRAwgQAAAAAHZAwAQAAAEAHJEwAAAAA0AEJEwAAAAB0QMIEAAAAAB2QMAEAAABAByRMAAAAANDB/wfVMdcQICtUowAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 1152x432 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"init()\n",
"aa_time_loop(time_steps)\n",
"vel_version1 = dh.gather_array(u.name, ghost_layers=False).copy()\n",
"plt.vector_field_magnitude(vel_version1[:, :, domain_size[2]//2, :])\n",
"plt.colorbar();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Version 2: Normal boundary handling"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0wAAAFlCAYAAADYoWhgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3df6xc5X3n8c9n7vUPfpPYhro21LR4uyF0QxTLRaKVklCI06Y12cLGNNsgLZL7CynZbXdFWkFblKyK1DZVGzYSKTQOagIRLRu3dUIIELGpIsfXhQSMg3LjknKxBXbMb+Ifd+a7f8xxOr7M3Ps9587cmbnzfklHnjnzned5Zs6Zc/09z3Oe44gQAAAAAOCNav1uAAAAAAAMKhImAAAAAOiAhAkAAAAAOiBhAgAAAIAOSJgAAAAAoAMSJgAAAADoYHwhK1u5cmWsW7duIasEBs4T+59LxS0543gq7tylL6frrifPkYypkYpb5ulkefnbF2Qjx+xUXE25uF7oX805LtHCSG6ZMmV2s95eyNbcSH833Xc02ch65H77S5z77UvSgeNnpeJee3l5Ku7i1eem6wYWo927dx+KiFX9bkcZ73nXafGDw/XK79/97aP3R8SmLjapJxY0YVq3bp0mJiYWskpg4PzHmz6Rijv3nc+m4j687sF03a81lqXiTqsdTcWtX3IwFXdGLX8wrSf/A3hGLfffz1O9JF13RjZRk6RxjXW17m4bc36QQT1y/5EuU2Y365XyiUs2iZ5Wbr99pXEsFbekxHczlmzjvtw5Cx2un5qK+/HxV3IFSrrl2V9Kxe3+ykWpuImb/3u6bmAxsv39frehrEOH69p5/9rK71+y+nsru9icnmFIHgAAAAB0sKA9TAAAAAAWiyg1EmBYkTABAAAAKC2UHw49zEiYAAAAAFTSSE4UNcy4hgkAAAAAOqCHCQAAAEBpoVA9GJIHAAAAAG1xDROArqufkjuw7P/mj6fiPvm//0u67sZY7t4uR1bm7l30wzfnRvVOn5YKkyQ1kkelxtJsXO77ztabvP9nUy1Xd7rMPt4JN7pct7N/X8v8HU7GZu/NWpvOfWhnbzNWYpj/+A9zdS9J3jZpyWu5L+eUEjeg3H9Z7j5jcebiv74BGFUhqU7CBAAAAADtjUIP05znNW0vt/1N29+yvcf2HxfrL7C90/Z3bd9jO3m+FwAAAACGQ2YgyFFJ746It0m6RNIm25dKulXSJyJivaQXJF3fu2YCAAAAGCQhqR5ReRkWcyZM0fRq8XRJsYSkd0u6t1i/TdJVPWkhAAAAgIHUmMcyLFLXMNkek7Rb0oWSbpP0PUkvRsR0ETIlaU1PWggAAABg4ISCSR9OiIi6pEtsny3pPklvaRfW7r22t0raKknnn39+xWYCAAAAGCgh1Rd/vpS6hulHIuJFSV+TdKmks22fSLjWStrf4T23R8SGiNiwatWq+bQVAAAAABZUZpa8VUXPkmyfIukXJO2V9LCkq4uw6yR9sVeNBAAAADBYQlzDdMJqSduK65hqkr4QEf9o+0lJd9v+mKRHJd3Rw3YCAAAAGChWvZ93VV8gcyZMEfFtSW9vs36fpI29aBSwmP3UtudScXv/54pU3L/+cv4WaOs//1oqbvmzL6fizp6u5yp++dW5Ywo+9ZRUXIyP5QpMTlsay5fl4pYk65Xyg56d+2MTtT7+UcrW3ch9306OeY8SH9n13PnK2g+P5+qu5Tagk/uYf/BiKk6SNJbbz+LM03LlJT/Ld37r7Fx5kk5ZnTtOrPmz1OXS0kfSVQMYEKH0YX+olbqGCQAAAABGSfK0DwAAAACcjCF5AAAAANBGiIQJAAAAADpqlLnYdEiRMAEAAAAobVR6mJj0AQAAAAA6oIcJAAAAQGkhqz4C/S8kTAAAAAAq4RomAAAAAGhjVK5hImECAAAAUIFVD4bkAeg2587ErNiZ+3kePz1/ZucHF5+eijv9wPJU3NLDR1JxPvOUVJwk1V4/livz1R/mCjx+PBf34su5eo8ly5MU9XoyMHJ1NxrJ4nLlqZGMK6OW2x+d/B24lv9D7PHkn7QlubjakiW58k7J/V5ixdm58iQ1lubaeOTHTk3FvfSTuc+y9HB+n/DUWam48RcPpcsEgEFEwgQAAACgtJDUGIFJHxb/JwQAAADQE3W58pJhe5Ptp2xP2r6xzevLbN9TvL7T9rpi/UbbjxXLt2y/v+U9T9t+vHhtYq420MMEAAAAoLSI3l7DZHtM0m2SrpA0JWmX7e0R8WRL2PWSXoiIC21vkXSrpA9IekLShoiYtr1a0rds/0NETBfve1dEpMYM08MEAAAAYBBtlDQZEfsi4pikuyVtnhGzWdK24vG9ki637Yh4vSU5Wq7mCMJKSJgAAAAAVNKQKy8JayQ90/J8qljXNqZIkF6StEKSbP+s7T2SHpf0my0JVEj6iu3dtrfO1QiG5AEAAAAorXkfpnn1v6yccQ3R7RFxe8vzdlnVzJ6ijjERsVPSW22/RdI221+KiCOSLouI/bbPkfSA7e9ExCOdGknCBAAAAKCCeV/DdCgiNszy+pSk81qer5W0v0PMlO1xSWdJOtwaEBF7bb8m6WJJExGxv1j/vO371Bz61zFhYkgeAAAAgNJOTCtedUnYJWm97QtsL5W0RdL2GTHbJV1XPL5a0kMREcV7xiXJ9k9I+mlJT9s+zfYZxfrTJF2p5gQRHdHDBAAAAGDgFDPc3SDpfkljku6MiD22b1Gzp2i7pDsk3WV7Us2epS3F239O0o22j0tqSPrtiDhk+ycl3VfcQH1c0uci4suztYOECVhgjTOWp+LO/P6xVFyM5zuKXzs395M/dkauzPFXx1Jxrufb2Fi+JBVXW740F/f60XTdKY1GOrQ2Xc8FRuWJexamPEly7n4Z6bqT5cWSEn+msm1Mapy6LBVXPzW3z5aqe1nut/XDVbnvp5H8Gt+0N7/vjB/N/RYap+W+RwDDqR7dPfbOFBE7JO2Yse7mlsdHJF3T5n13Sbqrzfp9kt5Wpg0kTAAAAABKC3m+kz4MBRImAAAAAJU0enjj2kFBwgQAAACgtC5MKz4UFv8nBAAAAICK6GECAAAAUFrIPZ/0YRCQMAEAAACoJHk/paFGwgQAAACgtAipPgKTPiz+TwgAAAAAFdHDBAAAAKACqyGuYQLQZWOHXk7F+ezlqbgoUfcZU8dyZSaPffXlY6m48dencwVKqh05no7NiOVLUnE+Xk8WmP/GY0nyEJsts0TdfePkzpONK2M8tz/GkmTcWK6NteONVFz91Pyf3OOn52KXvpKre9nLuX0nSmyXsaO5umsvvJouE8BwCY3GkDwSJgAAAACVjMJ9mEiYAAAAAJQWshojMK344k8JAQAAAKAiepgAAAAAVMKQPAAAAABoIyQ1mPQBAAAAANqx6kwrDgAAAABvNCo9TIv/EwIAAABARfQwAQAAAKiEIXkAuq5x6HAqbsk5Z6fijp21NF332OvTqThHpOIa47lOak83UnGS5OP1dGxKtup6MjD53ZSSLNONXFzUuv/HK19316vOy27D7H5bT37mRq5eT+f3nbEf5n4HS0qUmTF92lg6dulLx3OBz/+gYmsADLoIj8SQPBImAAAAAJXURyBhmvMT2j7P9sO299reY/vDxfo/sv2s7ceK5Rd731wAAAAAWDiZHqZpSb8bEf9i+wxJu20/ULz2iYj40941DwAAAMAgCkkNrmGSIuKApAPF41ds75W0ptcNAwAAADDIzJC8mWyvk/R2STuLVTfY/rbtO22/qcttAwAAADCgmvdhcuVlWKQTJtunS/o7SR+JiJclfUrST0m6RM0eqD/r8L6ttidsTxw8eLALTQYAAAAwCOqqVV6GRaqltpeomSz9bUT8vSRFxHMRUY+IhqRPS9rY7r0RcXtEbIiIDatWrepWuwEAAACg5+a8hsm2Jd0haW9E/HnL+tXF9U2S9H5JT/SmiQAAAAAGTWi4htZVlZkl7zJJvy7pcduPFet+X9K1ti9Rc/ji05J+oyctBAAAADCQGkM0tK6qzCx5X5fazhe4o/vNAQAAADAMIqQ6PUwAus3Ll+XiXnwtFbfs6PF03Y1Tl+bqnm6k4sZfPZaruMzJJ+cOvD42nStvul6i8oSIdKgbydgSZabqrXe3vGahye1Sz+072fJirMTOk/wefST5m1mS+xMZy3NxtaPJfVbS0tdzbWwk646x3Pe9/Pn88aT26pFcYPJ7BDCcRmFI3uLvQwMAAACAijjtAwAAAKC05qQPi7//ZfF/QgAAAAA9UZcrLxm2N9l+yvak7RvbvL7M9j3F6zttryvWb7T9WLF8y/b7s2XORA8TAAAAgNJCvb2GyfaYpNskXSFpStIu29sj4smWsOslvRARF9reIulWSR9Q85ZHGyJi2vZqSd+y/Q9Fs+cq8yT0MAEAAAAYRBslTUbEvog4JuluSZtnxGyWtK14fK+ky207Il6PiBOz7SxXM1HKlnkSEiYAAAAAFTSvYaq6SFppe6Jl2TqjgjWSnml5PlWsaxtTJEgvSVohSbZ/1vYeSY9L+s3i9UyZJ2FIHgAAAIBKGslrkTo4FBEbZnm9XeEz7yHRMSYidkp6q+23SNpm+0vJMk9CwgQAAACgtAW4ce2UpPNanq+VtL9DzJTtcUlnSTp8cjtjr+3XJF2cLPMkDMkDAAAAUMk8h+TNZZek9bYvsL1U0hZJ22fEbJd0XfH4akkPRUQU7xmXJNs/IemnJT2dLPMk9DABC60xa6/vj/jY8Vx5kStP6v4ZEtfrucBkWCnTuUJdb3S33hLfd6nYbupFve7yGcRkG13is8RYcg/PfpbkPlY7kiuujEi2sdblH3XtlfyHyR6jInnMA4CZihnubpB0v6QxSXdGxB7bt0iaiIjtku6QdJftSTV7lrYUb/85STfaPi6pIem3I+KQJLUrc7Z2kDABAAAAKK1549qeDslTROyQtGPGuptbHh+RdE2b990l6a5smbMhYQIAAABQyTwnfRgKJEwAAAAASuv1jWsHBZM+AAAAAEAH9DABAAAAqCQ5291QI2ECAAAAUF70ftKHQUDCBAAAAKC0EJM+AAAAAEBHo9DDtPgHHQIAAABARfQwAQus8eprqbixU0/JFbgk/zOuvXokFRe15NmisbFcXL2ei5PkeiMZ2OUzWtk2NiJfZvZ7zMrW3e16JSm7Xbpdd+S/b2dja7lzhZEtL7vvJOuVJI/nflt+dTpXYKMHv6vkPtF4/fV8mQCGyqhMK07CBAAAAKASEiYAAAAAaCPELHkAAAAA0NEozJLHpA8AAAAA0AE9TAAAAADKC65hAgAAAIC2mCUPAAAAAGYxCgkT1zABAAAAQAf0MAEAAAAojWnFAfRVHD2ainOtREdxLXdQcz1Z3nQ2sISI7pbX6HJ5ZdpX73Ldg15vL+p2mT/EybqT29CNZN3ZNpbYF6PRyFVdz8V1/XclKY7kjlFd/w0CGChBwgQAAAAA7Y3CfZhImAAAAACUFiMyrTiTPgAAAABAB/QwAQAAAKiEa5gAAAAAoC1myQMAAACAjuhhAgAAAIA2Qkz6AAAAAAAjjR4mAAAAAOVFT+6LPXBImICFVst1XcfRY12v2suXd7fA7FHSJbrrs7HZurvdxjJ/GRqNXFxtBDv7e/HdZLdhvcv7RPI3rUZ+33F2PxtLfj/JuuPI0Vx5knQ0GZv9fgAMJW5cCwAAAABthEZj0oc5T03ZPs/2w7b32t5j+8PF+jfbfsD2d4t/39T75gIAAADAwsn05U9L+t2IeIukSyX9ju2LJN0o6cGIWC/pweI5AAAAgJHQvA9T1WVYzJkwRcSBiPiX4vErkvZKWiNps6RtRdg2SVf1qpEAAAAABk9E9SXD9ibbT9metP2GDhrby2zfU7y+0/a6Yv0Vtnfbfrz4990t7/laUeZjxXLObG0odQ1T0YC3S9op6dyIOCA1k6q5KgIAAACwuPTyGibbY5Juk3SFpClJu2xvj4gnW8Kul/RCRFxoe4ukWyV9QNIhSb8cEfttXyzpfjU7fU74YERMZNqRnn7I9umS/k7SRyLi5RLv22p7wvbEwYMHs28DAAAAMMCaPUWuvCRslDQZEfsi4piku9Uc5daqddTbvZIut+2IeDQi9hfr90habntZlc+ZSphsL1EzWfrbiPj7YvVztlcXr6+W9Hy790bE7RGxISI2rFq1qkobAQAAAIyeNZKeaXk+pZN7iU6KiYhpSS9JWjEj5lclPRoRrfdD+JtiON5N9uz3kcjMkmdJd0jaGxF/3vLSdknXFY+vk/TFucoCAAAAsHjMc9KHlSdGohXL1hnFt0tkZl79NGuM7beqOUzvN1pe/2BE/Iykny+WX5/tM2auYbqsKORx248V635f0p9I+oLt6yX9m6RrEmUBAAAAWCTK3M+9jUMRsWGW16ckndfyfK2k/R1ipmyPSzpL0mFJsr1W0n2SPhQR3/v3Nsezxb+v2P6cmkP/PtupEXMmTBHxdbXP3CTp8rneDwAAAGBx6vGNa3dJWm/7AknPStoi6ddmxJwY9fYNSVdLeigiwvbZkv5J0kcj4p9PBBdJ1dkRcai47Oh9kr46WyNKzZIHoAvq9Vxco5EsLxknKY4fT8XNMZS3NTBdd9fVknPWZL/HrHmeSmsru0/08/vO6vb30+3tJ+X3nexn6UET047lftP5z5LcFyVFtszs/g1g6ITSkzdUKz9i2vYNas5wNybpzojYY/sWSRMRsV3NS4fusj2pZs/SluLtN0i6UNJNtm8q1l0p6TVJ9xfJ0piaydKnZ2sHCRMAAACAgRQROyTtmLHu5pbHR9Tm0qCI+Jikj3Uo9h1l2kDCBAAAAKCSHoy7GDgkTAAAAADKi55fwzQQSJgAAAAAVDMCXUzJK18BAAAAYPTQwwQAAACgEobkAQAAAEAHvbjbxqAhYQIAAABQWogeJgAAAABoLySRMAHotkje9d696OOenk6FhXMHP493/xASyc/d7e8nW29fZb+b5PYrV3Wfvp8S9aY/d6NRsTEd9OK7SX6WmM4dTxTJz9wo8VmSx7LsMQ8ABhUJEwAAAIBKhuF843yRMAEAAACohoQJAAAAANoxkz4AAAAAQEcj0MNU63cDAAAAAGBQ0cMEAAAAoLzgPkwAAAAA0NkIDMkjYQIAAABQ0eLvYeIaJgAAAADogB4mYKE5d54i6vVcedPT+arHkz/55KmUyNad/MzNunNnqqJfd8prDP7YgxiF8RFtRHY36/K+Y/dgn603cnGRjUvWnT3uqMQxqszvH8DwGYE/OSRMAAAAAKohYQIAAACANkISs+QBAAAAQHv9GiG/kBhYDAAAAAAd0MMEAAAAoJoR6GEiYQIAAABQDdcwAQAAAEB7pocJAAAAANoIjcSQPCZ9AAAAAIAO6GEChl2jkY+NbOxYsrzsaaUSbWz06TxO+rspU2by+/HiH/+9IPq070S22kaJ07DZ/TG7jyXrjno9Vx4ASJLMNUwAAAAA0NEIDMkjYQIAAABQzQgkTFzDBAAAAAAd0MMEAAAAoBp6mAAAAACgjVBz0oeqS4LtTbafsj1p+8Y2ry+zfU/x+k7b64r1V9jebfvx4t93t7znHcX6Sdt/ac8++xIJEwAAAIBKHNWXOcu2xyTdJum9ki6SdK3ti2aEXS/phYi4UNInJN1arD8k6Zcj4mckXSfprpb3fErSVknri2XTbO0gYQIAAABQTcxjmdtGSZMRsS8ijkm6W9LmGTGbJW0rHt8r6XLbjohHI2J/sX6PpOVFb9RqSWdGxDciIiR9VtJVszWChAkAAABAP6y0PdGybJ3x+hpJz7Q8nyrWtY2JiGlJL0laMSPmVyU9GhFHi/ipOco8CZM+AAAAAOiHQxGxYZbX211bNLNvatYY229Vc5jelSXKPAk9TAAAAAAq6eU1TGr2/pzX8nytpP2dYmyPSzpL0uHi+VpJ90n6UER8ryV+7RxlnoQeJmBAzTFhSzWN5NyfyaOYask2Rpk5RxslYruoVBsXUd2LSdRzcd3+bSWrLSW7T2R/0/XuNzJ7jGLvBha55Gx3Fe2StN72BZKelbRF0q/NiNmu5qQO35B0taSHIiJsny3pnyR9NCL++UfNjThg+xXbl0raKelDkv5qtkbM2cNk+07bz9t+omXdH9l+1vZjxfKLmU8MAAAAYJGYz4QPibMpxTVJN0i6X9JeSV+IiD22b7H9K0XYHZJW2J6U9D8knZh6/AZJF0q6qSVnOad47bck/bWkSUnfk/Sl2dqR6WH6jKRPqjmDRKtPRMSfJt4PAAAAAKVFxA5JO2asu7nl8RFJ17R538ckfaxDmROSLs62Yc4epoh4RMU4QAAAAAD4kd5OKz4Q5jPpww22v10M2XtT11oEAAAAYCj0eNKHgVA1YfqUpJ+SdImkA5L+rFOg7a0n5lY/ePBgxeoAAAAADBx6mNqLiOcioh4RDUmfVvMuvJ1ib4+IDRGxYdWqVVXbCQAAAAALrlLCZHt1y9P3S3qiUywAAACARWoEepjmnCXP9uclvVPSSttTkv5Q0jttX6LmR31a0m/0sI0AAAAABsywXYtU1ZwJU0Rc22b1HT1oCwAAAIBh0tsb1w6EzH2YAPRBRO6UTU8OU9HIxTXmM9FmBz0oMqUxAqfIFrta8teQ/G31VXZ/zP5WAaBXhuCQOl/9+q8JAAAAAAw8epgAAAAAVMI1TAAAAADQCQkTAAAAALQxIrPkcQ0TAAAAAHRADxMAAACAakagh4mECQAAAEA1JEwAAAAA0B7XMAEAAADACKOHCVhgrrm7BUb+1E4kY93IlljPhZX5zI0BP48T6S8H3eAS+0NjwE9z9mLfSX7m7G+/J8eTbh/zAGCBkTABAAAAqGbAz1V1AwkTAAAAgPJG5D5MJEwAAAAAqiFhAgAAAIAORiBhGvCrqwEAAACgf+hhAgAAAFCaxTVMAAAAANAZCRMAAAAAtDEis+RxDRMAAAAAdEAPEzDsGo18rJ0Ki1p3z6W4RBNVKxO8SDQW0em5Wm4fS4tFtD/0YDtHJMvMHifKHE8AQGJIHgAAAAB0RMIEAAAAAO2NwjVMJEwAAAAAqhmBhIlJHwAAAACgAxImAAAAAOXFPJcE25tsP2V70vaNbV5fZvue4vWdttcV61fYftj2q7Y/OeM9XyvKfKxYzpmtDQzJAwAAAFBJL69hsj0m6TZJV0iakrTL9vaIeLIl7HpJL0TEhba3SLpV0gckHZF0k6SLi2WmD0bERKYd9DABAAAAqKa3PUwbJU1GxL6IOCbpbkmbZ8RslrSteHyvpMttOyJei4ivq5k4zQsJEwAAAIBKHNWXhDWSnml5PlWsaxsTEdOSXpK0IlH23xTD8W6yZ79RJQkTAAAAgH5YaXuiZdk64/V2iczMVCsTM9MHI+JnJP18sfz6bMFcwwQssGjkTql4rBeVJwcaNxrdrbaWPzfj7lYt1WY9adRbyW29qHT7M/dz+2UlP3Nkf3+l6k7+YHpRd1L2mAdgSM3vJ34oIjbM8vqUpPNanq+VtL9DzJTtcUlnSTo8W6UR8Wzx7yu2P6fm0L/PdoqnhwkAAABAeb2fJW+XpPW2L7C9VNIWSdtnxGyXdF3x+GpJD8UsZ6hsj9teWTxeIul9kp6YrRH0MAEAAAAozWo/Hq5bImLa9g2S7pc0JunOiNhj+xZJExGxXdIdku6yPalmz9KWH7XPflrSmZKW2r5K0pWSvi/p/iJZGpP0VUmfnq0dJEwAAAAABlJE7JC0Y8a6m1seH5F0TYf3rutQ7DvKtIGECQAAAEA1I3CZIgkTAAAAgEp6eePaQUHCBAAAAKAaEiYAAAAA6GAEEiamFQcAAACADuhhAgAAAFBecA0TAAAAAHRGwgRg0M1yM+s3cKORC6x1ebRutl5JkazbTt4qr7F4juRltnW/pLdLVpntV+vPPtGT7ZL9zWTrTpY3DPsYgMEyCj1Mc/7PxPadtp+3/UTLujfbfsD2d4t/39TbZgIAAAAYODGPZUhkTuV+RtKmGetulPRgRKyX9GDxHAAAAAAWlTkTpoh4RNLhGas3S9pWPN4m6aoutwsAAADAgHNUX4ZF1WuYzo2IA5IUEQdsn9PFNgEAAAAYdEM2tK6qnt+HyfZW2xO2Jw4ePNjr6gAAAAAsFK5h6ug526slqfj3+U6BEXF7RGyIiA2rVq2qWB0AAAAALLyqCdN2SdcVj6+T9MXuNAcAAADAMLBG4xqmzLTin5f0DUk/bXvK9vWS/kTSFba/K+mK4jkAAACAUTICQ/LmnPQhIq7t8NLlXW4LAAAAgCHiEbjhddVZ8gD0WqP7B6BIDsJ1o5ELrCULLHMwTdYdybpt5+vuk1hEf2yyn6Un26XLv5mub5fs76pZeVfLTH+WHhx3ACxiQ9ZTVFXPZ8kDAAAAgGFFDxMAAACASoZp8oaqSJgAAAAAVEPCBAAAAADt0cMEAAAAAJ2MQMLEpA8AAAAA0AE9TAAAAADKC4bkAQAAAEBnJEwAAAAA8EYWPUwAeiEayTjn4lziUsRG7qgWySLdSH6WWv8ul4zIfWY7+X0Pg+x2KaOP27DbsvtEWvb7LlNvssyufxYAwBuQMAEAAACoZgRO3JAwAQAAAKiEIXkAAAAA0E6ISR8AAAAAoBP34LLZQbN4ruIFAAAAgC4jYQIAAABQTcxjSbC9yfZTtidt39jm9WW27yle32l7XbF+he2Hbb9q+5Mz3vMO248X7/lLzzFVLgkTAAAAgEoc1Zc5y7bHJN0m6b2SLpJ0re2LZoRdL+mFiLhQ0ick3VqsPyLpJkm/16boT0naKml9sWyarR0kTAAAAADKCzWnFa+6zG2jpMmI2BcRxyTdLWnzjJjNkrYVj++VdLltR8RrEfF1NROnH7G9WtKZEfGNaN7M7rOSrpqtESRMAAAAACqZZw/TStsTLcvWGcWvkfRMy/OpYl3bmIiYlvSSpBWzNHlNUc5sZZ6EWfKAYRclpqdxd8+RRPJmdW6UaGMt2cZsmcny0p9l9mHOWEyy+1j2po0lfgfZ/TGtkSyvzPEkqxdlAlgsDkXEhlleb/dHd+YBLRMzn3h6mAAAAABU1NtJH6YkndfyfK2k/Z1ibI9LOnw1Rl0AAAxhSURBVEvS4TnKXDtHmSchYQIAAABQmtXbSR8k7ZK03vYFtpdK2iJp+4yY7ZKuKx5fLemhmKWbPiIOSHrF9qXF7HgfkvTF2RrBkDwAAAAA5eUnb6hYfEzbvkHS/ZLGJN0ZEXts3yJpIiK2S7pD0l22J9XsWdpy4v22n5Z0pqSltq+SdGVEPCnptyR9RtIpkr5ULB2RMAEAAAAYSBGxQ9KOGetubnl8RNI1Hd67rsP6CUkXZ9tAwgQAAACgkuTQuqFGwgQAAACgGhImAAAAAGiPHiYAAAAAaCeUv8/bEGNacQAAAADogB4mYEBF8oyNa+1uWN2x0FxcI3kupUzd3ZadxrSR/My13Gee5dYOKIHvsUuyZ3azv/0SsscoAIvcCBwKSJgAAAAAVMI1TAAAAADQyQiMGCBhAgAAAFDJKPQwMekDAAAAAHRADxMAAACA8kJM+gAAAAAA7ViSuYYJAAAAADro/l0LBg7XMAEAAABAB/QwAQAAAKiEIXkAAAAA0A6TPgAYBtHIH6lcc3crT9YdJQb/upEcDF1jRHFHfDezy+5j2bOmyfKizFnYEr/rbipzPAEAKbhx7VxsPy3pFUl1SdMRsaEbjQIAAAAw+EbhxrXd6GF6V0Qc6kI5AAAAADBQGJIHAAAAoJoRGJI334HuIekrtnfb3tqNBgEAAAAYAiG5UX0ZFvPtYbosIvbbPkfSA7a/ExGPtAYUidRWSTr//PPnWR0AAACAgUEP0+wiYn/x7/OS7pO0sU3M7RGxISI2rFq1aj7VAQAAAMCCqpww2T7N9hknHku6UtIT3WoYAAAAgAEX81iGxHyG5J0r6T7bJ8r5XER8uSutAgAAADDwPAJD8ionTBGxT9LbutgWAAAAAMOEhAlA16UPLMnpY5wfWRuNXN2udb/uvsl+343kZ64NwWceVdltuJhE7jNnf/u9qBvAIhZK/3dlmPGXHwAAAAA6oIcJAAAAQGlWcA0TAAAAAHREwgQAAAAAHZAwAQAAAEAbTPoAAAAAAP1je5Ptp2xP2r6xzevLbN9TvL7T9rqW1z5arH/K9nta1j9t+3Hbj9memKsN9DABAAAAqKSXkz7YHpN0m6QrJE1J2mV7e0Q82RJ2vaQXIuJC21sk3SrpA7YvkrRF0lsl/bikr9r+DxFRL973rog4lGkHPUwAAAAAqomovsxto6TJiNgXEcck3S1p84yYzZK2FY/vlXS5bRfr746IoxHxr5Imi/JKI2ECAAAAUME8kqVmwrTS9kTLsnVGBWskPdPyfKpY1zYmIqYlvSRpxRzvDUlfsb27TZ1vwJA8YNhFiast3adzJI18d30km+hG8nPXuvyZs/X2om50R3b4SJltvViUOZ6ky1z8M2gBqOxQRGyY5XW3WTfzoNIpZrb3XhYR+22fI+kB29+JiEc6NYK/5gAAAADKC/V6SN6UpPNanq+VtL9TjO1xSWdJOjzbeyPixL/PS7pPcwzVI2ECAAAAUE1jHsvcdklab/sC20vVnMRh+4yY7ZKuKx5fLemhiIhi/ZZiFr0LJK2X9E3bp9k+Q5JsnybpSklPzNYIhuQBAAAAqKSXs+RFxLTtGyTdL2lM0p0Rscf2LZImImK7pDsk3WV7Us2epS3Fe/fY/oKkJyVNS/qdiKjbPlfSfc15ITQu6XMR8eXZ2kHCBAAAAKCaHl+nGBE7JO2Yse7mlsdHJF3T4b0fl/TxGev2SXpbmTYwJA8AAAAAOqCHCQAAAEB5oVIz4Q4rEiYAAAAAFaRnuxtqJEwAAAAAqiFhAgAAAIAOSJgA9E32AOR2N7KeZ9XJ8ciu5W6iIHd/fplIfj9uJNtYS7axzPfd7bqT3It9ott/ELPfTRnZNibr7v5nLlFeJNvYz2sHRuA/SQAgkTABAAAAqIJJHwAAAACgk0j3iA8zEiYAAAAA1YzA8FxuXAsAAAAAHdDDBAAAAKA8rmECAAAAgFmMwJA8EiYAAAAA1ZAwAQAAAEA7MRIJE5M+AAAAAEAH9DABAAAAKC8kNbgPE4BBV6orvLsHtagnA12ijckyXXMqrusDBVyiYz7ZxnTV7tNnLiG6PTSjzOxLfbp5YvRihqh+3QhyBIbWAOiyEThukDABAAAAqIaECQAAAADaiZG4DxOTPgAAAABAB/QwAQAAACgvpOjXNZcLiIQJAAAAQDUjMCSPhAkAAABANSMw6QPXMAEAAABAB/QwAQAAACgvghvXAgAAAEBHIzAkj4QJGCX9OqhFvftFLqYTWnYqbPH/SWpjBP4QA8AwixHoYZrXNUy2N9l+yvak7Ru71SgAAAAAgy6aJ7aqLkOicsJke0zSbZLeK+kiSdfavqhbDQMAAACAfpvPkLyNkiYjYp8k2b5b0mZJT3ajYQAAAAAGWIj7MM1hjaRnWp5PSfrZ+TUHAAAAwNBYVBcVtzefhKndVcpvSDFtb5W0VZLOP//8eVQHAAAAYFCEpBiBHqb5TPowJem8ludrJe2fGRQRt0fEhojYsGrVqnlUBwAAAGBgRDR7mKouQ2I+CdMuSettX2B7qaQtkrZ3p1kAAAAA0H+Vh+RFxLTtGyTdL2lM0p0RsadrLQMAAAAw0EZhSN68blwbETsk7ehSWwAAAAAMkyEaWleVYwFvGmX7oKTvL1iFaGelpEP9bgTegO0ymNgug4ntMpjYLoOJ7TKY2m2Xn4iIobrg3/aX1fwsVR2KiE3dak+vLGjChP6zPRERG/rdDpyM7TKY2C6Die0ymNgug4ntMpjYLsNlPpM+AAAAAMCiRsIEAAAAAB2QMI2e2/vdALTFdhlMbJfBxHYZTGyXwcR2GUxslyHCNUwAAAAA0AE9TAAAAADQAQnTiLB9je09thu2N8x47aO2J20/Zfs9/WrjqLK9qfjuJ23f2O/2jCrbd9p+3vYTLevebPsB298t/n1TP9s4amyfZ/th23uL49eHi/Vslz6yvdz2N21/q9guf1ysv8D2zmK73GN7ab/bOopsj9l+1PY/Fs/ZLn1m+2nbj9t+zPZEsY7j2BAhYRodT0j6z5IeaV1p+yJJWyS9VdImSf/H9tjCN280Fd/1bZLeK+kiSdcW2wQL7zNq/gZa3SjpwYhYL+nB4jkWzrSk342It0i6VNLvFL8Ptkt/HZX07oh4m6RLJG2yfamkWyV9otguL0i6vo9tHGUflrS35TnbZTC8KyIuaZlKnOPYECFhGhERsTcinmrz0mZJd0fE0Yj4V0mTkjYubOtG2kZJkxGxLyKOSbpbzW2CBRYRj0g6PGP1ZknbisfbJF21oI0acRFxICL+pXj8ipr/CVwjtktfRdOrxdMlxRKS3i3p3mI926UPbK+V9EuS/rp4brFdBhXHsSFCwoQ1kp5peT5VrMPC4PsfbOdGxAGp+Z93Sef0uT0jy/Y6SW+XtFNsl74rhn09Jul5SQ9I+p6kFyNiugjhWNYffyHpf0lqFM9XiO0yCELSV2zvtr21WMdxbIiM97sB6B7bX5X0Y21e+oOI+GKnt7VZx9SJC4fvH5iD7dMl/Z2kj0TEy82T5uiniKhLusT22ZLuk/SWdmEL26rRZvt9kp6PiN2233lidZtQtsvCuywi9ts+R9IDtr/T7wahHBKmRSQifqHC26YkndfyfK2k/d1pERL4/gfbc7ZXR8QB26vVPJuOBWR7iZrJ0t9GxN8Xq9kuAyIiXrT9NTWvMTvb9njRm8GxbOFdJulXbP+ipOWSzlSzx4nt0mcRsb/493nb96k5HJ/j2BBhSB62S9pie5ntCyStl/TNPrdplOyStL6YxWipmhNwbO9zm/Dvtku6rnh8naROPbXogeL6izsk7Y2IP295ie3SR7ZXFT1Lsn2KpF9Q8/qyhyVdXYSxXRZYRHw0ItZGxDo1/5Y8FBEfFNulr2yfZvuME48lXanmRFwcx4YIN64dEbbfL+mvJK2S9KKkxyLiPcVrfyDpv6k5I9VHIuJLfWvoCCrOBv6FpDFJd0bEx/vcpJFk+/OS3ilppaTnJP2hpP8r6QuSzpf0b5KuiYiZE0OgR2z/nKT/J+lx/fs1Gb+v5nVMbJc+sf2f1LxIfUzNE69fiIhbbP+kmhPXvFnSo5L+a0Qc7V9LR1cxJO/3IuJ9bJf+Kr7/+4qn45I+FxEft71CHMeGBgkTAAAAAHTAkDwAAAAA6ICECQAAAAA6IGECAAAAgA5ImAAAAACgAxImAAAAAOiAhAkAAAAAOiBhAgAAAIAOSJgAAAAAoIP/D1M7BqidCM9xAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 1152x432 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"ldc = create_lid_driven_cavity(domain_size, relaxation_rate=relaxation_rate, lid_velocity=lid_velocity)\n",
"ldc.run(time_steps)\n",
"vel_version2 = ldc.velocity[:, :, :, :]\n",
"\n",
"plt.vector_field_magnitude(vel_version2[:, :, domain_size[2]//2, :])\n",
"plt.colorbar();"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.9"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
%% Cell type:code id: tags:
```
python
from
lbmpy.session
import
*
from
lbmpy.boundaries.boundaries_in_kernel
import
update_rule_with_push_boundaries
from
lbmpy.macroscopic_value_kernels
import
macroscopic_values_getter
,
macroscopic_values_setter
from
collections
import
OrderedDict
from
time
import
perf_counter
```
%% Cell type:markdown id: tags:
# Version 1: compile-in boundaries
%% Cell type:code id: tags:
```
python
domain_size
=
(
32
,
32
,
32
)
relaxation_rate
=
1.8
time_steps
=
100
lid_velocity
=
0.05
```
%% Cell type:code id: tags:
```
python
dh
=
create_data_handling
(
domain_size
,
default_target
=
'
cpu
'
)
pdfs
=
dh
.
add_array
(
'
pdfs
'
,
values_per_cell
=
19
)
u
=
dh
.
add_array
(
'
u
'
,
values_per_cell
=
len
(
domain_size
))
```
%% Cell type:code id: tags:
```
python
boundaries
=
OrderedDict
((
((
0
,
1
,
0
),
UBB
([
lid_velocity
,
0
,
0
])),
((
1
,
0
,
0
),
NoSlip
()),
((
-
1
,
0
,
0
),
NoSlip
()),
((
0
,
-
1
,
0
),
NoSlip
()),
((
0
,
0
,
1
),
NoSlip
()),
((
0
,
0
,
-
1
),
NoSlip
()),
))
opt
=
{
'
symbolic_field
'
:
pdfs
,
'
cse_global
'
:
False
,
'
cse_pdfs
'
:
True
}
cr_even
=
create_lb_collision_rule
(
stencil
=
"
D3Q19
"
,
relaxation_rate
=
relaxation_rate
,
compressible
=
False
,
optimization
=
opt
)
cr_odd
=
create_lb_collision_rule
(
stencil
=
"
D3Q19
"
,
relaxation_rate
=
relaxation_rate
,
compressible
=
False
,
optimization
=
opt
)
update_rule_aa_even
=
update_rule_with_push_boundaries
(
cr_even
,
pdfs
,
boundaries
,
AAEvenTimeStepAccessor
,
AAOddTimeStepAccessor
.
read
)
update_rule_aa_odd
=
update_rule_with_push_boundaries
(
cr_odd
,
pdfs
,
boundaries
,
AAOddTimeStepAccessor
,
AAEvenTimeStepAccessor
.
read
)
getter_assignments
=
macroscopic_values_getter
(
update_rule_aa_even
.
method
,
velocity
=
u
.
center_vector
,
pdfs
=
pdfs
.
center_vector
,
density
=
None
)
getter_kernel
=
ps
.
create_kernel
(
getter_assignments
,
target
=
dh
.
default_target
).
compile
()
even_kernel
=
ps
.
create_kernel
(
update_rule_aa_even
,
target
=
dh
.
default_target
,
ghost_layers
=
1
).
compile
()
odd_kernel
=
ps
.
create_kernel
(
update_rule_aa_odd
,
target
=
dh
.
default_target
,
ghost_layers
=
1
).
compile
()
```
%% Cell type:code id: tags:
```
python
def
init
():
dh
.
fill
(
pdfs
.
name
,
0
,
ghost_layers
=
True
)
def
aa_time_loop
(
steps
=
100
):
assert
steps
%
2
==
0
,
"
Works only for an even number of time steps
"
dh
.
all_to_gpu
()
for
i
in
range
(
steps
//
2
):
dh
.
run_kernel
(
odd_kernel
)
dh
.
run_kernel
(
even_kernel
)
dh
.
run_kernel
(
getter_kernel
)
dh
.
all_to_cpu
()
```
%% Cell type:code id: tags:
```
python
init
()
aa_time_loop
(
time_steps
)
vel_version1
=
dh
.
gather_array
(
u
.
name
,
ghost_layers
=
False
).
copy
()
plt
.
vector_field_magnitude
(
vel_version1
[:,
:,
domain_size
[
2
]
//
2
,
:])
plt
.
colorbar
();
```
%% Output
%% Cell type:markdown id: tags:
# Version 2: Normal boundary handling
%% Cell type:code id: tags:
```
python
ldc
=
create_lid_driven_cavity
(
domain_size
,
relaxation_rate
=
relaxation_rate
,
lid_velocity
=
lid_velocity
)
ldc
.
run
(
time_steps
)
vel_version2
=
ldc
.
velocity
[:,
:,
:,
:]
plt
.
vector_field_magnitude
(
vel_version2
[:,
:,
domain_size
[
2
]
//
2
,
:])
plt
.
colorbar
();
```
%% Output
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment