From ef7c7239d9a6607e85ecb296bd2ee3ee0a690e28 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Wed, 11 Dec 2019 14:33:39 +0100
Subject: [PATCH] Re-added compiled-in boundaries with tests

are much faster in certain setups
---
 lbmpy/boundaries/boundaries_in_kernel.py      | 134 ++++++++++++++
 lbmpy_tests/test_compiled_in_boundaries.ipynb | 175 ++++++++++++++++++
 2 files changed, 309 insertions(+)
 create mode 100644 lbmpy/boundaries/boundaries_in_kernel.py
 create mode 100644 lbmpy_tests/test_compiled_in_boundaries.ipynb

diff --git a/lbmpy/boundaries/boundaries_in_kernel.py b/lbmpy/boundaries/boundaries_in_kernel.py
new file mode 100644
index 00000000..9be20626
--- /dev/null
+++ b/lbmpy/boundaries/boundaries_in_kernel.py
@@ -0,0 +1,134 @@
+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
diff --git a/lbmpy_tests/test_compiled_in_boundaries.ipynb b/lbmpy_tests/test_compiled_in_boundaries.ipynb
new file mode 100644
index 00000000..070c7c2b
--- /dev/null
+++ b/lbmpy_tests/test_compiled_in_boundaries.ipynb
@@ -0,0 +1,175 @@
+{
+ "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
+}
-- 
GitLab