From b1250a857b8ed06121f1cfbb77d89f4d80f75d05 Mon Sep 17 00:00:00 2001
From: Markus Holzer <markus.holzer@fau.de>
Date: Fri, 9 Jul 2021 08:14:40 +0000
Subject: [PATCH] FreeSlip

---
 lbmpy/advanced_streaming/indexing.py   |  28 +++
 lbmpy/boundaries/__init__.py           |   4 +-
 lbmpy/boundaries/boundaryconditions.py | 141 ++++++++++-
 lbmpy_tests/test_boundary_handling.py  | 116 ++++++++-
 lbmpy_tests/test_free_slip.ipynb       | 311 +++++++++++++++++++++++++
 5 files changed, 589 insertions(+), 11 deletions(-)
 create mode 100644 lbmpy_tests/test_free_slip.ipynb

diff --git a/lbmpy/advanced_streaming/indexing.py b/lbmpy/advanced_streaming/indexing.py
index 54249a48..99b330c2 100644
--- a/lbmpy/advanced_streaming/indexing.py
+++ b/lbmpy/advanced_streaming/indexing.py
@@ -231,3 +231,31 @@ class NeighbourOffsetArrays(CustomCodeNode):
         offset_symbols = NeighbourOffsetArrays._offset_symbols(dim)
         super(NeighbourOffsetArrays, self).__init__(code, symbols_read=set(),
                                                     symbols_defined=set(offset_symbols))
+
+
+class MirroredStencilDirections(CustomCodeNode):
+
+    @staticmethod
+    def mirror_stencil(direction, mirror_axis):
+        assert mirror_axis <= len(direction), f"only {len(direction)} axis available for mirage"
+        direction = list(direction)
+        direction[mirror_axis] = -direction[mirror_axis]
+
+        return tuple(direction)
+
+    @staticmethod
+    def _mirrored_symbol(mirror_axis):
+        axis = ['x', 'y', 'z']
+        return TypedSymbol(f"{axis[mirror_axis]}_axis_mirrored_stencil_dir", create_type(np.int64))
+
+    def __init__(self, stencil, mirror_axis, dtype=np.int64):
+        offsets_dtype = create_type(dtype)
+
+        mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(mirror_axis)
+        mirrored_directions = [stencil.index(MirroredStencilDirections.mirror_stencil(direction, mirror_axis))
+                               for direction in stencil]
+        code = "\n"
+        code += _array_pattern(offsets_dtype, mirrored_stencil_symbol.name, mirrored_directions)
+
+        super(MirroredStencilDirections, self).__init__(code, symbols_read=set(),
+                                                        symbols_defined={mirrored_stencil_symbol})
diff --git a/lbmpy/boundaries/__init__.py b/lbmpy/boundaries/__init__.py
index 1d196589..6dd746dd 100644
--- a/lbmpy/boundaries/__init__.py
+++ b/lbmpy/boundaries/__init__.py
@@ -1,8 +1,8 @@
 from lbmpy.boundaries.boundaryconditions import (
     UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow,
-    ExtrapolationOutflow, NeumannByCopy, NoSlip, StreamInConstant)
+    ExtrapolationOutflow, NeumannByCopy, NoSlip, StreamInConstant, FreeSlip)
 from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
 
-__all__ = ['NoSlip', 'UBB', 'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
+__all__ = ['NoSlip', 'FreeSlip', 'UBB', 'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
            'FixedDensity', 'DiffusionDirichlet', 'NeumannByCopy',
            'LatticeBoltzmannBoundaryHandling', 'StreamInConstant']
diff --git a/lbmpy/boundaries/boundaryconditions.py b/lbmpy/boundaries/boundaryconditions.py
index d863e7b6..4a7a9e21 100644
--- a/lbmpy/boundaries/boundaryconditions.py
+++ b/lbmpy/boundaries/boundaryconditions.py
@@ -1,3 +1,5 @@
+from warnings import warn
+
 from lbmpy.advanced_streaming.utility import AccessPdfValues, Timestep
 from pystencils.simp.assignment_collection import AssignmentCollection
 from pystencils import Assignment, Field
@@ -5,7 +7,7 @@ from lbmpy.boundaries.boundaryhandling import LbmWeightInfo
 from pystencils.data_types import create_type
 from pystencils.sympyextensions import get_symmetric_part
 from lbmpy.simplificationfactory import create_simplification_strategy
-from lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays
+from lbmpy.advanced_streaming.indexing import NeighbourOffsetArrays, MirroredStencilDirections
 from pystencils.stencil import offset_to_direction_string, direction_string_to_offset, inverse_direction
 
 import sympy as sp
@@ -107,8 +109,126 @@ class NoSlip(LbBoundary):
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
         return Assignment(f_in(inv_dir[dir_symbol]), f_out(dir_symbol))
 
+
 # end class NoSlip
 
+class FreeSlip(LbBoundary):
+    """
+    Free-Slip boundary condition, which enforces a zero normal fluid velocity $u_n = 0$ but places no restrictions
+    on the tangential fluid velocity $u_t$.
+
+    Args:
+        stencil: LBM stencil which is used for the simulation
+        normal_direction: optional normal direction. If the Free slip boundary is applied to a certain side in the
+                          domain it is not necessary to calculate the normal direction since it can be stated for all
+                          boundary cells. This reduces the memory space for the index array significantly.
+        name: optional name of the boundary.
+    """
+
+    def __init__(self, stencil, normal_direction=None, name=None):
+        """Set an optional name here, to mark boundaries, for example for force evaluations"""
+        self.stencil = stencil
+
+        if normal_direction and len(normal_direction) - normal_direction.count(0) != 1:
+            raise ValueError("It is only possible to pre specify the normal direction for simple situations."
+                             "This means if the free slip boundary is applied to a straight wall or side in the "
+                             "simulation domain. A possible value for example would be (0, 1, 0) if the "
+                             "free slip boundary is applied to the northern wall. For more complex situations "
+                             "the normal direction has to be calculated for each cell. This is done when "
+                             "the normal direction is not defined for this class")
+
+        if normal_direction:
+            self.mirror_axis = normal_direction.index(*[dir for dir in normal_direction if dir != 0])
+
+        self.normal_direction = normal_direction
+        self.dim = len(stencil[0])
+
+        if name is None and normal_direction:
+            name = f"Free Slip : {offset_to_direction_string([-x for x in normal_direction])}"
+
+        super(FreeSlip, self).__init__(name)
+
+    def init_callback(self, boundary_data, **_):
+        if len(boundary_data.index_array) > 1e6:
+            warn(f"The calculation of the normal direction for each cell might take a long time, because "
+                 f"{len(boundary_data.index_array)} cells are marked as Free Slip boundary cells. Consider specifying "
+                 f" the normal direction beforehand, which is possible if it is equal for all cells (e.g. at a wall)")
+
+        dim = boundary_data.dim
+        coords = [coord for coord, _ in zip(['x', 'y', 'z'], range(dim))]
+
+        boundary_cells = set()
+
+        # get a set containing all boundary cells
+        for cell in boundary_data.index_array:
+            fluid_cell = tuple([cell[coord] for coord in coords])
+            direction = self.stencil[cell['dir']]
+            boundary_cell = tuple([i + d for i, d in zip(fluid_cell, direction)])
+            boundary_cells.add(boundary_cell)
+
+        for cell in boundary_data.index_array:
+            fluid_cell = tuple([cell[coord] for coord in coords])
+            direction = self.stencil[cell['dir']]
+            ref_direction = direction
+            normal_direction = [0] * dim
+
+            for i in range(dim):
+                sub_direction = [0] * dim
+                sub_direction[i] = direction[i]
+                test_cell = tuple([x + y for x, y in zip(fluid_cell, sub_direction)])
+
+                if test_cell in boundary_cells:
+                    normal_direction[i] = direction[i]
+                    ref_direction = MirroredStencilDirections.mirror_stencil(ref_direction, i)
+
+            ref_direction = inverse_direction(ref_direction)
+            for i, cell_name in zip(range(dim), self.additional_data):
+                cell[cell_name[0]] = -normal_direction[i]
+            cell['ref_dir'] = self.stencil.index(ref_direction)
+
+    @property
+    def additional_data(self):
+        """Used internally only. For the FreeSlip boundary the information of the normal direction for each pdf
+        direction is needed. This information is stored in the index vector."""
+        if self.normal_direction:
+            return []
+        else:
+            data_type = create_type('int32')
+            wnz = [] if self.dim == 2 else [('wnz', data_type)]
+            data = [('wnx', data_type), ('wny', data_type)] + wnz
+            return data + [('ref_dir', data_type)]
+
+    @property
+    def additional_data_init_callback(self):
+        if self.normal_direction:
+            return None
+        else:
+            return self.init_callback
+
+    def get_additional_code_nodes(self, lb_method):
+        if self.normal_direction:
+            return [MirroredStencilDirections(self.stencil, self.mirror_axis)]
+        else:
+            return []
+
+    def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
+        if self.normal_direction:
+            normal_direction = self.normal_direction
+            mirrored_stencil_symbol = MirroredStencilDirections._mirrored_symbol(self.mirror_axis)
+            mirrored_direction = inv_dir[sp.IndexedBase(mirrored_stencil_symbol, shape=(1,))[dir_symbol]]
+        else:
+            normal_direction = list()
+            for i, cell_name in zip(range(self.dim), self.additional_data):
+                normal_direction.append(index_field[0](cell_name[0]))
+            normal_direction = tuple(normal_direction)
+
+            mirrored_direction = index_field[0]('ref_dir')
+
+        return Assignment(f_in(inv_dir[dir_symbol]), f_in[normal_direction](mirrored_direction))
+
+
+# end class FreeSlip
+
 
 class UBB(LbBoundary):
     """Velocity bounce back boundary condition, enforcing specified velocity at obstacle
@@ -173,12 +293,11 @@ class UBB(LbBoundary):
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
         vel_from_idx_field = callable(self._velocity)
         vel = [index_field(f'vel_{i}') for i in range(self.dim)] if vel_from_idx_field else self._velocity
-        direction = dir_symbol
 
         assert self.dim == lb_method.dim, \
             f"Dimension of UBB ({self.dim}) does not match dimension of method ({lb_method.dim})"
 
-        neighbor_offset = NeighbourOffsetArrays.neighbour_offset(direction, lb_method.stencil)
+        neighbor_offset = NeighbourOffsetArrays.neighbour_offset(dir_symbol, lb_method.stencil)
 
         velocity = tuple(v_i.get_shifted(*neighbor_offset)
                          if isinstance(v_i, Field.Access) and not vel_from_idx_field
@@ -193,7 +312,7 @@ class UBB(LbBoundary):
         c_s_sq = sp.Rational(1, 3)
         weight_of_direction = LbmWeightInfo.weight_of_direction
         vel_term = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) * weight_of_direction(
-            direction, lb_method)
+            dir_symbol, lb_method)
 
         # Better alternative: in conserved value computation
         # rename what is currently called density to "virtual_density"
@@ -205,12 +324,13 @@ class UBB(LbBoundary):
             density_equations = cqc.output_equations_from_pdfs(pdf_field_accesses, {'density': density_symbol})
             density_symbol = lb_method.conserved_quantity_computation.defined_symbols()['density']
             result = density_equations.all_assignments
-            result += [Assignment(f_in(inv_dir[direction]),
-                                  f_out(direction) - vel_term * density_symbol)]
+            result += [Assignment(f_in(inv_dir[dir_symbol]),
+                                  f_out(dir_symbol) - vel_term * density_symbol)]
             return result
         else:
-            return [Assignment(f_in(inv_dir[direction]),
-                               f_out(direction) - vel_term)]
+            return [Assignment(f_in(inv_dir[dir_symbol]),
+                               f_out(dir_symbol) - vel_term)]
+
 
 # end class UBB
 
@@ -259,6 +379,7 @@ class SimpleExtrapolationOutflow(LbBoundary):
 
         return Assignment(f_in.center(inv_dir[dir_symbol]), f_out[tangential_offset](inv_dir[dir_symbol]))
 
+
 # end class SimpleExtrapolationOutflow
 
 
@@ -406,6 +527,7 @@ class ExtrapolationOutflow(LbBoundary):
 
         return AssignmentCollection(boundary_assignments, subexpressions=subexpressions)
 
+
 # end class ExtrapolationOutflow
 
 
@@ -459,6 +581,7 @@ class FixedDensity(LbBoundary):
         return subexpressions + [Assignment(f_in(inv_dir[dir_symbol]),
                                             2 * eq_component - f_out(dir_symbol))]
 
+
 # end class FixedDensity
 
 
@@ -493,6 +616,7 @@ class DiffusionDirichlet(LbBoundary):
         return [Assignment(f_in(inv_dir[dir_symbol]),
                            2 * w_dir * self.concentration - f_out(dir_symbol))]
 
+
 # end class DiffusionDirichlet
 
 
@@ -516,6 +640,7 @@ class NeumannByCopy(LbBoundary):
         return [Assignment(f_in(inv_dir[dir_symbol]), f_out(inv_dir[dir_symbol])),
                 Assignment(f_out[neighbour_offset](dir_symbol), f_out(dir_symbol))]
 
+
 # end class NeumannByCopy
 
 
diff --git a/lbmpy_tests/test_boundary_handling.py b/lbmpy_tests/test_boundary_handling.py
index 13f6880e..97a504bf 100644
--- a/lbmpy_tests/test_boundary_handling.py
+++ b/lbmpy_tests/test_boundary_handling.py
@@ -2,13 +2,23 @@ import numpy as np
 import pytest
 
 from lbmpy.boundaries import NoSlip, UBB, SimpleExtrapolationOutflow, ExtrapolationOutflow,\
-    FixedDensity, DiffusionDirichlet, NeumannByCopy, StreamInConstant
+    FixedDensity, DiffusionDirichlet, NeumannByCopy, StreamInConstant, FreeSlip
 from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
 from lbmpy.creationfunctions import create_lb_function, create_lb_method
 from lbmpy.geometry import add_box_boundary
 from lbmpy.lbstep import LatticeBoltzmannStep
 from lbmpy.stencils import get_stencil
 from pystencils import create_data_handling, make_slice
+from pystencils.slicing import slice_from_direction
+from pystencils.stencil import inverse_direction
+
+
+def mirror_stencil(direction, mirror_axis):
+    for i, n in enumerate(mirror_axis):
+        if n != 0:
+            direction[i] = -direction[i]
+
+    return tuple(direction)
 
 
 @pytest.mark.parametrize("target", ['cpu', 'gpu', 'opencl'])
@@ -105,6 +115,110 @@ def test_simple(target):
     assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 8] == 5))
 
 
+def test_free_slip_index_list():
+    stencil = get_stencil('D2Q9')
+    dh = create_data_handling(domain_size=(4, 4), periodicity=(False, False))
+    src = dh.add_array('src', values_per_cell=len(stencil), alignment=True)
+    dh.fill('src', 0.0, ghost_layers=True)
+    method = create_lb_method(stencil='D2Q9', method='srt', relaxation_rate=1.8)
+
+    bh = LatticeBoltzmannBoundaryHandling(method, dh, 'src', name="bh")
+
+    free_slip = FreeSlip(stencil=stencil)
+    add_box_boundary(bh, free_slip)
+
+    bh.prepare()
+    for b in dh.iterate():
+        for b_obj, idx_arr in b[bh._index_array_name].boundary_object_to_index_list.items():
+            index_array = idx_arr
+
+    # normal directions
+    normal_west = (1, 0)
+    normal_east = (-1, 0)
+    normal_south = (0, 1)
+    normal_north = (0, -1)
+
+    normal_south_west = (1, 1)
+    normal_north_west = (1, -1)
+    normal_south_east = (-1, 1)
+    normal_north_east = (-1, -1)
+
+    for cell in index_array:
+        direction = stencil[cell[2]]
+        inv_dir = inverse_direction(direction)
+
+        boundary_cell = (cell[0] + direction[0], cell[1] + direction[1])
+        normal = (cell[3], cell[4])
+        # the data is written on the inverse direction of the fluid cell near the boundary
+        # the data is read from the mirrored direction of the inverse direction where the mirror axis is the normal
+        assert cell[5] == stencil.index(mirror_stencil(list(inv_dir), normal))
+
+        if boundary_cell[0] == 0 and 0 < boundary_cell[1] < 5:
+            assert normal == normal_west
+
+        if boundary_cell[0] == 5 and 0 < boundary_cell[1] < 5:
+            assert normal == normal_east
+
+        if 0 < boundary_cell[0] < 5 and boundary_cell[1] == 0:
+            assert normal == normal_south
+
+        if 0 < boundary_cell[0] < 5 and boundary_cell[1] == 5:
+            assert normal == normal_north
+
+        if boundary_cell == (0, 0):
+            assert cell[2] == cell[5]
+            assert normal == normal_south_west
+
+        if boundary_cell == (5, 0):
+            assert cell[2] == cell[5]
+            assert normal == normal_south_east
+
+        if boundary_cell == (0, 5):
+            assert cell[2] == cell[5]
+            assert normal == normal_north_west
+
+        if boundary_cell == (5, 5):
+            assert cell[2] == cell[5]
+            assert normal == normal_north_east
+
+
+def test_free_slip_equivalence():
+    # check if Free slip BC does the same if the normal direction is specified or not
+
+    stencil = get_stencil('D2Q9')
+    dh = create_data_handling(domain_size=(4, 4), periodicity=(False, False))
+    src1 = dh.add_array('src1', values_per_cell=len(stencil), alignment=True)
+    src2 = dh.add_array('src2', values_per_cell=len(stencil), alignment=True)
+    dh.fill('src1', 0.0, ghost_layers=True)
+    dh.fill('src2', 0.0, ghost_layers=True)
+
+    shape = dh.gather_array('src1', ghost_layers=True).shape
+
+    num = 0
+    for x in range(shape[0]):
+        for y in range(shape[1]):
+            for direction in range(shape[2]):
+                dh.cpu_arrays['src1'][x, y, direction] = num
+                dh.cpu_arrays['src2'][x, y, direction] = num
+                num += 1
+
+
+    method = create_lb_method(stencil='D2Q9', method='srt', relaxation_rate=1.8)
+
+    bh1 = LatticeBoltzmannBoundaryHandling(method, dh, 'src1', name="bh1")
+    free_slip1 = FreeSlip(stencil=stencil)
+    bh1.set_boundary(free_slip1, slice_from_direction('N', dh.dim))
+
+    bh2 = LatticeBoltzmannBoundaryHandling(method, dh, 'src2', name="bh2")
+    free_slip2 = FreeSlip(stencil=stencil, normal_direction=(0, -1))
+    bh2.set_boundary(free_slip2, slice_from_direction('N', dh.dim))
+
+    bh1()
+    bh2()
+
+    assert np.array_equal(dh.cpu_arrays['src1'], dh.cpu_arrays['src2'])
+
+
 def test_exotic_boundaries():
     step = LatticeBoltzmannStep((50, 50), relaxation_rate=1.8, compressible=False, periodicity=False)
     add_box_boundary(step.boundary_handling, NeumannByCopy())
diff --git a/lbmpy_tests/test_free_slip.ipynb b/lbmpy_tests/test_free_slip.ipynb
new file mode 100644
index 00000000..4ff77fab
--- /dev/null
+++ b/lbmpy_tests/test_free_slip.ipynb
@@ -0,0 +1,311 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this test case a domain is tested with shear inflow profile at west, NoSlip walls at south, top and bottom, a FreeSlip wall north and an Outflow boundary at the eastern border.\n",
+    "\n",
+    "Due to the FreeSlip wall at the northern boundary an almost zero velocity gradient should occure near the north boundary. This is tested here."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "import pystencils as ps\n",
+    "from pystencils.slicing import slice_from_direction, make_slice\n",
+    "\n",
+    "from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling\n",
+    "from lbmpy.boundaries.boundaryconditions import FreeSlip, NoSlip, ExtrapolationOutflow, UBB\n",
+    "\n",
+    "\n",
+    "from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule\n",
+    "from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments\n",
+    "from lbmpy.stencils import get_stencil\n",
+    "\n",
+    "import lbmpy.plot as plt\n",
+    "\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gpu = False\n",
+    "target = 'cpu'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "stencil = get_stencil('D3Q27')\n",
+    "domain_size = (30, 15, 30)\n",
+    "dim = len(domain_size)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "dh = ps.create_data_handling(domain_size=domain_size, periodicity=(False, False, False))\n",
+    "\n",
+    "src = dh.add_array('src', values_per_cell=len(stencil), alignment=True)\n",
+    "dh.fill('src', 0.0, ghost_layers=True)\n",
+    "dst = dh.add_array('dst', values_per_cell=len(stencil), alignment=True)\n",
+    "dh.fill('dst', 0.0, ghost_layers=True)\n",
+    "\n",
+    "velField = dh.add_array('velField', values_per_cell=dh.dim, alignment=True)\n",
+    "dh.fill('velField', 0.0, ghost_layers=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "method = create_lb_method(stencil='D3Q27', method='srt', relaxation_rate=1.8)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Initialisation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "init = pdf_initialization_assignments(method, 1.0, (0, 0, 0), src.center_vector)\n",
+    "\n",
+    "ast_init = ps.create_kernel(init, target=dh.default_target)\n",
+    "kernel_init = ast_init.compile()\n",
+    "\n",
+    "dh.run_kernel(kernel_init)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Update Rules"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "update = create_lb_update_rule(lb_method=method,\n",
+    "                               output={'velocity': velField},\n",
+    "                               optimization={\"symbolic_field\": src,\n",
+    "                                             \"symbolic_temporary_field\": dst},\n",
+    "                               kernel_type='stream_pull_collide')\n",
+    "\n",
+    "ast_kernel = ps.create_kernel(update, target=dh.default_target)\n",
+    "kernel = ast_kernel.compile()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Boundary Handling"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def velocity_info_callback(boundary_data, **_):\n",
+    "    boundary_data['vel_1'] = 0\n",
+    "    boundary_data['vel_2'] = 0\n",
+    "    u_max = 0.1\n",
+    "    x, y = boundary_data.link_positions(0), boundary_data.link_positions(1)\n",
+    "    dist = (domain_size[1] - y) / domain_size[1]\n",
+    "    boundary_data['vel_0'] = u_max * (1 - dist)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAABTgAAAKzCAYAAAA6OWJuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAB7CAAAewgFu0HU+AABz5UlEQVR4nOzdd5hdVb0/4O/KzCSTnpAEUoiEmgACQiAKUqVZQPmBV1FUQBAQFRVEuNdy1avXCgoqIhYQxQsqKJh7kd5rCEVagnRIIJWUSZ2yfn/MmXAI00smO3nf5znPWfvstddaJ5xnA5+svVbKOQcAAAAAQBH16e0BAAAAAAB0loATAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFFZlbw+A16WU+kXETqXDeRFR34vDAQAAYMNXERGjSuVHc86renMwAJ0h4Fy/7BQR03p7EAAAAGyU9oiIB3p7EAAd5RF1AAAAAKCwzOBcv8xrKtx///0xZsyY3hwLAAAAG7hXXnklpkyZ0nQ4r7W6AOsrAef6Zc2am2PGjInNN9+8N8cCAADAxsU+EEAheUQdAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwurRgDOltGlK6bCU0rdSStemlOanlHLpdUk7rp9QVr+9r+e7MN7ne7oPAAAAAKD7VPZw+3N6uP3mzOyFPgEAAACAXtDTAWe5FyNiRkQc0oFrZkXETu2o9+8R8dFS+XcdHFdzro6Ir7ZyfnU39AEAAAAAdFFPB5zfiohpETEt5zwnpTQhIp5r78U559qIeKy1OimliojYv3S4NCL+2qmRvtGinHOr/QIAAAAAva9HA86c83/2ZPslB0XE2FL5LznnFeugTwAAAABgPbAh7KL+ibJydzyeDgAAAAAURKEDzpTS4Ig4onT4fETc3muDAQAAAADWuXW5yVBP+GBEDCiVf59zzt3U7r4ppYcjYuuIqIjG3eDvj4j/iYirO9tPSmnzNqqM7ky7AAAAALCxKnrAWf54+qXd2O6Wax1PKL0+FBF3pZQ+nHOe1Yl2X+riuAAAAACAMoUNOFNKb4mI/UqHd+ecn+6GZldHxDURcX007t6+OCKGRcSeEfHpiBgfEe+MiBtSSnvmnBd3Q58AAAAAQCcVNuCMiI9FRCqVu2v25pSc86JmPr81pfSziPhLRBwSEdtHxH9GxOkdbH98G+dHR8S0DrYJAAAAAButIgecHy+9r4qIK7qjwRbCzaZzS1NKH4qIZyNik4g4KaV0ds55dQfaf7m18yml1k4DAAAAAGsp5C7qKaUpETGpdHhNa8Fkdyo9kn556XBgROy+LvoFAAAAAJpXyIAzem5zofZ4oqw8bh33DQAAAACUKVzAmVKqioijS4dzI+If63gIeR33BwAAAAC0oHABZ0S8LyJGlMp/zDnXreP+dygrz17HfQMAAAAAZYoYcJY/nv67ddlxSmlovD57dHlEPLAu+wcAAAAA3qhQAWdKaZNonMEZEfFozvnhDlx7a0opl14Tmjn/7pRS/1auHxQRf4rXZ4/+Jue8qt2DBwAAAAC6XWVPNp5S2jsitin7aGRZeZuU0nHl9XPOl7TR5NER0bdU7u7Zm2dHxGUppasi4s6IeCYiaiJiaETsFRGnRMRbSnVnRsQ3url/AAAAAKCDejTgjIgTI+LYFs69s/Qqd0kb7TU9nl4fEZd1flgt2iQax3xiK3Vui4hjcs4Le6B/AAAAAKADejrg7DYppW0j4u2lwxtyzq92cxdfiogDI2LPiJgYjbNNh0XjWpuzI+K+iPifiLg+52wndQAAAABYD/RowJlzPi4ijuumtv4VEakL1+/fxvkHwqZBAAAAAFAohdpkCAAAAACgnIATAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAAqrsrcHQPMuOfuuGD5oVG8PAwAAgA3YazXzensIAF1mBicAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMKq7O0B0Ly9Vh4doyvWXX/vHTN13XUGAABAiwZvf/Y666t2YW3EZeusO4AeYQYnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIXVowFnSmnTlNJhKaVvpZSuTSnNTynl0uuSdrZxXNk1bb2O66ZxjyyN+Z8ppSWl1z9Ln43ojj4AAAAAgK6r7OH25/Rw+90upfT2iPhbRIxe69ROpdeJKaUjcs73r+uxAQAAAABv1NMBZ7kXI2JGRBzShTYOjYjZrZx/uQttR0ppfET8PSJGRURdRJwbEVNLpw+LiNMjYkxE/D2lNDnn3KX+AAAAAICu6emA81sRMS0ipuWc56SUJkTEc11o76mc8/PdMbAWfCcaw82IiI/mnP9cdu6OlNL0iLgiIjaNiG9HxHE9OBYAAAAAoA09ugZnzvk/c85Tc87r/aPqKaXREXFM6fC6tcLNiIjIOf8pIq4rHX68dA0AAAAA0Evsov6698frfx4Xt1LvktJ7n9I1AAAAAEAvEXC+bu+y8m2t1Cs/984eGgsAAAAA0A7rcpOh7nBxSmliRIyMiCUR8XRE3BgRv8g5z+pi2zuU3hfnnF9tqVLO+ZWU0pKIGBIR23ekg5TS5m1U8cg7AAAAAHRA0QLO/cvKI0qvt0fEGSmlL+Scf9mFtpvCx/bsjP5SROwYEeM72MdLHawPAAAAALSiKAHnsxFxVUTcE6+HhFtFxFER8cGIqI6IC1NKOed8USf7GFx6r2lH3WWl90Gd7AsAAAAA6AZFCDj/GhG/yznntT6fFhFXpJQOi8bwsyoifpxSuqa1R8xbUV16X92OuqtK7/072EdbMz5HR+P3AgAAAADaYb3fZCjnvLiZcLP8/NSI+FbpcEBEnNDJrlaW3vu2o26/0vuKjnSQc365tVdEdCaYBQAAAICN1nofcLbTRRHRFILu18k2lpbe2/PY+cDSe3seZwcAAAAAesgGEXDmnOdGxILS4bhONtO0uVBbO51HvP6ouU2DAAAAAKAXbRABZ0mLj7G30xOl96EppdEtVUopjYmIIaXDJ7vYJwAAAADQBRtEwJlSGhURI0uHszvZzJ1l5dYecy8/d1cn+wIAAAAAusEGEXBGxEkRkUrl2zrZxjUR0VAqH99KveNK7w2lawAAAACAXrJeB5wppQkppV3bqHNYRHy9dLgiIi5uod6tKaVcek1Y+3zO+dWIuKx0eGhK6YPNtPFvEXFo6fD3pWsAAAAAgF5S2ZONp5T2johtyj4aWVbeJqV0XHn9nPMlazUxISJuSSndExF/j4hHImJu6dxWEfHB0qtp9uaXcs6zujDkr0TEuyNiVET8T0pp94iYWjp3WEScUSrPi4ivdqEfAAAAAKAb9GjAGREnRsSxLZx7Z+lV7pIW6u5ZerVkeUR8Med8UYdGt5ac80sppcMj4m8RMToiziq9yr0aEUfknF8OAAAAAKBX9XTA2VXTI+Jj0Rhu7h4RY6JxFmhlRLwWEY9HxE0R8euc89yWGumInPN9KaWdIuLzEXFENM4ijYh4LiKujoif5JwXdEdfAAAAAEDX9GjAmXM+Ll7flKcz1y+NxnUxL2urbjva2r8DdedHxNdKLwAAAABgPbVebzIEAAAAANAaAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFg9GnCmlDZNKR2WUvpWSunalNL8lFIuvS5pZxsDUkpHppR+kVKallJ6LaVUm1JakFK6J6X0jZTS6G4a7/Nl42vt9Xx39AcAAAAAdE1lD7c/pysXp5R2joi7ImJQM6c3iYh3lF5fTCmdlHO+oiv9AQAAAADF0tMBZ7kXI2JGRBzSgWuGxOvh5l0RMTUiHoiIBRExKiKOjIhPlepdllJaknO+thvGenVEfLWV86u7oQ8AAAAAoIt6OuD8VkRMi4hpOec5KaUJEfFcB65viIg/RcQ3c85PNHP++pTStRHx14ioiIifppS2zTnnLo57Uc75sS62AQAAAAD0sB4NOHPO/9nF6++OiLvbqHN1SumqiDgqIraOiF0j4sGu9AsAAAAAFMOGsov6LWXlrXttFAAAAADAOrWhBJz9ysr1vTYKAAAAAGCd2lACzv3Kyk92Q3v7ppQeTiktTSktTyk9l1K6IqV0REopdUP7AAAAAEA3WJe7qPeIlNIuEfG+0uGjOefuCDi3XOt4Qun1oYi4K6X04ZzzrI42mlLavI0qozvaJgAAAABszAodcKaU+kXEr6NxB/WIiK90scnVEXFNRFwfEY9FxOKIGBYRe0bEpyNifES8MyJuSCntmXNe3MH2X+ri+AAAAACAMoUOOCPiZxGxe6n8u5zz37vY3pSc86JmPr81pfSziPhLRBwSEdtHxH9GxOld7A8AAAAA6ILCBpwppX+PiBNLh9Mi4jNdbbOFcLPp3NKU0oci4tmI2CQiTkopnZ1zXt2BLsa3cX50NH4XAAAAAKAdChlwppROjoj/Lh3OiIj35pyX9XS/OefFKaXLI+LUiBgYjbNH7+7A9S+3dt7+RQAAAADQMYXbRT2l9JGIuKB0+EJEHJxznr8Oh/BEWXncOuwXAAAAAFhLoQLOlNL7I+LSaBz3KxFxYFuzIntAXsf9AQAAAAAtKEzAmVI6MCL+FI2P1S+Ixpmbz/TCUHYoK8/uhf4BAAAAgJJCBJwppb0i4uqI6BcRiyPi0Jzz470wjqERcXTpcHlEPLCuxwAAAAAAvG69DzhTSm+LiP+Nxk19lkXE+3LO0zvRzq0ppVx6TWjm/LtTSv1buX5QNM4gHVH66Dc551UdHQcAAAAA0H16dBf1lNLeEbFN2Ucjy8rbpJSOK6+fc75kreu3jojrImJY6aOvRsTilNJbW+l2bs55bieGe3ZEXJZSuioi7oyIZyKiJiKGRsReEXFKRLylVHdmRHyjE30AAAAAAN2oRwPOiDgxIo5t4dw7S69yl6x1vE9EbFp2/ON29PnN6Hz4uEk0jvnEVurcFhHH5JwXdrIPAAAAAKCb9HTAWSRfiogDI2LPiJgYjbNNh0XjWpuzI+K+iPifiLg+52wndQAAAABYD/RowJlzPi4ijuvC9ZfEm2d1drat/ds4/0DYNAgAAAAACmW932QIAAAAAKAlAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYfVowJlS2jSldFhK6VsppWtTSvNTSrn0uqQT7b0npfTXlNLLKaVVpfe/ppTe083jHpBS+nJKaVpKaWFKaVlKaUZK6ZyU0hbd2RcAAAAA0HmVPdz+nO5oJKXUJyIuiogT1jo1rvQ6IqX064g4Oefc0MW+tomI/4uIbdc6NbH0OjGldEzOeWpX+gEAAAAAum5dPqL+YkRc38lrvxOvh5sPRcRHImJK6f2h0ucnRsS3uzLAlNLgiPjfeD3c/FVEHBgRe0XEVyKiJiKGRMQVKaW3daUvAAAAAKDrenoG57ciYlpETMs5z0kpTYiI5zrSQEppu4j4UunwgYjYN+e8onQ8LaV0TUTcFhG7R8SZKaXf5pyf7uR4z4yI7UrlL+ecf1h27p6U0q2lvgZExE8iYv9O9gMAAAAAdIMencGZc/7PnPPUnHNXHlX/QrwexH6uLNxs6mN5RHyudFgZEV/sTCcppaqIOK10+GREnLN2nZzz3RHxm9LhfimlPTrTFwAAAADQPdbrXdRTSikiPlA6nJFzvre5eqXPZ5YOP1C6rqMOiIihpfLvWlnL85Ky8v/rRD8AAAAAQDdZrwPOiNgyIsaWyre1Ubfp/LiImNCJvvZupq3mPBARy0vld3aiHwAAAACgm/T0GpxdtUNZeUYbdcvPbx8dXOuzvX3lnOtSSk9HxM6lftotpbR5G1VGd6Q9AAAAANjYre8BZ3kg+HIbdV8qK4/vQl/Lcs6L2tHXzhExKqXUL+e8qp19vNR2FQAAAACgvdb3R9QHl5Vr2qi7rKw8qAt9tdVPd/QFAAAAAHSD9X0GZ3VZeXUbdctnUfbvQl9t9dOVvtqaWTo6IqZ1oD0AAAAA2Kit7wHnyrJy3zbq9isrr+hCX2310+m+cs6tPmbfuc3fAQAAAGDjtb4/or60rNzWo+ADy8rtecy8pb7a88h5V/sCAAAAALrB+h5wls94bGsH8vLHvzuzmU9TXwNTSsPa2de8DmwwBAAAAAB0s/U94HyirDypjbrl55/sqb5SSpURsXUX+gEAAAAAusn6HnA+FxGzS+X92qi7b+l9VkQ834m+7iwrt9bX7vH6I+p3daIfAAAAAKCbrNcBZ845R8TVpcNJKaV3NFev9HnTrMurS9d11K0RsbhUPja1vOPPcWXlv3aiHwAAAACgm6zXAWfJTyKivlT+aUqpf/nJ0vFPS4d1pfpvklK6JKWUS6/91z6fc14dEeeXDrePiC8108aeEXFC6fC2nPO0jnwRAAAAAKB7VfZk4ymlvSNim7KPRpaVt0kpHVdeP+d8ydpt5JyfSin9MCLOjsbHw+9KKX0/Ip6JxrUwz4qIXUvVf5hz/lcXhvzDiPhwRGwXET9IKW0TEZdHxIqIOCAi/iMa/8xWRMQXutAPAAAAANANejTgjIgTI+LYFs69s/Qqd0kLdb8SEZtGxCejMcy8vJk6v4mIr3Z8iK/LOS9NKb0vIv4vIraNiJNKr3JLIuKYnPPDXekLAAAAAOi6IjyiHjnnhpzzCRHxvmhck3N2RKwuvV8dEe/NOZ+Yc27ohr6ejsYQ9ayIeCAiFkXE8oiYGRE/joidc85Tu9oPAAAAANB1PTqDM+d8XLxxU56utvd/0Ti7skfHknNeFhE/KL0AAAAAgPVUIWZwAgAAAAA0R8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFgCTgAAAACgsCp7ewAAAADAhmf69OnjI2KviHhHRIyPiAG9OyKgF9VExGMRcWNE3D958uS67mxcwAkAAAB0m+nTp6eIODEiTo6I1KdPn/59+vSpDk+Rwsasvr6+flLO+f0R8fD06dNPmzx58vLualzACQAAAHSn/4yIwyorK0dWVlaOiIg+ffr0aejTp09DROReHhuwjuWcU0NDQ0VlZWVqaGhYtnr16j4Rcf706dNP6a6ZnAJOAAAAoFtMnz59m4g4rKqqanRFRcXwTTbZZMHQoUMX9+/ff1VKqbeHB/SShoaGtHjx4sFz5swZHRGbr169uiEi9oiIe7qjfdPDAQAAgO5yZERUVFRUDBs1atTcMWPGzB0wYIBwEzZyffr0ycOHD1+y2WabvdqnT5+BKaW+EXFQt7XfXQ0BAAAAG71dKyoqBkdE2mSTTV7r7cEA65ehQ4cuTSnlioqKgRGxU3e1u94HnCmlW1NKuYOv/TvRzzd6sn0AAADYCAxJKVVWVFTUVVZWNvT2YID1S58+fXKfPn3qI6IiIgZ1W7vd1dB6pCEi/tXbgwAAAICNUIqISCnZTAhoVk/cH4qwydDxETGwjTo7RMQVpfJNOedZXeyzrSmyz3WxfQAAAACgG6z3AWfOuc0wMaX08bLDS7uhz8e62gYAAAAA0PMK/4h6SqlPRBxTOqyJiKt6cTgAAAAAwDpU+IAzIg6MiHGl8l9yzst7czAAAAAAwLqzIQScnygrd/nxdAAAAIANyfnnnz8ipTQ5pTR55syZfXt7PL1typQpE1NKk6dMmTKxt8dC91jv1+BsTUppUET8v9LhCxFxaze1e31EvC0ihkXEooh4IiL+ERG/zDm/1oV2N2+jyujOtg0AAABFcvC5t22/cPnqqt4eR0/bZEDf2htO3+/J3h4HbMgKHXBGxFHx+g7rf8g5d9c28weXlUdFxH6l11kppeNyzld3st2XujwyAAAA2AAsXL66akHNhh9wAj2v6AFndz+e/mhE/C0i7o+I2RFRFRETo3ETo0OicUbnlSmlw3PO13ZDfwAAAACsQ/fff//M3h4D3auwAWfpce/9S4f35pyf6mKTP8k5f6OZz++LiEtTSidHxIURURERv04pbZ1zXtnBPsa3cX50REzrYJsAAAAAsNEqbMAZER+L1zdJ+l1XG8s5L2rj/C9TSntExAkRMTYaH4+/rIN9vNza+ZRSR5oDAAAAgI1ekXdR/3jpfVVEXLGO+vxlWXm/ddQnAAAAQIvmzZtXceqpp47bcsstd6yurt5tk0022WWvvfba7re//e3w9rYxc+bMvieccML4bbbZZseBAwfu2r9//1232GKLt370ox/d4v777+/f2rVNO7SffvrpYyMi/v73vw8+6KCDtt500013rq6u3m2rrbba8cwzzxyzZMmSN+RQV1xxxdD99ttvm6Z6W2+99Y7//u//PnrlypUtzgBbuXJl+uMf/zj0E5/4xFve+ta3bj9kyJC3VVZW7jZs2LC37bzzzpNOP/30sa+88kqrE/pa20V95syZfZu+z/nnnz8iIuKvf/3rkHe9613bjBw5cpe+ffvuNm7cuJ2OOeaYtzzzzDPWkF1PFHIGZ0pp94jYoXQ4tSs7m3fQE2XlceuoTwAAAIBmPfjgg9Xvfve7t5s3b96asG3VqlWV99xzz+B77rln8LXXXrtgn332WdpaGz/72c9GnHHGGVusXr36DcHiiy++2O/FF1/s96c//WnkmWeeOeu73/3uq22N5z/+4z9Gf+973xtXvg/0c889V/2jH/1o7E033TT01ltvfWrQoEENJ5xwwvhLLrlk0/Jrn3322ervfe974+66667BN998878qK98cWx1zzDFbXHXVVSPW/nzx4sUVjz766MBHH3104MUXXzzqiiuuePqQQw5Z1tZ42/KZz3xm3AUXXDC6/LPZs2f3/eMf/zjq2muvHX7jjTfO3G233Tq6hCHdrJABZ7xxc6EuP57eAd21SzsAAABAlyxcuLDP+973vm2bws33ve99r33iE5+YP2bMmLonnnii+qc//elmf/nLX0bMmDGjxRmYl19++dDTTjttQs45BgwY0HDyySfPOfTQQ5dUVlbmO+64Y9B55503etGiRZXf+973xg0bNqz+rLPOmtdSWzfeeOOQRx99dODb3va2ZaeccsrcHXbYYeXcuXMrzz///E1vv/32oQ899NDAr371q6M32WST+ksuuWTTfffdd/EnP/nJ+VtvvfXqF154oeqHP/zhmEceeWTgHXfcMeTcc88d9eUvf/lNfdXV1aXNN9981Xve855FU6ZMWbbllluurqqqys8++2zfG264Ycif//znkYsWLao8+uijt3n00UcfHzduXF1n/3wvueSSUQ899NDAPfbYo+aEE06Yt8MOO6xcuHBhxe9+97sRf/3rX0e89tprlZ/85CcnPPzwwzM62wfdI5Un6kWQUqqKiFkRMSoi5kXE2Jxzp3+sHex793h9E6Bf55w/1c3tbx4RL0VE3HRUxOiB3dl66947Zuq66wwAAIAWDd7+7HXWV+3C2ph5+poNpce3tXdEW6ZPn/5/lZWVO1RXVw+bOHHi063VnfztG3ZeULN6g3/Ed8SgvrXTv3rwP3ui7ZNPPnnziy66aLOIiLPPPvtNMyxXrVqVDjzwwG3uuuuuIU2fzZgx49GJEyeubjo/fvz4nebNm1c1YMCAhhtuuGHGXnvttaK8jaeeeqrv3nvvPWnevHlV1dXVDc8+++yjY8aMeUMOk1Ka3FQ+9NBDX5s6deqz5bMv6+rqYvfdd5/0yCOPDBw4cGBDXV1dOuaYY+b95je/eam8naVLl/aZNGnSjrNnz+673XbbrZg5c2b5k7QREfH444/323777Vf16dP8qov3339//wMOOGDS8uXL+5x22mmvnHfeebPXrjNlypSJ06ZNG7THHnvUrL2j+syZM/tOmjRpp6bjo48+ev5ll132wtr9HX300VtcccUVIyMi7rzzzife+c53vuHPjZbNnDlzm5UrVy6qq6t7YvLkye/tjjaLuAbne6Ix3IyI+OO6CjdLTi4r37YO+wUAAABYY+XKlenyyy8fGRGx3XbbrfjOd77zpsfH+/Xrl3/3u989X1lZ2ezstt///vfDmmZ/fuELX3hl7XCz1Pbqb37zmy+X+uxzwQUXvOnx8CbV1dUNv/vd715Y+9HyysrKOP744+dFRCxbtqzP8OHDay+44II3hemDBw9u+NCHPrQgIuKpp57qv2DBgoq16+y4444thpsREVOmTFlx9NFHz4+IuPbaa4e1WLEdRo0aVfvb3/72xeb6+/d///c1f9633HLL4K70Q9cVMeAsfzz90vZckFI6LqWUS69vNHN+p5TSNm20cVJEnFg6fDUi/trO8QIAAAB0qzvvvHPAkiVLKiIiPvKRjyxoKfTbeuuta/fee+8lzZ276aabhkREpJTiM5/5zPyW+jr++ONfGzRoUH1ExC233DKkpXrvfOc7l2y22Wb1zZ2bPHnymvD0Pe95z6J+/fo1G7q+7W1vW95UnjlzZt+W+moyb968iscff7zfAw88UD1t2rTqadOmVQ8bNqwuIuKZZ57pv2rVqhY3LGrLe9/73tf69+/f7Dh32WWXVQMGDGiIiHj22Wf7dbYPukeh1uBMKQ2PiMNKh4/lnB/spqYnR8SvU0q3RMS1EfFoRCyIxj+fSRFxTEQcUqpbHxEn5Zy7vFAtAAAAQGc8/PDDa9bVfMc73tFqRjF58uRlt95669C1P29am3PcuHGrxo4d2+ITstXV1XmHHXZYfv/99w9+6qmnWlzPc5tttlnV0rlNNtlkTfvbbbddi5vyDB8+fE1A2hTgru3+++/v/8Mf/nCzW2+9dcj8+fNbXOagoaEh5s+fX9HZdTgnTZrU6uZBQ4YMqVu+fHnfmpqaIk4g3KAUKuCMiA9HRFMq3q7Zmx1QEREHlV4tWRARJ+Sc/97NfQMAAAC028KFC9dkOmPGjKltre5mm23W7PlFixZVRESMGDGizQBw0003rY2IWLJkSYtZUtOMxuZUVFS0q175TNS6uro3zb788Y9/PPLMM898S319fbtmZi5btqzT4WNr44x4faztHQs9p2gJ88dL7/URcVk3tvt/EXFCRPw6IqZHxMsRsSIiVkbE7Gic1fn5iNgq53x1N/YLAAAA0CUpdS1f6+r168pDDz1U3RRubrLJJnVf+9rXXr7jjjuefPXVVx9euXLlgznn6Tnn6T/+8Y+fb7qmaJtr0zmFmsGZc35nJ6+7JCIuaeX83Ij4bekFAAAAsF4bPnz4mlmXs2fPrtp5551bfDx8zpw5zT7GPWzYsPqIiPnz57eZD82dO7cqovGx7I6Ptnv86le/GlFfX58qKirixhtvnLnrrrs2+wh5+exWNg5Fm8EJAAAAsNF729vetmbTnnvvvXdga3UffPDBZs9PmjRpRUTErFmz+s2ePbvFUHDVqlXpiSeeGBDRuGN750bcdU1rhk6cOHF5S+FmRMvflw2XgBMAAACgYPbee+/lQ4YMqY+IuOKKK0Y0NDS/XORzzz1Xdeeddza78/mBBx64JKLxMe4LLrhgREt9XXLJJcNramoqIiIOOOCAZndkXxea1uRcvnx5i3nWCy+8UHXTTTe9aUMlNmwCTgAAAICC6d+/f/7Qhz40P6JxZuPXv/71zdauU1tbG8cdd9wWtbW1zS6y+bGPfWzRqFGjaiMizjvvvDH333//m3ZIf/rpp6u+9rWvbR4RUV1d3XDqqacu6N5v0n5bbbXVyoiIF198sfqGG2540yzNpUuX9vnQhz605cqVK+VdGxn/wAEAAAAK6Hvf+94rTTukf+c739n88MMP3/Ivf/nLkDvvvHPARRddNHy33Xbb/vbbbx+64447Lm/u+urq6nz++ee/kFKKmpqaine9612TzjzzzDE33HDDwJtvvnngN7/5zU3f/va37zBv3ryqiIhvfOMbL48ZM6bX1uA8/vjjF0RENDQ0xFFHHbXt2WefPfraa68ddMsttwz4/ve/P2qnnXba4f777x+822671fTWGOkdFl0FAAAAKKARI0bUT5069an3vOc9282fP79q6tSpm0ydOnWT8jpHHXXUgn333Xfp5z//+QnNtXH00Ucvnjdv3vNf+tKXtli2bFmfH/3oR2N/9KMfjS2vU1FREWeeeeass846a14Pfp027bfffsvPOOOM2eecc87YpUuXVnz/+98f9/3vf/8NdT71qU/Neetb37riwQcfHNRLw6QXmMEJAAAAUFC77777yscee+zxU0455dUttthiVd++ffOwYcPq3v72ty+98MILn/vLX/7yfFttfO5zn1vwz3/+87Hjjz9+7lZbbbWyf//+DdXV1Q3jx49fdfTRR8+/++67n/jud7/76jr4Om360Y9+9Mrll1/+9Dvf+c4lQ4YMqa+qqsqbbbZZ7SGHHLLoqquu+tdFF130cm+PkXUv5Zx7ewyUpJQ2j4iXIiJuOipi9Drc8+u9Y6auu84AAABo0eDtz15nfdUurI2Zp89sOhyfc+5SODR9+vT/q6ys3KG6unrYxIkTn26t7uRv37DzgprVVV3prwhGDOpbO/2rB/+zt8cB64uZM2dus3LlykV1dXVPTJ48+b3d0aZH1AEAAIB1bpMBfWt7ewzrwsbyPaE3CTgBAACAde6G0/d7srfHAGwYrMEJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACisyt4eAAAAALDxeeY/3rp93dJ5Vb09jp5WOXhU7db//diTPdH2UUcdNeGqq64aERExY8aMRydOnLi6rWvGjRu30+zZs/uOHTt29axZsx4tP5dSmtzSdf369cvDhg2r22GHHZYfeeSRr5188skLqqpa/sfX1E9z56qqqvLQoUPrtt122xWHH374os9+9rMLBg8e3NDW2KElAk4AAABgnatbOq+qfsncDT7g3FCsWrUqzZkzp2rOnDlDb7nllqG/+MUvNv3HP/7xr/Hjx9d1tK3a2to0f/78qvnz51fdc889Qy644ILNpk6d+q9ddtllVU+MnQ2fgBMAAACANXbcccflF1988XPlny1ZsqTikUce6X/RRRdtOnPmzP5PPPHEgCOOOGLr6dOnz2ytrVGjRtX+7//+71Plny1fvrzP448/Xn3xxRePeuCBBwa9/PLL/Q4//PBtZ86c+Xj//v1zT3wnNmwCTgAAAADWGDBgQMMee+yxcu3PDzzwwGUnnXTSgp133nmH5557rvrBBx8cdOONNw486KCDlrXUVlVVVW6urf3222/5ySefvHCvvfba7v777x/80ksv9bv00kuHn3zyyQu7+/uw4bPJEAAAAADtMmjQoHziiSfObTq+9957B3a2rYqKivjCF74wp+l42rRpA7o6PjZOAk4AAAAA2m2rrbZas5nRqlWrUlfa2nbbbdesu7lq1So5FZ3ihwMAAABAuz3//PNrdkd/y1ve0ubO7a155pln+nVXW2y8BJwAAAAAtEtNTU369a9/vWlERP/+/RsOP/zwJZ1tq6GhIc4777zNIiJSSnHkkUcu6qZhspGxyRAAAAAAayxfvrzPtGnTqss/q6mp6fPwww8P+NWvfjXqmWeeqU4pxde//vWXR48eXd9aW7W1tWnttlasWNHniSeeqL7kkktG3nfffYMjIk455ZRXJ0+e/KbNiKA9BJwAAAAArPH4448PmDJlyo4tnX/nO9+55Kyzznr18MMPX9pWW/Pmzatqra2ddtpp2Ze+9KVXP/GJTyzq5HDBI+oAAAAAtN999903+Kc//emmTz/9dFVX23r88ccHXnjhhaPuv//+/t0xNjZOAk4AAAAA1thjjz1qcs7Ty18rV658cMaMGY9+97vffXHQoEH1N9xww7C99tpr+4ceeqi6tbbGjh27eu22Vq9ePf3ZZ5/9589+9rPnRo8evfqee+4ZcuCBB0687rrrBq2r78iGRcAJAAAAUEAppTXlnHO7rmlvvbX169cvT5w4cfXZZ58974YbbphZWVmZ582bV3XiiSdu0dG2qqqqYsstt6z9zGc+s/Duu++eMXTo0PqampqKT37yk1vW1tZ2anxs3AScAAAAAAVUXV3d0FRetmxZuzKeFStW9ImIGDBgQENbdVuy++67r9xvv/0WR0Q8+OCDg/75z3/262xbW2yxRe2RRx65ICJi9uzZff/+978P6WxbbLwEnAAAAAAFtMkmm9Q1lWfNmtXmepgrVqxIS5curYyIGDp0aF1b9Vuz3Xbbrdnx/MEHH+zS+pmTJk1a09YjjzxiLU46TMAJAAAAUEA777zziqbytGnTBrRV/9577+1fX18fERE77LDDijaqt6qurm7N8/G1tbWptbrtaKvZdqG9BJwAAAAABXTooYcuraioyBERV1111SYNDa0/df673/1uRFP5oIMOWtKVvh9++OGBTeUJEyas7kpb06dPX9PW+PHju9QWGycBJwAAAEABjR8/vu4973nPaxERTzzxxICvfOUro1uqe8011wy+7LLLRkU07mz+0Y9+dFFn+7388suHTps2bVBExLBhw+r233//ZZ1t68477xwwderUTSIiqqqq8uGHH96l4JWNU2VvDwAAAACAzvnZz3728t133z1k4cKFld/73vfG3XHHHYM/8pGPLNh+++1XVVVV5RdeeKHv3//+96FXXnnliPr6+tSnT5+48MILn6+sbDkSWr58eZ9p06ZVl3+2evXq9OKLL/adOnXq0CuuuGJk0+df+9rXZlVVtbz8Z21tbVq7rbq6uvTKK69UXX/99UN+97vfjVq9enWKiDjllFNeHTduXJfWBmXjJOAEAAAAKKgtttii9uabb55x5JFHbvPss89W33XXXUPuuuuuZnciHzx4cP1FF1303OGHH760tTYff/zxAVOmTNmxtTqVlZX5rLPOmnX66afPb63evHnzqtpqK6UUxx133Nyf/OQns1urBy0RcAIAAAAU2C677LJqxowZj//617/e5Jprrhn2z3/+c+DChQsr6+vr09ChQ+u23XbblQcddNDiz3/+8/NHjBhR35k+KioqYtCgQfVvectbVu69995LTz311PlvfetbV3WmrT59+sSAAQPqx48fv2qPPfaoOfnkk+fvtddeXdr0iI2bgBMAAACg4KqqquLTn/70wk9/+tMLO9tGznl6d41n1qxZj3ZXW9AWAScAAACwzlUOHlXb22NYFzaW7wm9qRABZ0opt7PqbTnn/buhv49ExPERsXNEDIuIORFxR0T8POd8T1fbBwAAgI3d1v/92JO9PQZgw9CntwewPkkp9U8p/W9E/DEiDo6IzSKiX0S8JSKOiYg7U0r/2YtDBAAAAADKFGIGZ5lfRMQFrZxf1sX2fxsR7y2Vb4mI8yJidkTsFBH/ERFbR8Q3Ukqv5Jwv6mJfAAAAAEAXFS3gnJtzfqwnGk4pvSsiji4d/j0i/l/OuWlnsWkppWsiYno0zub8fkrpzznn13piLAAAAABA+3hE/XVfKr3XRcSpZeFmRETknOdHxFmlw2ERceK6GxoAAAAA0BwBZ0SklAZHxIGlwxtzzi+3UPWqiFhSKv+/Hh8YAAAAANAqAWejPSKib6l8W0uVcs6rI+LepmtSSlU9PTAAAAAAoGVFCzj/LaX0REppeUppaUrpXyml36WUDuhiuzuUlWe0UbfpfGVEbNvFfgEAAACALijaJkM7rHW8Ten1iZTS3yLiuJzz4k60u3lZuaXH05u8VFYeHxFPtLeTlNLmbVQZ3d62AAAAAIDiBJzLI+KaiLgpGmdQ1kTEqIjYLyJOiYgREXFERFydUjo451zbwfYHl5Vr2qi7rKw8qIP9vNR2FQAAAACgvYoScI7LOS9q5vMbUko/jYhrI2LXaAw8Px0R53ew/eqy8uo26q4qK/fvYD8AAAAAQDcqRMDZQrjZdG5OSumD0TizsyoiPhcdDzhXlpX7tlirUb+y8ooO9jO+jfOjI2JaB9sEAAAAgI1WIQLOtuScn00p3RAR742IbVJKY3POszvQxNKycluPnQ8sK7f1OPsb5JxbXd8zpdSR5gAAAABgo1e0XdRbU77Zz7gOXlsePLa1EVD5LExragIAAABAL9qQAs7chWvLw9FJbdRtOl8XEf/qQp8AAAAAQBdtSAHnDmXljjyeHtG47mXT5kL7tVQppdQ3It7RdE0ndmsHAAAAALrRBhFwppS2jIiDS4fP5JxndeT6nPPSiLipdHhQSqmlx9SPjIghpfJfOzxQAAAAAKBbrfcBZ0rp8JRSi5shpZQ2i4gr4/Xdzy9ops5xKaVcen2jhaZ+VHqvjIifp5Qq1mpjZER8v3S4KCJ+3e4vAQAAAAD0iPU+4IyIn0bECyml81NKH0kp7ZlSeltK6aCU0rcj4rGI2LVU986I+HlnOsk53xwRl5cO3x8RN6SU3p9S2j2ldHxE3BsRbymdPyvn/FqnvxEAAAAAb1BTU5POOOOMMRMnTtyhf//+u6aUJqeUJn/yk58cHxFx/vnnj2j6bObMmX3bao+NR4szI9czYyPic6VXS66MiBNzzqu60M8no/ER9PdGxAGlV7mGiPivnPNFXegDAAAANnpHXH3E9q+tfK2qt8fR04ZXD6/92wf+9uS66GvBggUVF1100Sb/+Mc/hj799NP9Fy5cWFlZWZlHjBhRt/POOy87/PDDFx1//PGvVVauf3HQqlWr0r777jvxoYceGtjbY1lfzZw5s++kSZN26sg1Bx544KIbb7zxmZ4a0/pi/ftFv9mx0bjxz54RsVVEjIzGELImIl6KiLsj4nc553u62lHOeUVEvC+l9NGIOC4idomIYRExJyLuiIifdUc/AAAAsLF7beVrVQtXLtzgA8515Zxzzhn57W9/e9yiRYvelPXU1NRUvPDCC/3+/ve/b/L9739/5c9//vMXDj300JqeHlNKaXJExBe/+MVXzj333FY3hL744ouHN4WbRx111ILjjjtuwWabbVYbEbHZZpvV9fRYKbb1PuDMOd8WEbd1sY1LIuKSDtT/Y0T8sSt9AgAAAKwLJ5100ua/+tWvNouIqKioyIcddthr73//+xdtueWWq1avXt3niSeeqP7Tn/60yb333jv4mWeeqX7/+9+/3YUXXvjc8ccfv94sv3fTTTcNiYgYOXJk7eWXX/78+jjLdH1y4IEHLvrud7/b5ibbw4YNa1gX4+ltfi0AAAAABfW9731vVFO4udlmm9VeddVV/9prr71WlNc59NBDa774xS/Ov/DCCzc57bTTJqxevTqdcsopW06cOHHl2nV7yyuvvFIVETF+/PjVws22DR06tH6PPfZY2dvjWF8UYZMhAAAAANby1FNP9f3GN74xPiKif//+Dddff/3M1gLLU045ZeF55533fETE6tWr03HHHbdlQ8P6McFv9erVfSIiKisrc2+PheIRcAIAAAAU0Pe///3NVq1alSIiTj/99Nk777xzmxsvf/rTn164zz77LImI+Ne//tX/8ssvH7p2nXHjxu2UUpp81FFHTWitraOOOmpCSmnyuHHj3rDxTdP1Tcc//vGPxzTtft70OuqooybMnDmzb9PxtGnTBkVETJs2bVB5vbXbbo/Zs2dXnnbaaWO33377HQYPHvy2fv367TZu3LidjjjiiC2vu+66QS1d09TnD37wg1Gtfd/ynd3X9h//8R+jU0qTKysrd1u8eLHcbR3xBw0AAABQMA0NDXHllVeOiIiorq5u+OIXvzi/vdd+/vOfn9NUvuSSS0b2xPh6y1VXXTVk4sSJb/3pT386ZsaMGf1ramoqVq9enWbPnt336quv3uTd7373xE984hNvqa+vf8N1Y8eOrdt6661XRkTcfvvtg5tr+957713z+d13391snTvvvHNwRMSOO+64fOjQoWumxzaFvuXBL93HogYAAAAABTN9+vTqxYsXV0RETJ48uWbEiBH1bV3T5AMf+MCS6urqhpUrV/Z54IEHmp3R2BXXXnvtU6tWrUpTpkzZMSLimGOOmff5z39+bnmdkSNH1o8dO7bu/vvvfzwi4vjjj9/y8ccfH7Djjjsuv/jii59rqtevX792P7J+99139z/66KO3qa2tTZWVlfnYY4+de8QRRywePHhw/bRp0wb85Cc/GTNr1qy+v//970cNHDiw/he/+MUbNunZc889lz7zzDPV991335v+TJ566qm+s2fP7lt23H/OnDkVm2222Zo/99ra2njwwQcHRUTstddeS9s7brpOwAkAAABQMA888MCApvIuu+yyvCPXVlZWxsSJE1c88sgjA1977bXK559/vmrChAm13TW2tR+V33TTTeta2hCn6fMBAwY0NL13dvOck08+eUJtbW2qqKiIK6644ukjjzxySdO5/fbbb/mxxx772p577jnpmWeeqb7oootGn3DCCQt23333lWV1lv7hD38YNX/+/KqHHnqoetddd11z7vrrrx8cEbHNNtusXLFiRZ9Zs2b1vf766wd//OMfX9RU58477xywbNmyPhERBxxwQI8GnIsXL66YNm1adVv1Jk6cuHrIkCHrx0KrPUjACQAAAFAw8+fPX5PpjB49usPh5MiRI9dcM3fu3MruDDh7wy233DLgscceGxAR8eEPf3heebjZZNSoUfUXXHDB84ceeuikhoaGOO+88zb9/e9//2LT+UMPPXRNKHnDDTcMLg84b7vttsERjbM8V6xY0ecvf/nLiFtuueUNAedNN900OCKioqIiDj744Joe+aKv9zVsypQpw9qq9/e///2pww47bIOfTWoNTgAAAICCWbp0aUVTedCgQR2eoTdw4MA11yxatKjw+dB11103pKl88sknt7ge6SGHHLJsq622WhkRcccddwwpPzd+/Pi6Lbfcstl1OO+9995BEY0zM/fbb7+lEW9eh7Np/c3tt99++fDhw9/wz2TWrFmP5pyn55ynd+b70TozOAEAAAAKZvDgwWvWfqypqelwQNn0KHVExLBhwwr/CPPjjz/ePyKiqqoq77nnnq0+sr/rrrsue/bZZ6tffPHFfitXrkzV1dVr1vncc889lz733HNvWIfz6aefrnr55Zf7pZTi0EMPXdr0Z1e+DmddXV1Mnz59na2/eeSRRy648sorn+/pfoqi8Ak9AAAAwMZm5MiRdU3lV199taqj18+fP3/NNZtuumlda3WLYNGiRZUREUOHDq2rqmr9j2OzzTarjYjIOce8efMqys81zc5sWocz4vX1N7feeuuVY8eOrdt2221Xjxs3bnXOec25u+++e0BNTU1FRM+vv8mbCTgBAAAACmby5Mkrmsr//Oc/B7RWd211dXUxc+bM/hERw4cPryv6+pvlUkpduv7d7373G9bhjHj9cfU999xzzbl3vOMdSyMibrnllsER63b9Td5MwAkAAABQMLvvvvuKoUOH1kdEPPDAA4MWLFhQ0dY1Ta6++uohK1eu7BMRsccee7wpjGsKCRsaWn9yffny5etNrjRs2LC6iMaZnLW1ree1c+bMqYpo/J6jRo2qLz/3lre8pW6LLbZYFfF6sHnPPfcMjnjjzMy11+G84447BkdETJo0afmIESPe0CY9b735IQIAAADQPn369ImjjjpqQUTEypUr+5x33nkj23vtz372s02byscee+ybNuQZOHBgfUTE4sWLW9275dlnn61u/4h71o477rgiIqK2tjbdc889rc5offjhhwdGRLzlLW9ZVb7+ZpOmNTTvu+++Qc8991zViy++uGb9zaY6TeWnnnqq/yuvvFL5wAMPDIp44yxP1h0BJwAAAEABffnLX57Tt2/fHBFxzjnnjH3sscf6tXXNRRddNPzWW28dGhGx7bbbrjj66KMXr11n/PjxqyIiHnvssQEtzeJ84IEHqp966qn+rfXVr1+/HBGxatWqrj033g6HHnrokqbyRRdd1GLYe+ONNw585plnqiMi9tlnnyXN1dl3333XrMN53nnnjYp4ff3Npjrbbbfd6rFjx67OOcf3v//9TZt2tbf+Zu8QcAIAAAAU0MSJE1d//etffzmi8XHxQw45ZLt77rmnxdDx17/+9fDPfvazW0Y07jZ+ySWXPNenz5ujoX322WdpRMS8efOqLrrook3WPv/aa6/1OeGEEya0Nb6RI0fWRkQ899xzbQavXXXAAQcs33HHHZdHRFx++eUjr7766sFr11mwYEHFqaeeukVE4wzY0047bW5zbZWvw/mb3/xm04jmZ2Y2rcPZVKdPnz5xyCGHNLv+5rhx43ZKKU1OKU3u+LejLa1ONQYAAABg/fWVr3xl7jPPPNPv4osv3vSVV17pu88++2x/+OGHL3z/+9+/eKuttlq1evXq9Pjjj1f/6U9/GtG0lmTfvn3zhRde+Nxee+21ork2P/WpTy380Y9+NLampqbiC1/4woSnn36633vf+94lKaV8//33D7zgggs2mzNnTt/tt99++ZNPPtni4+CTJ0+umTVr1iY33XTTsB/+8Icj999//5r+/fvniIjhw4fXjxs3rlt3b7/oooue33///bevra1N//Zv/7btcccdN/cDH/jAosGDBzdMmzZtwE9+8pPRL7/8cr+IiJNOOunVPfbYY2Vz7UyYMKH2LW95y6oXX3yxX2s7o++3335Lr7rqqhFNdSZOnLh85MiR62T9zcWLF1dMmzatzSUCKioqYrfddmv2e25IBJwAAAAABfbb3/72pUmTJq387//+73GLFy+u+Nvf/jbib3/724jm6m611VYrf/7zn7/w7ne/u8WdvseOHVt33nnnvXDSSSdttWrVqnTOOeeMPeecc8Y2na+urm74xS9+8dzUqVOHthZwnnXWWa/+4x//GL569er05S9/eYvyc0ceeeSCK6+88vlOfN0W7bXXXisuv/zyp4899titampqKn71q19t9qtf/Wqztet9/OMfn/ezn/1sVmtt7bnnnktffPHFfhGNmxGVr7/ZZO3P9txzz3W2e/pNN900bMqUKcPaqjdo0KD6pUuXPtzzI+pdHlEHAAAAKLgvf/nL855++ulHv/Od77y4zz77LBk9evTqfv365QEDBjSMHz9+1WGHHbbwl7/85bMzZ858vLVws8knP/nJ12644YYnDz744EXDhw+vq6qqyqNHj1595JFHLrjjjjuePP74419rq4299tprxc033/zkYYcdtnDMmDGrm9YL7UlHHnnkkpkzZz722c9+9tVJkyatGDRoUH3fvn3zmDFjVr///e9f+I9//GPmpZde+mJFReubzjetwxnx5vU3m0ycOHH12LFjVzcdW3+z96Sce/y3RTullDaPiJciIm46KmL0wHXX93vHTF13nQEAANCiwdufvc76ql1YGzNPn9l0OD7n/HJX2ps+ffr/VVZW7lBdXT1s4sSJT7dW94irj9j+tZWvVXWlvyIYXj289m8f+NuTvT0OWF/MnDlzm5UrVy6qq6t7YvLkye/tjjY9og4AAACsc0I/oLt4RB0AAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFgCTgAAAAC61VFHHTUhpTR53LhxOzV3fty4cTullCYfddRRE9bx0NgAVfb2AAAAAICNzzOHHb59/WuvVfX2OHpaxfDhtVtP/fuTvT0O2JAJOAEAAIB1rv6116rqFyzY4ANOoOd5RB0AAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAABXPAAQdsk1KavMsuu0xq7vzUqVMHp5Qmp5QmDx069G319fVvqvPiiy9WNtX5wQ9+MKrp8/r6+rjmmmsGn3TSSZvvtttuk4YPH75LZWXlboMHD37bpEmTdjjppJM2/9e//tW3B78edIiAEwAAAKBg9t5776UREY8//viAxYsXvynfueWWWwY1lZcsWVJx33339V+7znXXXTe4qXzwwQcvbSqfeeaZYz/wgQ9s96tf/Wqzhx56aOCiRYsq6+vrU01NTcXMmTP7/+pXv9ps55133vHSSy8d1u1fDDpBwAkAAABQMAcddNDSiIj6+vp0ww03DFr7/J133jm4/PjGG28cvHadW2+9dXBExIgRI+p23XXXlU2f19XVxahRo2o/9rGPzfv5z3/+3PXXXz/jjjvuePL3v//9M6eccsqrAwYMaFi5cmWfk046aasHH3ywuvu/HXSMgBMAAACgYN75zncuHzhwYENExM033/yG8HLFihXpkUceGRQRccABByyOiLj99tvfFHDee++9gyMi3v72ty8t//wzn/nM/JdeeunR3//+9y+eeuqpCw8++OBle++99/KPfexji37xi1/M+uc///nYpptuWrtq1ar0rW99a0xPfUdoLwEnAAAAQMFUVlbG5MmTl0ZE3HXXXW8IL2+99daBq1atSoMGDar/4he/OCciYtq0aYPL1+GcNWtW5bPPPlsdEbHPPvu8IeCcOHHi6n79+uWW+t56661rP/vZz74aEXHTTTcNbWho6LbvBZ0h4AQAAAAooL333rsm4s3rcDbN6Nx9991rDjrooJrq6uqGtdfhbGn9zeYsXLiwz4wZM/o+8MAD1dOmTaueNm1a9YABAxoiImpqaipmzJhhwyF6lYATAAAAoIAOPPDAZtfhbFp/c5999lnav3//vMsuuyyLeOM6nLfccsvgiIjhw4fXTZ48eWWs5amnnup77LHHjh83btxOI0aM2HX77bffaY899thxypQpO06ZMmXHL33pS1s01Z0zZ05lz31LaFshAs6U0u4ppa+nlK5PKb2cUlqVUqpJKT2VUro4pbR3N/XzjZRSbudr/+7oEwAAAKAz9tlnn2VNMymbZm2uXLkyPfzwwwMjXt+IqGnH9fJ1OFtafzMi4k9/+tOQXXbZZcdLL71009mzZ7c5O3P58uWFyJfYcK33P8CU0u0RMS0ivhkRB0fEuIjoGxEDI2LbiDguIu5IKf0upWRKNAAAALBRqKqqit12260m4vV1OG+//fYBK1eu7DNo0KD6vfbaa3lExAEHHLA04vV1OF955ZXKZ555ptn1N1955ZXKE088cauVK1f2GTBgQMPpp58++8Ybb5wxa9asR1asWPFgznl6znn61Vdf/VTTNTnntK6+MzSnCFOIx5beZ0fEnyPijoh4MSIqImLPiDgjGkPPT0REVUR8tJv63amN8891Uz8AAAAAnbL33nsvvfPOO4c0rcN50003rVl/s7KyMfZ517vetaxfv365aR3Op556ql/OjXsIrb3+5u9///vhS5curYiIuOyyy54+4ogjml2fc8GCBUXIlNhIrPczOCNiRkR8OCLeknP+Qs75ypzztJzzvTnnH0fE2yKi6W8NPpJS2rc7Os05P9bGa1l39AMAAADQWWuvw1m+/mZTndI6nDURjetw3nrrrYMjIoYNG/am9Tcff/zx6oiIoUOH1rcUbkZETJs2bUD3fxvonPU+4Mw5H5Zz/lPOub6F8/OjcRZnkw+um5EBAAAA9K599tlnef/+/RsiIq677rohDz300KCI19ffbFK+Dufdd989OCJiypQpNX36vDEaqqurSxERq1atSvX1zUYxsXTp0j5/+ctfRnT3d4HOWu8Dzna6pay8da+NAgAAAGAd6tevX951111rIiIuv/zykStWrHjD+ptNDjjggJqIiHvvvXfI008/3T/izetvRkRsu+22qyIiVq5c2ec3v/nN8LXP19XVxUc/+tEt5s2bV9UT3wc6Y0MJOPuVlZv/6wUAAACADVDT7MyampqKiDeuv9nkXe96V03fvn3zsmXL+rS0/mZExMc//vGFffv2zRERp5122pannnrquKuvvnrw7bffPuCnP/3piLe97W3bT506dZOmzY1gfbChBJz7lZWf7I4GU0rXp5TmppRWl95vTSmdnVJ6099edKDNzVt7RcTo7hg7AAAAsPFoWoezSXMzMwcMGJB32WWXNfuJDB06tH6PPfZYsXa9rbfeuvYHP/jBC3369IlVq1alX/ziF6OPOOKI7fbbb7/tTzvttAmPP/74gPe9732vff3rX5/dM98GOq7wAWdKqU9EnF320Z+6qemDI2JUNO7MPioaQ9TvRsSzKaUPdLLNl9p4TevimAEAAICNzL777ru8urq6oel47fU3mzTN9IyI2GOPPZauvf5mk89//vML/vGPf8w46KCDFg0fPryusrIyjxo1qnafffZZ8qtf/erZqVOnPltRUdHt3wM6q7LtKuu9L0bElFL5qpzz9C6292hE/C0i7o+I2dEYcE6MiGMi4pCIGBYRV6aUDs85X9vFvgAAAGCjVDF8eG1vj2FdWBffs7q6Oq9YseKhtuqde+65s88999x2zbw8+OCDlx188MHPtHT+sMMOW9paBnPllVc+HxHPt3R+1qxZj7ZnHNAehQ44U0r7RcT3SodzI+LTXWzyJznnbzTz+X0RcWlK6eSIuDAiKiLi1ymlrXPOKzvQ/vg2zo8OszgBAADYCGw99e/dssQcQGEDzpTSjhHx12j8Disj4t9yznO70mbOeVEb53+ZUtojIk6IiLERcVREXNaB9l9u7XxKqb1NAQAAAABR0DU4U0pbRsT1ETE8GndNPzrnfPs66v6XZeX9WqwFAAAAAPS4wgWcKaWxEXFjNM6gzBHxyZzz1etwCE+Ulcetw34BAAAAgLUUKuBMKY2MiBsiYqvSR5/LOV+6joeR13F/AAAAAEALChNwppSGRsR1EbFD6aOzc84/74Wh7FBWbtfOYwAAAABAzyhEwJlSGhAR/xsRu5U++k7O+fu9NJyTy8q39dIYAAAAAIAoQMCZUuobjbulv7P00Xk55692op3jUkq59PpGM+d3Silt00YbJ0XEiaXDV0vjAgAAAAB6SWVvD6Ad/iciDimVb46I36SU3tpK/dU556c60c/kiPh1SumWiLg2Ih6NiAXR+Gc0KSKOKRtHfUSclHNe1ol+AAAAAIBuUoSA88iy8rsi4p9t1H8hIiZ0sq+KiDio9GrJgog4Ief89072AQAAAAB0kyIEnOvK/0XECRGxZ0TsGhGbRcSIiEgRsTAiHomIf0TEJTnnJb01SAAAAADgdet9wJlzTt3UziURcUkr5+dGxG9LLwAAAACgANb7TYYAAAAAAFoi4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAoCSlNDmlNPn0008fu/a5qVOnDm46P3Xq1MG9MT7erLK3BwAAAABsfP7nm/dtv6JmdVVvj6On9R/Ut/Yj//n2J3ui7alTpw4+/PDDt2tv/fPOO+/50047bUFPjGV98swzz1T9/Oc/H3XbbbcNfvbZZ6uXLl1aUVlZGUOHDq0bP378qp122mn5vvvuW/P+979/yYgRI+p7e7y97aijjppw1VVXjWg6vuKKK/71oQ99aElr16SUJkdEHHnkkQuuvPLK53t4iG0ScAIAAADr3Iqa1VUrltZu8AEn69Y555wz8qtf/er4lStXvuGp5fr6+pg7d27V3Llzq6ZPnz7okksu2fR973vfa1OnTn22t8a6vvqv//qvcW0FnOsbAScAAABAwR1zzDHzPv/5z89trc6WW25Zu67G0xt++ctfbvKlL31pi4iIfv365Q9+8IPz3/3udy/ZYostVuec46WXXqqaNm3awOuvv37ok08+OaAzfRx22GFLc87Tu3fk65fHHntswGWXXTb0mGOOWdzbY2kvAScAAABAwW266aZ1e+yxx8reHkdvqauri69+9aubR0QMHDiw4aabbprx9re/fcXa9Uqh3ewHH3yw+sEHH+y/zge6nhs2bFjdokWLKr/97W+P+8hHPrK4T59ibN9TjFECAAAAQAtuueWWgfPnz6+KaJzN2ly4WW633XZbeeKJJ762bkZXHJ/97GdfjYiYMWNG/0svvXRYLw+n3QScAAAAABuR008/fWzTTuAREQsWLKg488wzx2y//fY7DB48+G0ppcnnn3/+iLWv+/3vfz/sPe95z1ZjxozZqV+/frsNHjz4bW9961u3P+OMM8bMmzevoj19d0cbzXnuuef6NpW32WabHpvJ2tYu6kcdddSElNLkcePG7VQaV9WJJ564+YQJE97av3//XYcPH77L/vvvv81f/vKXIT01xq748pe/PG/EiBF1ERHf+c53xtXXF2MPJgEnAAAAwEbq0Ucf7bfzzjvv8KMf/WjsjBkz+tfU1LwpZJw3b17Fnnvuud0nPvGJrf/xj38Mf/XVV/uuXr061dTUVDz++OMDzj333LGTJk1660033TSwpX66o43W9OvXLzeVn3zyyfXi0fPbb799wO67777Db37zm81eeOGFfitXruyzaNGiyttuu23ov/3bv237qU99avOutH/++eePaApbTz/99LHdMebBgwc3nHbaaa9ERDz99NPVv/71rzfpjnZ7moATAAAAYCP1wQ9+cOu5c+dWHXvssXP/+te/PnX77bc/+ctf/vLZHXbYYWVExIoVK9L++++/3b333ju4oqIijjjiiAW//OUvn73ppptm/OMf/5h51llnzRo2bFjdwoULK4888shtn3rqqb5r99EdbbTl7W9/+/Km8h//+MdR11xzzZtmV65LK1as6PORj3xk65qamopTTz311WuvvXbmzTff/OS3v/3tl0aNGlUbEfHrX/96s//6r//atDfH2ZwvfelL85rG+L3vfW9MXV1dbw+pTTYZAgAAACi4uXPnVk6bNq26pfNjx46tGzdu3JuSqqeffrr/n//8538deeSRS5o+22effdaEhV/+8pfHPvHEEwMGDx5c/7//+79PlZ+LiDj00ENrPvnJTy7ce++9J82bN6/qS1/60rhrrrnmufI63dFGWyZNmrT6gAMOWHzLLbcMXbVqVfrABz6w3Vvf+tblBx100OI999yzZp999lk+ZsyYdZbUvfbaa5VLly7N11xzzVPvec97apo+P+CAA5Yfc8wxr73jHe/Yfs6cOVXf/e53x51wwgkLx44du96kiAMGDMinn376K//+7//+lueff776wgsvHPHZz352QW+PqzVmcAIAAAAU3GWXXTZqypQpO7b0Ouecc5qdKXjUUUfNLw83yy1evLjPJZdcMioi4uyzz561djDZZLvttlt9xhlnvBIRce211w5fsmRJn+5so70uu+yy59/61reuaf+xxx4b8JOf/GTMhz/84W3Hjh27y4QJE9567LHHjr/zzjsHdLTtzvjoRz86rzzcbDJhwoTab3/72y9FNM70vPDCC9+03mlv+8IXvjB/zJgxqyMifvCDH4ypra3t7SG1SsAJAAAAsJH62Mc+trClc9dee+3gpjU5P/axj7W64/iBBx64NCKirq4u3XXXXWsCxO5oo73GjBlT98ADD8z40Y9+9MIOO+zwpiD1hRde6HfppZduus8++2x/xBFHbNmZELUjPvWpT7U46/HjH//4osGDB9dHRNxyyy2d2nDotNNOW5Bznp5znn7uuefO7uw4m1NdXZ2bAueXXnqp309/+tOR3dl+d/OIOgAAAEDBffGLX3ylMyHXHnvssaKlc9OmTVsTMm6xxRa7tLfNWbNmVXVnGx3Rr1+/fMYZZ8w/44wz5j///PNVN9xww6AHHnhg4PTp0wc+8sgjA+vq6lJExNVXX73JoYceWnXHHXc8VVnZ/fFYVVVVfsc73tHsbNWmce6www7L77vvvsEzZ85cLzZFWttpp502/yc/+cnol19+ud8555wz5tRTT11QXV2d275y3TODEwAAAGAjNWrUqBbXfpw7d26nQsbly5evyZu6o43OmjBhQu2nPvWp1375y1++/MADD8x8+eWXH/nMZz7zap8+jU3fe++9gy+66KIe2SV86NChdW0Fp00b+SxZsmS9nIBYVVUVX/7yl1+JiJg9e3bf8847b72dxble/gECAAAA0PNaC+Hq6+vXlO+8884n+vbt267Ze1tuueWaBRu7o43ustlmm9X/7Gc/m5VzjgsuuGB0RMSVV145/NRTT23xMf3OSil1d5O94tOf/vSCc845Z8wLL7zQ79xzzx1z2mmnze/fv/96N4tTwAkAAADAm4wYMWJNOjl69Oi6rbfeusOhY3e00d0+85nPzGsKOF944YUWd57vikWLFlXW1dW1GiDPmzevKiJiyJAh680O6murrKyMs88+e/anP/3pLefOnVv1ox/9aNTXvva1ub09rrV5RB0AAACAN9ltt93WrCF58803D+qtNrrbhAkT1oSsPTXTsra2Nt17770tbpRUW1sbTz755ICIiO22267FdVDXB5/61KcWbr311isjIs4///zRNTU16930VAEnAAAAAG9y+OGHL6murm6IiPjFL36xWUNDQ6+00R4dafeOO+4Y2FQeP378qh4ZUET85je/GdHSud///vfDlyxZUhERccABByzpqTF0h4qKivjKV74yOyJi/vz5VT/4wQ827e0xrU3ACQAAAMCbjBw5sv64446bGxHx0EMPDTzxxBPHl6+pubaXXnqp8txzz33DRjTd0UZ7/PnPfx763ve+d6u77rqr1R3J58yZU/HFL35xfNPx4YcfvqijfbXXH/7wh1HXXXfdm2atvvjii5Vf/epXN4+IqK6ubjjllFMWdKb9888/f0RKaXJKafLpp58+tqvjbc2xxx772sSJE1dERPz0pz8d3ZN9dYY1OAEAAABo1rnnnjv77rvvHvzPf/5z4MUXX7zp3XffPfjYY4+dN3ny5OWDBg1qWLBgQeWjjz5affPNNw+5/fbbh2633XYrTj/99Pnd3UZbGhoa4tprrx1+7bXXDp84ceKKgw8+ePGUKVOWjRs3rrZfv34Nr776atXtt98+6A9/+MOohQsXVkZE7Ljjjss/+9nPdqif9ho+fHhd//79Gz7wgQ9sd+KJJ845/PDDF1dXVzfcfffdA3/84x+PaVp/86yzzpo9bty49XYNziZ9+vSJr371q7M//vGPb71o0aL1Lk9c7wYEAAAAwPqhf//++dZbb33q6KOP3vL6668fNnPmzP7/8R//8ZaW6g8aNOhN0zO7o422jBgxoq5///4NK1as6DNz5sz+M2fObHUm51577bXkyiuvfLaqqqqjXbVL//79G/7nf/7nmQ984APb/vznPx/985///E2zHo877ri53/jGN+b0yAB6wMc+9rFF3/ve95Y//vjjLa4t2lsEnAAAAAC0aPjw4Q3XXXfdM9ddd92giy++eMR99903aN68eX1XrlyZBg0a1DB+/PhVu+6667LDDjts0ZFHHtnsepLd0UZrDjnkkGVz5859+Jprrhlyyy23DH7kkUcGPv/88/0WL15cWV9fH4MGDWoYN27cql122WXZRz/60dcOO+ywpV3/k2ndvvvuu3zatGlPfOc73xl90003DZ07d27f/v371++0007LP/vZz8750Ic+tF6vvdmcr3/967M+/OEPb9vb41hbyjn39hgoSSltHhEvRUTcdFTE6IFtXNCN3jtm6rrrDAAAgBYN3v7sddZX7cLamHn6zKbD8Tnnl7vS3vTp0/+vsrJyh+rq6mETJ058urW6//PN+7ZfUbO6Z6bPrUf6D+pb+5H/fPuTvT0O1o2jjjpqwlVXXTVi7Nixq2fNmvVob49nfTRz5sxtVq5cuaiuru6JyZMnv7c72jSDEwAAAFjnhH5Ad7GLOgAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwrLJEAAAAAB0gyuvvPL5iHi+l4ex0TGDEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAgO6SIyJyzqm3BwKsn3ri/iDgBAAAALrLkpxzXX19fWVdXZ3MAXiDhoaG1NDQUBER9RFR013tutkAAAAA3eWh+vr6pRGRFy5cOLy3BwOsXxYvXjw455zq6+uXRcSj3dWugBMAAADoLldFRH19ff2iefPmbfrKK69sunz58n45594eF9CLGhoa0muvvTZkzpw5oxsaGpblnFdHxI3d1X5ldzUEAAAAbNwmT5789PTp06fW1tYelnOuW7hw4YiFCxeO6NOnT0NKqSGlJOmEjUzOOTU0NFSU3petXr365Yh4OCKmdVcfAk4AAACgO30zIl6uq6s7pa6ubkGfPn369+nTpzo8RQobs/r6+vqmmZsPR8RpkydPruuuxpNp4uuPlNLmEfFSRMTNW20do6uqenlEAAAAbMhera2Ndz37TNPh+Jzzy93V9vTp08dHxJ6l1+YRMbC72gYKpyYa19y8MSKmdWe4GWEGJwAAANADJk+e/FI0TuL5U2+PBdiwmR4OAAAAABSWgBMAAAAAKCwBJwAAAABQWIULOFNKW6SUzkkpzUgpLUspLUwpTUspnZlSGtCN/bwnpfTXlNLLKaVVpfe/ppTe0119AAAAAABdU6hNhlJKh0fEHyJiSNnHAyJi99LrxJTS+3LOT3ehjz4RcVFEnLDWqXGl1xEppV9HxMk554bO9gMAAAAAdF1hZnCmlHaNiCuiMdysiYivRMReEXFgRPyqVG27iPjflNLgLnT1nXg93HwoIj4SEVNK7w+VPj8xIr7dhT4AAAAAgG5QpBmc50VE/4ioi4hDcs73lJ27OaX0r4j4QTSGnGdExDc62kFKabuI+FLp8IGI2DfnvKJ0PC2ldE1E3BaNs0XPTCn9tiuzRQEAAACArinEDM6U0pSI2Kd0+Ju1ws0m50TEk6Xy51NKVZ3o6gvxeuj7ubJwMyIics7LI+JzpcPKiPhiJ/oAAAAAALpJIQLOiDiirHxxcxVK62FeWjocFhEHdKSDlFKKiA+UDmfknO9toZ97I2Jm6fADpesAAAAAgF5QlIBz79L7soiY3kq928rK7+xgH1tGxNhm2mmtn3ERMaGD/QAAAAAA3aQoa3BuX3p/Oudc10q9Gc1c0147tNBOe/p5rj0dpJQ2b6PK6Pa0AwAAAAA0Wu8DzpRSdUSMLB2+3FrdnPNrKaVlETEwIsZ3sKvy8LHVfiLipbJyR/p5qe0qAAAAAEB7FeER9cFl5Zp21F9Weh/Ug/0sKyt3tB8AAAAAoJus9zM4I6K6rLy6HfVXld7792A/q8rKHemnrdmeoyNiWgfaAwAAAICNWhECzpVl5b7tqN+v9L6iB/vpV1Zudz8551YffbchOwAAAAB0TBEeUV9aVm7P4+ADS+/teZy9s/0MLCt3tB8AAAAAoJus9wFnznllRCwoHba6C3lKaXi8Hj52dEOf8tmVbe12Xv6ouY2DAAAAAKCXrPcBZ8kTpfdtUkqtPVY/qaz8ZCf7WLud7u4HAAAAAOgmRQk47yy9D4yIya3U26+sfFcH+3guImY3005z9i29z4qI5zvYDwAAAADQTYoScP6trHx8cxVSSn0i4hOlw0URcUtHOsg554i4unQ4KaX0jhb6eUe8PoPz6tJ1AAAAAEAvKETAmXO+PyLuKB2ekFLas5lqZ0TE9qXyeTnn2vKTKaX9U0q59Lqkha5+EhH1pfJPU0r912qjf0T8tHRYV6oPAAAAAPSSQgScJZ+PiBURURkR16eU/j2l9I6U0gEppV9GxA9K9Z6KiHM600HO+amI+GHpcPeIuCul9OGU0u4ppQ9H42Pvu5fO/zDn/K/OfhkAAAAAoOta27BnvZJzfqgUMv4hIoZExH83U+2piHhfznlpF7r6SkRsGhGfjIhdI+LyZur8JiK+2oU+AAAAAIBuUKQZnJFz/ntE7BwRP47GMHN5NK63+UBEnBURu+acn+5iHw055xMi4n3RuCbn7IhYXXq/OiLem3M+Mefc0JV+AAAAAICuS/bIWX+klDaPiJciIm7eausYXVXVyyMCAABgQ/ZqbW2869lnmg7H55xf7s3xAHRGoWZwAgAAAACUE3ACAAAAAIVVmE2GNjbb3nZrbL755r09DAAAADZgg19+OWL8+N4eBkCXmMEJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYAk4AAAAAoLAEnAAAAABAYQk4AQAAAIDCEnACAAAAAIUl4AQAAAAACkvACQAAAAAUloATAAAAACgsAScAAAAAUFgCTgAAAACgsAScAAAAAEBhCTgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwKnt7ALxBRVPhlVde6c1xAAAAsBFY6/89K1qqB7A+Sznn3h4DJSml3SNiWm+PAwAAgI3SHjnnB3p7EAAd5RF1AAAAAKCwzOBcj6SU+kXETqXDeRFR34HLR8frsz/3iIhXu3FobBz8hugOfkd0ld8QXeU3RHfwO6KrivQbqoiIUaXyoznnVb05GIDOsAbneqT0L5JOPQ6QUio/fDXn/HK3DIqNht8Q3cHviK7yG6Kr/IboDn5HdFUBf0Mv9PYAALrCI+oAAAAAQGEJOAEAAACAwhJwAgAAAACFJeAEAAAAAApLwAkAAAAAFJaAEwAAAAAoLAEnAAAAAFBYKefc22MAAAAAAOgUMzgBAAAAgMIScAIAAAAAhSXgBAAAAAAKS8AJAAAAABSWgBMAAAAAKCwBJwAAAABQWAJOAAAAAKCwBJwAAAAAQGEJOAEAAACAwhJwbgBSSluklM5JKc1IKS1LKS1MKU1LKZ2ZUhrQ2+Nj/ZVSyu183drbY2XdSyltmlI6LKX0rZTStSml+WW/iUs60d57Ukp/TSm9nFJaVXr/a0rpPT0wfNYD3fEbSikd14F71XE9+43oDSml3VNKX08pXV92/6hJKT2VUro4pbR3B9tzL9rIdMdvyL1o45VSGpJSOrr0/1u3pZSeTiktTimtTinNTSndmlL6ckppRDvb2yul9IeU0gsppZUppVdTStellD7S098FYENW2dsDoGtSSodHxB8iYkjZxwMiYvfS68SU0vtyzk/3xviAQpvTHY2klPpExEURccJap8aVXkeklH4dESfnnBu6o0/WG93yG2LjlVK6PSL2aeZU34jYtvQ6LqV0aUR8Kue8upW23Is2Qt35G2KjNSUi/qeFc6MiYr/S68yU0sdyzte11FBK6RsR8bV440SjzSLikIg4JKV0TER8MOe8sjsGDrAxEXAWWEpp14i4IiL6R0RNRHw3Im4pHR8dEZ+KiO0i4n9TSrvnnJf21lhZ7/0iIi5o5fyydTUQ1lsvRsSMaPwP8I76TrweKDwUET+IiGciYuuI+HJE7BoRJ0bEvIj4jy6PlPVVV35DTQ6NiNmtnH+5C22zfhpbep8dEX+OiDui8bdUERF7RsQZ0RhOfiIiqiLio6205V60cerO31AT96KNz0vR+P9Z00vlV6IxpNw8Ij4YEUdGxMiIuCalNCXn/MjaDaSUTo6I/ywdPhMR/x0Rj0bjb/TzEXFARLwvIn4b7fsdAlAm5Zx7ewx0UtnfSNdFxL4553vWOn9mNP7He0TEN3PO31i3I2R9l1JqugH4ffAmKaVvRsS0iJiWc56TUpoQEc+VTv8u53xcO9rYLiIej8a/UHsgGu9VK8rOD4iI26JxxnldRGxvxvmGo5t+Q8dFxMWlwy1zzs93/0hZX6WUpkbEpRFxZc65vpnzIyPirmj8C92IiP1yzrc3U8+9aCPVjb+h48K9aKOUUqpo7rezVp0jIuKvpcO/5pyPXOv8JhHxbEQMjcaAfXLOeX55H6XrDy99dEDO+dZu+QIAGwlrcBZUSmlKvP64zW/WDjdLzomIJ0vlz6eUqtbJ4IANQs75P3POU3POXXnM+Avx+tMCnysPFEp9LI+Iz5UOKyPii13oi/VMN/2G2IjlnA/LOf+ppXChFBCcUfbRB1to6gvhXrRR6sbfEBuptsLNUp2/RcTM0mFzSyKcGI3hZkTEWeXhZlkfp0ZEU19ndmqwABsxAWdxHVFWvri5CqX1oy4tHQ6LxsceANaJlFKKiA+UDmfknO9trl7p86b/KfhA6TqA9rqlrLz12ifdi2iHVn9D0E5Ny4FVN3PuiNL7koi4qrmLc84vR8SNpcMDU0qDu3V0ABs4AWdxNe32uCwa14JpyW1l5Xf23HAA3mTLeH3ts9taq1h2flxETOipAQEbpH5l5eZmWrkX0Za2fkPQqpTSxIh4W+lwxlrn+kbjRkUREfe0sZFV0z2oXzQumQFAOwk4i2v70vvTOee6VuqV/wt2+xZrsbH7t5TSEyml5SmlpSmlf6WUfpdSMuuXrtihrDyjxVpvPu9eRUsuTinNTimtTinNTyndm1L6dkppXG8PjF61X1n5yWbOuxfRlrZ+Q2tzLyJSSgNSStumlE6PxmCyaRmMn6xVdbto3NQqwj0IoMcIOAsopVQdjbv0RbSxS2PO+bV4fQfs8T05Lgpth2j8j6j+ETEoIraJxt1Eb04p/TWlNLS1i6EFm5eV29pR9qWysnsVLdk/IsZE407HIyLi7RHxlYh4urQ7LRuZlFKfiDi77KM/NVPNvYgWtfM3tLb9w71oo5RSOi6llEsbdS6LiKeicd+DzUpVvhcRf1zrMvcggHWgsu0qrIfK12OpaUf9ZRExMBqDKyi3PCKuiYibovFvjGsiYlQ0zmQ4JRr/o/2IiLg6pXRwzrm2l8ZJMXXkXrWsrOxexdqejcY1y+6J1//nb6uIOCoaNwSpjogLU0o553xR7wyRXvLFeP3Rz6tyzs0t2+NeRGva8xtq4l5ESx6OiJNyztOaOeceBLAOCDiLqXzh6tbWcGmyqvTevwfGQrGNyzkvaubzG1JKP42IayNi12gMPD8dEeevw7FRfB25V60qK7tXUe6vEfG7nHNe6/NpEXFFSumwaAwcqiLixymla3LOr67rQbLupZT2i8bZUhERc6Px31PNcS+iWR34DUW4F9HobxHxQKncPxo3pfpQRPy/iPiflNIXcs5T17rGPQhgHfCIejGtLCv3bUf9poXTV/TAWCiwFsLNpnNzonE2QtOszc+tizGxQenIvap8gwf3KtbIOS9uJlAoPz81Ir5VOhwQESesk4HRq1JKO0Zj4FQZjfeaf8s5z22hunsRb9LB35B7ERHR+N/OOefHSq9pOefLc85HRuPSTltF41NPx611mXsQwDog4CympWXl9jy6MLD03p7H2WGNnPOzEXFD6XCblNLY1urDWjpyrxpYVnavoqMuioim4GG/1ipSfCmlLSPi+ogYHo07Xh+dc769lUvci3iDTvyG2su9aCOVc/59RPw5Gv//+mcppU3KTrsHAawDAs4CyjmvjIgFpcPNW6ubUhoer/+L8qXW6kILnigr2x2UjihfSL/Ve1W8cSF99yo6pDTrqunfi+5TG7DSX7TdGBFjozFI+mTO+eo2LnMvYo1O/obaxb1oo9f0OxoYEe8u+9w9CGAdEHAWV1PotE1KqbW1VCeVlZ/swfGw4WrxcSxoQ3k4PqnFWm8+715FZ7hXbeBSSiOj8amCrUoffS7nfGk7LnUvIiK69BvqCPeijde8svIWZeWnonGmcIR7EECPEXAW152l94ERMbmVeuWPx9zVc8NhA7ZDWXl2r42CInouXv/NtPWo3r6l91kR8XxPDYgNU0ppVESMLB26T22AUkpDI+K6eP3fSWfnnH/ezsvdi+jqb6i9fbgXbdzKZ+2uebw857w6Iu4vHe6ZUmptHc6me9SqeH0zIwDaQcBZXH8rKx/fXIWUUp9oXPA6ImJRRNzSs0NiQ1Nao+rg0uEzOedZvTkeiqW0GUPT41qTUkrvaK5e6fOmGQtXt7aJA7TgpIhIpfJtvTkQul9KaUBE/G9E7Fb66Ds55++393r3Irr6G+oA96KN27+VlR9d69zfSu9DIuLI5i5OKW0e/7+9+wmxqorjAP49WJCLVgpSBEoF0aJFYAuzyMCVZptaRC0Ki1y4KGlRRAUtyqKNYkR/DA0CoZAWRS3KNmlCLapFG1EqaiVBRqJl0a/FvTaDjI6jz5m58z4fGJhz332H8+Bw3rnfe859ydq+uK+q/pjqPACmJuAcqKr6KskXffHh1tqqKU57IsmN/f/bq+rvKc5hTLXWNpzr8QattWVJ9mbi1x5fm5WGsdBsy8S2rB2ttcWTX+zLO/riP/35kCRpra1ord08zTl3JXmuL55MsuuSN4xZ0690+iDJ6v7Q9qp65gKq2hZj0VgaRR8yFo231tpDrbUrpjlnS5J1ffGHTFynnbYzye/9/y+11pac8f5F6ebai/pDr1xUowHGUHNzerj6idaBJIvTbYN4Md0qzcVJ7kt3Fznpnvuy0l1AJmut/Zjk8nQh5sF0W/FOpttatSbJpkxss9qfZG1V/TXb7WTutNZuS3L9pENLMzHhPpBusv6/qtp9lnq2JnmqL36T5OUkR5Jcl+TJJKcvGrdW1dOjaDvzw8X2odbamnTfaweTfJjkuyRH+5evTXJv/3d6xdTmqnIzZgFpre3NxGqnz5M8nnM/4/BUVR06S13GojE0ij5kLBpv/Zz5ynRz5v3pxo3j/bGbkjyQiQD9VJL1VfXZFPVsSvJ6XzyS5IV0Kz2vTtcv7+xf21NV91+CjwKwoAk4B661tiHJu+m2O0zlULov2cOz1yqGoJ+sLZ/uvHSTuUeq6tglbRDzTmttd5IHz/f8qmpTHe8fl/FWko3nePvbSR6tqn9n0kbmt4vtQ5NChemcSLKlqt6cQfMYgNbaTCeqP1XVirPUZSwaQ6PoQ8ai8TaDOfMvSTZW1afnqOv5JM9mIgw/08dJ7qmqP2faToBxJ+BcAFpry5M8lmR9kmvS3Tk8nOT9JK9W1Yk5bB7zVGvtjnQPMl+VbvXB0nRB+fEkPyf5Msk7VXVwzhrJnBpVwDmpvnXpVpbfkq6//Zrk6yRvVNUnF95S5qsRBJxXJrk73Ti1MslV6frOZUl+S/J9kn1JdlbV0bDgjDLgnFSnsWiMjCjgNBaNsdbaDemus1an25WwLMmSdDufjib5NslHSd47n+uu1tqtSTYnub2v61i6VcG7qmrP6D8BwHgQcAIAAAAAg+VHhgAAAACAwRJwAgAAAACDJeAEAAAAAAZLwAkAAAAADJaAEwAAAAAYLAEnAAAAADBYAk4AAAAAYLAEnAAAAADAYAk4AQAAAIDBEnACAAAAAIMl4AQAAAAABkvACQAAAAAMloATAAAAABgsAScAAAAAMFgCTgAAAABgsAScAAAAAMBgCTgBAAAAgMEScAIAAAAAgyXgBAAAAAAGS8AJAAAAAAyWgBMAAAAAGCwBJwAAAAAwWAJOAAAAAGCwBJwAAAAAwGAJOAEAAACAwfoP5zIh3PcBOzMAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 1200x800 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "bh = LatticeBoltzmannBoundaryHandling(method, dh, 'src', name=\"bh\")\n",
+    "\n",
+    "inflow = UBB(velocity_info_callback, dim=dim)\n",
+    "outflow = ExtrapolationOutflow(stencil[4], method)\n",
+    "wall = NoSlip(\"wall\")\n",
+    "freeslip = FreeSlip(stencil, (0, -1, 0))\n",
+    "\n",
+    "bh.set_boundary(inflow, slice_from_direction('W', dim))\n",
+    "bh.set_boundary(outflow, slice_from_direction('E', dim))\n",
+    "bh.set_boundary(wall, slice_from_direction('S', dim))\n",
+    "bh.set_boundary(wall, slice_from_direction('T', dim))\n",
+    "bh.set_boundary(wall, slice_from_direction('B', dim))\n",
+    "bh.set_boundary(freeslip, slice_from_direction('N', dim))\n",
+    "\n",
+    "plt.figure(dpi=200)\n",
+    "plt.boundary_handling(bh, make_slice[:, :, domain_size[2]//2])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def timeloop(timeSteps):\n",
+    "    for i in range(timeSteps):\n",
+    "        bh()\n",
+    "        dh.run_kernel(kernel)\n",
+    "        dh.swap(\"src\", \"dst\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "timeloop(500)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vel_profile = dh.gather_array('velField')[-2, :, domain_size[2] // 2, 0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "0.07442458175252796"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "np.max(vel_profile)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[<matplotlib.lines.Line2D at 0x7f6797ef7e80>]"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkLElEQVR4nO3dd3yV9d3/8deHBMKeCSAzQALIHgEsKq3iwNFSb62CC1etWm5tq7WOOop3h7Wtt62r3oIDqaComDpKXXUiEJAVZlgmYSWElYSQ9fn9kWN/MQ1wkCTXycn7+Xjw4FzX9TV5B5P3ufK9lrk7IiISvRoFHUBERGqXil5EJMqp6EVEopyKXkQkyqnoRUSiXGzQAaqKj4/3xMTEoGOIiNQrS5YsyXX3hOq2RVzRJyYmkpaWFnQMEZF6xcy2Hm6bpm5ERKKcil5EJMqp6EVEopyKXkQkyqnoRUSinIpeRCTKqehFRKJcxJ1HLyISaUrLysnac5AtuwvYuruQ3fmHauXz9O3civOHdKnxj6uiFxEBSr4q89yCfxf65twCtu4uIGvPQUrLv/7sDrOaz3D+kC4qehGR41FcWk7mnkK27i5gc27F31t2F7Ilt4DsvQcpq1TmLeNi6dmhOQO7tuG8ISeQ2KEFifEt6NmhOQkt47DaaPpaoqIXkahTXu6kb9vPoi15bM7NZ+vuQrbsLiB7z0Eq75i3ioslMb4FQ7q1YeKwLvTs0ILEDs1JjG9BhxZN6lWZH4mKXkSiwo59RXy8IYePN+TySUYueQXFALRqGkuv+BYM796OC4Z1rSjz+IpCbx9FZX4kKnoRqZcOFpexaEseH6/P4aMNOazfmQ9AfMs4vtM3gVP7xjO2TzwdW9WvaZbaoKIXkXrB3Vmz/cC/99oXbcmjuLScJrGNGJ3YngtHdOPU5AT6d25Fo0YNu9irUtGLSMTKOXCITzJy+Hh9Lh9tyCU3dFpj304tueKknozrm8DoxPY0axITcNLIpqIXkYhRVFLGkq17+GhDDh+tz2XN9v0AtG/RhFOS4jk1OZ5TkxPo3KZpwEnrFxW9iARqd/4h/r58Gx+sy2Hh5t0UlZTTOMYY2bMdPz+7H+OSExjYpbWmY46Dil5E6lxZufPR+hzmLM7k3TU7KS13eie0YNKoHpyaHM9JvTvQIk71VFPC+pc0swnAI0AM8LS7/67K9jjgeWAksBu4xN23mNllwM8rDR0CjHD3ZTWQXUTqma27C3g5LYu5S7LYsb+I9i2acNXYRC4e1Z2+nVoFHS9qHbXozSwGeAw4E8gCFptZqruvrjTsWmCPuyeZ2STgQSrKfhYwK/RxBgPzVPIiDUtRSRlvr9rOnMWZfL4pj0YG3+6bwH3fHcD4EzvRJFb3Vqxt4ezRjwYy3H0TgJnNBiYClYt+InB/6PVc4FEzM3evfHOIycDs404sIhHP3VmZvY+X0jJ5fdk2DhSV0qN9c247qy8XjuzGCW2aBR2xQQmn6LsCmZWWs4Axhxvj7qVmtg/oAORWGnMJFW8I/8HMrgeuB+jRo0dYwUUk8uwpKGbesmzmLM5k7Y4DxMU24tzBJ3BxSnfG9GqvA6oBqZOjHWY2Bih091XVbXf3p4CnAFJSUry6MSISmcrLnU8ycpmTlsk76TspLitnSLc2PPD9QXxvaBfaNGscdMQGL5yizwa6V1ruFlpX3ZgsM4sF2lBxUPYrk4AXjyOniESYzLxC5i6pOLCavfcgbZs35tIxPbg4pTsDurQOOp5UEk7RLwaSzawXFYU+Cbi0yphUYAqwALgIeP+r+XkzawRcDJxaU6FFJBhFJWX8c/VOXlqcyacbK2ZmT0mK585z+3PmgE7ExeoK1Uh01KIPzblPBeZTcXrlDHdPN7NpQJq7pwLTgZlmlgHkUfFm8JVxQOZXB3NFpP7ZX1TCzAVbeebTzeTmF9O1bTNuGZ/MRSO70a1d86DjyVHY10+MCV5KSoqnpaUFHUNEgNz8Q8z4ZDMzF2zlwKFSxvVN4Ien9uLkPvE6sBphzGyJu6dUt02XnonIf8jaU8hTH21izuJMisvKOXfQCdz4nT4M6tom6GjyDajoReTfNuw8wBMfbiR12TbM4ILhXfnRt/vQJ6Fl0NHkOKjoRYRlmXt5/IMM/rl6J80ax3DltxL54bheurApSqjoRRood+ezjbt5/F8ZfJqxm9ZNY7n59CSuOrkX7Vs0CTqe1CAVvUgDU17uvLNmJ49/kMHyrH10bBXHXef259IxPWmpO0ZGJf1fFWkgSsrKSV22jSc/3MiGXfn0aN+cX18wiAtHdKNpY53/Hs1U9CJRrqikjJfSMvnrh5vI3nuQ/p1b8cikYZw3+ARiY3TnyIZARS8Spape5DSyZzumTRzI6f07YqZz4BsSFb1IlDlQVMLTH29mxiebOXColG/3TeCm7/RhdK/2KvgGSkUvEiUOlZbxwudf8tgHGeQVFDNhYGemnp6ki5xERS9S35WVO699kc3D76wne+9BTk7qwO1n92do97ZBR5MIoaIXqafcnXfX7OKh+WtZvzOfwV3b8OCFQzglOT7oaBJhVPQi9dCizXk8+I+1LNm6h17xLXjs0hGcM6izbjQm1VLRi9Qja7bv56H563h/7S46tY7jNxcM5gcp3Wis0yTlCFT0IvVAZl4hf3pnPfOWZdMqLpZfTOjPVWMTadZEFzrJ0anoRSJYzoFDPPZBBrMWbqWRGT8a14cbv92HNs31HFYJn4peJAIdKCrh/z7ezNMfb+JQaTkXp3TnlvHJdG7TNOhoUg+p6EUiSNVz4c8bfAK3ntWX3rofvBwHFb1IBKh6LvwpSfHcPqEfQ7q1DTqaRIGwit7MJgCPUPFw8Kfd/XdVtscBzwMjgd3AJe6+JbRtCPBXoDVQDoxy96Ka+gJE6jOdCy914ahFb2YxwGPAmUAWsNjMUt19daVh1wJ73D3JzCYBDwKXmFks8AJwhbsvN7MOQEmNfxUi9VD6tn1M+/tqFm7Oo3d8Cx6/rOJceN2PRmpaOHv0o4EMd98EYGazgYlA5aKfCNwfej0XeNQqvlvPAla4+3IAd99dQ7lF6q3c/EP88Z/rmb34S9o2a8wDEwcyeXQP3TJYak04Rd8VyKy0nAWMOdwYdy81s31AB6Av4GY2H0gAZrv77487tUg9VFxaznOfbeHP723gYEkZV4/txS3jk3WqpNS62j4YGwucAowCCoH3zGyJu79XeZCZXQ9cD9CjR49ajiRSt9yd99bs4tdvrWFzbgGn9Uvg7vMGkNRRZ9JI3Qin6LOB7pWWu4XWVTcmKzQv34aKg7JZwEfungtgZm8BI4CvFb27PwU8BZCSkuLH/mWIRKb1Ow/wwBur+XhDLn0SWvDM1aM4rV/HoGNJAxNO0S8Gks2sFxWFPgm4tMqYVGAKsAC4CHjf3b+asrndzJoDxcC3gYdrKrxIpNpbWMzD76znhYVf0qJJDPeeP4ArvtVT96SRQBy16ENz7lOB+VScXjnD3dPNbBqQ5u6pwHRgppllAHlUvBng7nvM7E9UvFk48Ja7v1lLX4tI4ErKypn1+VYefncDB4pKuGxMT356Zl/at2gSdDRpwMw9smZKUlJSPC0tLegYIsfsw/U5PPDGajJ25XNyUgfuPX8g/Tq3CjqWNBCh458p1W3TlbEix2lTTj6/fnMN763dRc8OzXnqipGcOaCTzoeXiKGiF/mG9h0s4S/vbeC5BVuIi43hznP6c9XJicTF6tbBEllU9CLHqKzcmb34S/74z/XsKSzmkpTu3HpWPxJaxQUdTaRaKnqRY7Bg426mvbGaNdv3M7pXe+49fwCDurYJOpbIEanoRcKQvfcg//PGat5etYNu7ZrpvjRSr6joRY6guLSc6Z9s5s/vbQDgtrP6ct2pvWnaWPPwUn+o6EUOY8HG3dzz+ioyduVz9sBO3PvdgXRt2yzoWCLHTEUvUsWuA0X89q21vPZFNt3bN2PGVSmc3r9T0LFEvjEVvUhIWbkza+FWHpq/jkMl5dx8ehI3nZakaRqp91T0IsCyzL38ct5KVmXv55SkeKZNHKjntErUUNFLg7a3sJiH5q/jb4u+JKFlHI9eOpzzBp+gs2kkqqjopUFyd+YuyeK3b69l38ESrjm5Fz85I5lWTfUQEIk+KnppcNbu2M8981axeMseRvZsxwMTBzGgS+ugY4nUGhW9NBj5h0p55N31zPh0C62bxvL7C4dw0chuNGqkaRqJbip6iXruzlsrd/DAG6vZsb+IyaN7cPvZ/Wine8RLA6Gil6i2ObeAe19fxccbchlwQmsev3wEI3q0CzqWSJ1S0UtUKiop4/F/beTJf20kLrYR9393AJef1JNYPcpPGiAVvUSdD9fncM+8VXyZV8jEYV24+9wT6di6adCxRAKjopeosa+whAfeXM3cJVn0TmjB364bw9ik+KBjiQRORS9R4b01O7nrtZXk5hfz49P6cPP4ZD3pSSQkrAlLM5tgZuvMLMPM7qhme5yZzQltX2hmiaH1iWZ20MyWhf48WcP5pYHbW1jMz+Ys49rn0mjXvAnzbjqZn5/dXyUvUslR9+jNLAZ4DDgTyAIWm1mqu6+uNOxaYI+7J5nZJOBB4JLQto3uPqxmY4vA/PQd/HLeKvYUFHPz+GSmnpZEk1gdbBWpKpypm9FAhrtvAjCz2cBEoHLRTwTuD72eCzxqulmI1JK8gmLuT00ndfk2TjyhNc9ePYqBXfQ4P5HDCafouwKZlZazgDGHG+PupWa2D+gQ2tbLzL4A9gO/dPePq34CM7seuB6gR48ex/QFSMPy9srt3PP6KvYdLOGnZ/TlptP60FinTIocUW0fjN0O9HD33WY2EphnZgPdfX/lQe7+FPAUQEpKitdyJqmHcvMPcd/r6by5cjuDurbmhevG0L+z7k8jEo5wij4b6F5puVtoXXVjsswsFmgD7HZ3Bw4BuPsSM9sI9AXSjje4NAzuzhsrtnNfajr5RaX8/Ox+XD+ut/biRY5BOEW/GEg2s15UFPok4NIqY1KBKcAC4CLgfXd3M0sA8ty9zMx6A8nAphpLL1Et58Ah7pm3in+k72BotzY89IOh9O3UKuhYIvXOUYs+NOc+FZgPxAAz3D3dzKYBae6eCkwHZppZBpBHxZsBwDhgmpmVAOXADe6eVxtfiEQPdyd1+TbuS02nsLiMO87pz3Wn9NLtC0S+IauYXYkcKSkpnpammZ2Gatf+Iu56bRXvrtnJ8B5teeiioSR11CP9RI7GzJa4e0p123RlrEQEd+fVpdn86u/pHCot55fnncjVJ/ciRveKFzluKnoJ3I59Rdz56go+WJdDSs92/P6iIXowt0gNUtFLYNydl5dk8cAbqykpK+fe8wcwZWyi9uJFapiKXgKxr7CE2+Yu553VOxndqz2/v3AIifEtgo4lEpVU9FLnlmfu5cd/W8rO/UX88rwTuebkXnpuq0gtUtFLnXF3nv1sC795aw0dWzXl5RvGMqx726BjiUQ9Fb3Uif1FJfxi7greXrWDM07syB9+MJS2zfVwbpG6oKKXWrcqex83zVpK9t6D3H3uiVx3ai90c1ORuqOil1rj7rzw+VYeeGMNHVo24aUfncTInu2DjiXS4KjopVYcKCrhzldX8saK7XynXwJ/ungY7VtoqkYkCCp6qXGrt+3nx39bypd5hdw+oR83jOujs2pEAqSilxrj7sxenMl9qem0a96YF394EqN7aapGJGgqeqkRBYdKufu1lcxbto1Tk+N5+JJhxLeMCzqWiKCilxqwbscBbpq1hM25Bdx6Zl9+fFqSpmpEIoiKXo7LS2mZ3Pv6Klo1bcwL141hbJ/4oCOJSBUqevlGCotLuWdeOq8szWJsnw48Mmk4Ca00VSMSiVT0cswydh3gxheWkpGTzy3jk7l5fLLuOCkSwVT0ckxeXZrF3a+tokVcDDOvGcMpyZqqEYl0KnoJS1FJGfe9ns6ctEzG9GrPnycPp1PrpkHHEpEwqOjlqDLzCvnh82ms3XGAqacl8ZMzkvWgbpF6JKyfVjObYGbrzCzDzO6oZnucmc0JbV9oZolVtvcws3wzu62GcksdWb1tP//1xGds23uQZ68exW1n91PJi9QzR/2JNbMY4DHgHGAAMNnMBlQZdi2wx92TgIeBB6ts/xPw9vHHlbr02cZcLvnrAmIbGXNvHMt3+nUMOpKIfAPh7JqNBjLcfZO7FwOzgYlVxkwEngu9nguMt9B9aM3s+8BmIL1GEkudeHPFdq6asZjObZryyo1j6dupVdCRROQbCqfouwKZlZazQuuqHePupcA+oIOZtQR+AfzqSJ/AzK43szQzS8vJyQk3u9SS5z7bwtQXlzK0exvm3jCWLm2bBR1JRI5DbU+23g887O75Rxrk7k+5e4q7pyQkJNRyJDkcd+eh+Wu5LzWdM0/sxMxrx9CmeeOgY4nIcQrnrJtsoHul5W6hddWNyTKzWKANsBsYA1xkZr8H2gLlZlbk7o8eb3CpWSVl5dz56krmLsni0jE9eGDiIF0EJRIlwin6xUCymfWiotAnAZdWGZMKTAEWABcB77u7A6d+NcDM7gfyVfKRp7C4lB/PWsoH63L4yRnJ3DI+WY/6E4kiRy16dy81s6nAfCAGmOHu6WY2DUhz91RgOjDTzDKAPCreDKQeyCso5ppnF7Miay+/uWAwl47pEXQkEalhVrHjHTlSUlI8LS0t6BgNQmZeIVNmLCJ770H+PHk4Zw/sHHQkEfmGzGyJu6dUt01XxjZQq7ftZ8ozizhUUsYL141hVKKeBCUSrVT0DdBnG3P50fNLaNk0llk6R14k6qnoG5g3V2znp3OW0bNDc567ZrTOkRdpAFT0Dcizn27mV2+sJqVnO/7vyhTaNm8SdCQRqQMq+gag4kKodTz+r42cOaATf5k8nKaNY4KOJSJ1REUf5SpfCDV5dA8emDhQd58UaWBU9FFMF0KJCKjoo1ZeQTFXP7uYlVl7+fUFg7hsTM+gI4lIQFT0UeirC6Gy9h7kictH6kIokQZORR9lKl8INUsXQokIKvqosnrbfi55agEt43QhlIj8fyr6KLF1dwFXzlhEy7hYXr7hW3Rr1zzoSCISIXSeXRTYtb+IK6YvorS8nJnXjlbJi8jXqOjruX0HS7hyxiJy8w/xzFWjSOqo6RoR+ToVfT12sLiM655bzMacfJ68fCTDe7QLOpKIRCDN0ddTJWXlTP3bUtK27uHPk4Yzrq+etSsi1dMefT1UXu784pUVvLd2F9MmDuK7Q7sEHUlEIpiKvp5xd37z1hpeXZrNT8/oyxUn6YpXETkyFX0988SHG3n6k81M+VZPbh6fFHQcEakHwip6M5tgZuvMLMPM7qhme5yZzQltX2hmiaH1o81sWejPcjO7oIbzNygvLvqS3/9jHd8b2oX7vjtQNygTkbActejNLAZ4DDgHGABMNrMBVYZdC+xx9yTgYeDB0PpVQIq7DwMmAH81Mx0A/gb+sWo7d7+2knF9E/jDD4bSqJFKXkTCE84e/Wggw903uXsxMBuYWGXMROC50Ou5wHgzM3cvdPfS0PqmgNdE6Ibms4xcbn5xGUO7t+XJy0fQJFYzbiISvnAaoyuQWWk5K7Su2jGhYt8HdAAwszFmlg6sBG6oVPz/ZmbXm1mamaXl5OQc+1cRxVZm7eOHz6eRGN+cZ64aRfMm+oVIRI5Nre8auvtCdx8IjALuNLOm1Yx5yt1T3D0lIUHng39lY04+U55ZRNvmTXj+mjF6xquIfCPhFH020L3ScrfQumrHhObg2wC7Kw9w9zVAPjDom4ZtSLbvO8iV0xdhwMxrR9O5zX+8P4qIhCWcol8MJJtZLzNrAkwCUquMSQWmhF5fBLzv7h76b2IBzKwn0B/YUiPJo9iegmKunL6IfQdLePbq0fROaBl0JBGpx4464evupWY2FZgPxAAz3D3dzKYBae6eCkwHZppZBpBHxZsBwCnAHWZWApQDN7l7bm18IdGisLiUa55bzNbdhTx7zSgGd2sTdCQRqefMPbJOhElJSfG0tLSgYwSiuLSc655P45MNOTx+2UgmDNIjAEUkPGa2xN1TqtumUzgiRHm5c+vLy/lofQ6/+6/BKnkRqTE6ITsCuDv3/z2dvy/fxu0T+jFpdI+gI4lIFFHRR4BH3tvA8wu28sNTe3Hjt/sEHUdEooyKPmAzF2zhf9/dwIUjunHXuSfq/jUiUuNU9AFKXb6Ne1PTOePEjjx44WCVvIjUChV9QD5an8PP5ixjVM/2PHrpCGJj9L9CRGqH2iUAm3ML+PGspSR1bMn/TUmhaeOYoCOJSBRT0dexwuJSbnxhCTExxv9dmUKbZo2DjiQiUU7n0dchd+fu11axbucBnrlqFN3bNw86kog0ANqjr0MzP9/Ka19UPOv1O/06Bh1HRBoIFX0dWbI1j2l/X83p/Tsy9TQ961VE6o6Kvg7kHDjETbOW0qVtMx6+eJgeAygidUpz9LWstKyc/35xKXsLS3j1plG0aa6DryJSt1T0teyh+ev4fFMef/zBUAZ20S2HRaTuaeqmFr29cjt//WgTl5/UgwtHdgs6jog0UCr6WpKxK5/bXl7O0O5tuef8AUHHEZEGTEVfCwoOlXLDC0uIaxzDE5eNIC5WV76KSHBU9DXM3bn9lRVsysnnL5OH06Vts6AjiUgDp6KvYdM/2cybK7bz87P7c3JSfNBxRERU9DVp4abd/PbttZw1oBM3fLt30HFERIAwi97MJpjZOjPLMLM7qtkeZ2ZzQtsXmlliaP2ZZrbEzFaG/j69hvNHjF37i5j64hf0aN+cP1w8VPeWF5GIcdSiN7MY4DHgHGAAMNnMqp5Gci2wx92TgIeBB0Prc4HvuvtgYAows6aCR5KSsnJumrWU/KJSnrx8JK2b6qIoEYkc4ezRjwYy3H2TuxcDs4GJVcZMBJ4LvZ4LjDczc/cv3H1baH060MzM4moieCT5zVtrSNu6h99dOJh+nVsFHUdE5GvCKfquQGal5azQumrHuHspsA/oUGXMhcBSdz9U9ROY2fVmlmZmaTk5OeFmjwivL8vmmU+3cPXJiUwcVvWfRUQkeHVyMNbMBlIxnfOj6ra7+1PunuLuKQkJCXURqUas23GAO15ZSUrPdtx17olBxxERqVY4RZ8NdK+03C20rtoxZhYLtAF2h5a7Aa8BV7r7xuMNHCn2F5Vw4wtLaBEXy2OXjaCxnvkqIhEqnHZaDCSbWS8zawJMAlKrjEml4mArwEXA++7uZtYWeBO4w90/raHMgXN3bntpOVvzCnns0uF0at006EgiIod11KIPzblPBeYDa4CX3D3dzKaZ2fdCw6YDHcwsA/gZ8NUpmFOBJOBeM1sW+lPvH6305Ieb+Ofqndx5Tn/G9K56KEJEJLKYuwed4WtSUlI8LS0t6BiH9WlGLldMX8g5g0/g0cnDdb68iEQEM1vi7inVbdPE8jHYtvcgN7/4Bb0TWvL7C4eo5EWkXlDRh+lQaRk3zVpKUUkZT14+khZxemaLiNQPaqswPfDGapZl7uWJy0aQ1LFl0HFERMKmPfowvLIkixc+/5IfjevNOYNPCDqOiMgxUdEfxept+7nrtZWc1Ls9Pz+7X9BxRESOmYr+CA4Wl3Hz7C9o06wxf5k8glhdFCUi9ZDm6I/gt2+vIWNXPjOvHU1Cq6i7F5uINBDaRT2MD9bu4vkFW7nm5F6cmlx/7r8jIlKVir4aufmH+Pnc5fTv3IrbJ2heXkTqN03dVOHu/GLuCvYXlfLCdWNo2jgm6EgiIsdFe/RVzFr4Je+t3cUvJvSnf+fWQccRETluKvpKMnbl8z9vrubU5HiuHpsYdBwRkRqhog8pLi3nJ3O+oFnjGP7wg6E0aqT72IhIdNAcfcjD765nVfZ+nrx8pO4vLyJRRXv0wOebdvPkhxuZNKo7EwZ1DjqOiEiNavBFv+9gCT+bs4ye7Ztzz/kDgo4jIlLjGvTUjbvzy3mr2HngEK/cOFa3HhaRqNSg9+jnLcvm78u38ZPxyQzr3jboOCIitaLBFn1mXiH3zksnpWc7bjotKeg4IiK1pkEWfVm587OXluHAw5cMI0anUopIFAur6M1sgpmtM7MMM7ujmu1xZjYntH2hmSWG1ncwsw/MLN/MHq3h7N/YE//KYPGWPUybOJDu7ZsHHUdEpFYdtejNLAZ4DDgHGABMNrOqp6dcC+xx9yTgYeDB0Poi4B7gthpLfJyWZ+7lf9/dwPlDTuCC4V2DjiMiUuvC2aMfDWS4+yZ3LwZmAxOrjJkIPBd6PRcYb2bm7gXu/gkVhR+4gkOl/GTOMjq2iuPX3x+MmaZsRCT6hVP0XYHMSstZoXXVjnH3UmAf0CHcEGZ2vZmlmVlaTk5OuP/ZMfufN1ezZXcBf7x4GG2aN661zyMiEkki4mCsuz/l7inunpKQUDsP+ZifvoMXF2Vy/bjefKtP2O9BIiL1XjhFnw10r7TcLbSu2jFmFgu0AXbXRMCasGt/EXe8soKBXVpz65l6kIiINCzhFP1iINnMeplZE2ASkFplTCowJfT6IuB9d/eai/nNlZc7t81dQWFxGY9MGkaT2Ij4JUZEpM4c9Zp/dy81s6nAfCAGmOHu6WY2DUhz91RgOjDTzDKAPCreDAAwsy1Aa6CJmX0fOMvdV9f4V3IYzy3Ywkfrc3hg4kCSOraqq08rIhIxwrq5i7u/BbxVZd29lV4XAT84zH+beBz5jsu6HQf47dtrOb1/Ry4/qWdQMUREAhW18xiHSsu4ZfYXtIqL5cELh+hUShFpsKL2do0P/WMda3ccYPqUFBJaxQUdR0QkMFG5R//Jhlye/mQzl5/Ug/Endgo6johIoKKu6PcUFHPry8vok9CCu8/Vg0RERKJq6sbdueu1leQVFDN9yiiaNYkJOpKISOCiao/+5SVZvL1qBz87sx+DurYJOo6ISESImqLfklvAr1LTOal3e64f1zvoOCIiESNqit4MRvRsx58u1oNEREQqi5o5+p4dWjDz2jFBxxARiThRs0cvIiLVU9GLiEQ5Fb2ISJRT0YuIRDkVvYhIlFPRi4hEORW9iEiUU9GLiEQ5i5BHu/6bmeUAW4/jQ8QDuTUUp7bVp6xQv/Iqa+2pT3nrU1Y4vrw93T2hug0RV/THy8zS3D0l6BzhqE9ZoX7lVdbaU5/y1qesUHt5NXUjIhLlVPQiIlEuGov+qaADHIP6lBXqV15lrT31KW99ygq1lDfq5uhFROTronGPXkREKlHRi4hEuagpejObYGbrzCzDzO4IOs+RmFl3M/vAzFabWbqZ3RJ0pqMxsxgz+8LM3gg6y9GYWVszm2tma81sjZl9K+hMh2NmPw19D6wysxfNrGnQmSozsxlmtsvMVlVa197M3jGzDaG/2wWZ8SuHyfpQ6PtghZm9ZmZtA4z4NdXlrbTtVjNzM4uvic8VFUVvZjHAY8A5wABgspkNCDbVEZUCt7r7AOAk4McRnhfgFmBN0CHC9AjwD3fvDwwlQnObWVfgZiDF3QcBMcCkYFP9h2eBCVXW3QG85+7JwHuh5UjwLP+Z9R1gkLsPAdYDd9Z1qCN4lv/Mi5l1B84CvqypTxQVRQ+MBjLcfZO7FwOzgYkBZzosd9/u7ktDrw9QUURdg011eGbWDTgPeDroLEdjZm2AccB0AHcvdve9gYY6sligmZnFAs2BbQHn+Rp3/wjIq7J6IvBc6PVzwPfrMtPhVJfV3f/p7qWhxc+BbnUe7DAO828L8DBwO1BjZ8pES9F3BTIrLWcRwcVZmZklAsOBhQFHOZL/peIbrzzgHOHoBeQAz4Smmp42sxZBh6qOu2cDf6Biz207sM/d/xlsqrB0cvftodc7gE5BhjkG1wBvBx3iSMxsIpDt7str8uNGS9HXS2bWEngF+Im77w86T3XM7Hxgl7svCTpLmGKBEcAT7j4cKCBypha+JjS3PZGKN6cuQAszuzzYVMfGK87PjvhztM3sbiqmTGcFneVwzKw5cBdwb01/7Ggp+myge6XlbqF1EcvMGlNR8rPc/dWg8xzBycD3zGwLFVNip5vZC8FGOqIsIMvdv/oNaS4VxR+JzgA2u3uOu5cArwJjA84Ujp1mdgJA6O9dAec5IjO7CjgfuMwj+8KhPlS86S8P/bx1A5aaWefj/cDRUvSLgWQz62VmTag4oJUacKbDMjOjYg55jbv/Keg8R+Lud7p7N3dPpOLf9X13j9i9TnffAWSaWb/QqvHA6gAjHcmXwElm1jz0PTGeCD1wXEUqMCX0egrweoBZjsjMJlAx7fg9dy8MOs+RuPtKd+/o7omhn7csYEToe/q4REXRhw62TAXmU/GD8pK7pweb6ohOBq6gYu94WejPuUGHiiL/DcwysxXAMOA3wcapXui3jrnAUmAlFT+PEXXJvpm9CCwA+plZlpldC/wOONPMNlDxW8nvgsz4lcNkfRRoBbwT+jl7MtCQlRwmb+18rsj+TUZERI5XVOzRi4jI4anoRUSinIpeRCTKqehFRKKcil5EJMqp6EVEopyKXkQkyv0/BIv5SXuiOBEAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.plot(vel_profile)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "At the FreeSlip wall the gradient of the velocity should be almost zero"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.testing.assert_almost_equal(np.gradient(vel_profile)[-1], 0)"
+   ]
+  }
+ ],
+ "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.8.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
-- 
GitLab