{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lbmpy.session import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 02: Setting up geometry and boundary conditions\n",
    "\n",
    "In this tutorial you will learn how to set up boundaries for your LBM simulation. \n",
    "We begin with simple setup for parametric geometries, then we'll see how to set up spatially and temporally varying velocity boundary conditions.\n",
    "\n",
    "## Geometry Setup\n",
    "\n",
    "Lets start with a 3D simulation of a pipe. The flow should be along the x axis, with a circular cross section \n",
    "in the y-z plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "domain_size = (64, 16, 16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc1 = LatticeBoltzmannStep(domain_size=domain_size, method='srt', relaxation_rate=1.9, force=(1e-6, 0, 0),\n",
    "                           periodicity=(True, False, False), optimization={'target': 'cpu'})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For now the flow should be driven by an external force which we already specified. Additionally we have enabled perdiodicity in x direction. Now the circular pipe needs to be created. Therefor we use a mask callback function that gets x,y,z coordiante arrays that have the same shape as the domain. The callback function has to return a boolean array, where the boundary is set in cells that are marked with True."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def pipe_geometry_callback(x, y, z):\n",
    "    radius = domain_size[1] / 2\n",
    "    y_mid = domain_size[1] / 2\n",
    "    z_mid = domain_size[2] / 2\n",
    "    return (y - y_mid) ** 2 + (z - z_mid) ** 2 > radius ** 2\n",
    "\n",
    "\n",
    "wall = NoSlip()\n",
    "sc1.boundary_handling.set_boundary(wall, mask_callback=pipe_geometry_callback)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we plot a cross section through our domain to see if the boundary setup worked. A slice has to be passed in to specify which plane of the 3D domain we want to plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/oAAAFlCAYAAABbSZMeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAbTUlEQVR4nO3df7DddX3n8dc7iQmiiCBX0JgQIxIFbMpyzYi6Wn+WItUqywjWFa2aKeM2stZWjTM4dGYztut0tbtd1xRQnEVsR7G11LplXQvutEovQoSI/NiUQigxFxArkEZDPvtHLjaNCbn3nJPc5MPjMcPc+/2e7/l+38l87x2e+Z7zPdVaCwAAANCHObM9AAAAADA6Qh8AAAA6IvQBAACgI0IfAAAAOiL0AQAAoCNCHwAAADoyb38e7KijjmpLlizZn4cEHocmJydnewSAA87Y2NhsjwA8Dlx33XX3ttb8wpll+zX0lyxZkomJif15SOBxaO3atbM9AsABZ+XKlbM9AvA4UFX/MNsz4KX7AAAA0BWhDwAAAB0R+gAAANARoQ8AAAAdEfoAAADQEaEPAAAAHRH6AAAA0BGhDwAAAB0R+gAAANARoQ8AAAAdEfoAAADQEaEPAAAAHRH6AAAA0BGhDwAAAB0R+gAAANCRebM9AACs2bBwtkdghlYvvXu2RwAA9sAVfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjgh9AAAA6MheQ7+qLqmqzVV10y7rf6Oqbqmq9VX1e/tuRAAAAGC6pnNF/zNJTtt5RVW9Iskbkvxca+3EJB8b/WgAAADATO019Ftr1yS5f5fV5yX5aGtt69Q2m/fBbAAAAMAMDfoe/eOT/Nuq+lZVXV1VLxzlUAAAAMBg5g3xvCOSvCjJC5P8SVUtba21XTesqpVJVibJ4sWLB50TAAAAmIZBr+hvTHJF2+HaJNuTHLW7DVtra1tr46218bGxsUHnBAAAAKZh0ND/0ySvTJKqOj7J/CT3jmooAAAAYDB7fel+VV2e5BeSHFVVG5N8JMklSS6Z+si9Hyc5d3cv2wcAAAD2r72GfmvtnD089NYRzwIAAAAMadCX7gMAAAAHIKEPAAAAHRH6AAAA0JG9vkcfgP6s2bBwtkfgIHegnUOrl9492yMAwAHDFX0AAADoiNAHAACAjgh9AAAA6IjQBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjgh9AAAA6IjQBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjgh9AAAA6IjQBwAAgI4IfQAAAOiI0AcAAICOzJvtAQCYnjUbFs72CHDAGtXPx+qld49kPwAwm1zRBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjuw19KvqkqraXFU37eax91dVq6qj9s14AAAAwExM54r+Z5KctuvKqlqU5DVJ7hzxTAAAAMCA9hr6rbVrkty/m4f+S5LfTtJGPRQAAAAwmIHeo19Vr09yd2tt3TS2XVlVE1U1MTk5OcjhAAAAgGmacehX1aFJPpzkguls31pb21obb62Nj42NzfRwAAAAwAwMckX/OUmenWRdVd2R5FlJvl1Vx4xyMAAAAGDm5s30Ca21G5M8/dHlqdgfb63dO8K5AAAAgAFM5+P1Lk/yt0mWVdXGqnrnvh8LAAAAGMRer+i31s7Zy+NLRjYNAAAAMJSB7roPAAAAHJiEPgAAAHRE6AMAAEBHZnzXfQBmZs2GhbM9AjBNo/x5Xb307pHtCwBmwhV9AAAA6IjQBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjgh9AAAA6IjQBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjgh9AAAA6IjQBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjgh9AAAA6IjQBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoyF5Dv6ouqarNVXXTTuv+c1V9r6q+U1Vfqqqn7tsxAQAAgOmYzhX9zyQ5bZd1VyU5qbX2c0luTfKhEc8FAAAADGCvod9auybJ/bus+6vW2rapxW8medY+mA0AAACYoVG8R//Xkvzlnh6sqpVVNVFVE5OTkyM4HAAAALAnQ4V+VX04ybYkl+1pm9ba2tbaeGttfGxsbJjDAQAAAHsxb9AnVtW5Sc5I8qrWWhvdSAAAAMCgBgr9qjotyQeSvLy19vBoRwIAAAAGNZ2P17s8yd8mWVZVG6vqnUn+W5LDklxVVTdU1f/Yx3MCAAAA07DXK/qttXN2s/rifTALAAAAMKRR3HUfAAAAOEAIfQAAAOiI0AcAAICODPzxegC9W7Nh4WyPABzERvU7ZPXSu0eyHwAeP1zRBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjgh9AAAA6IjQBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjgh9AAAA6IjQBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjgh9AAAA6IjQBwAAgI4IfQAAAOiI0AcAAICOCH0AAADoiNAHAACAjuw19KvqkqraXFU37bTuyKq6qqpum/p6xL4dEwAAAJiO6VzR/0yS03ZZ98EkX2utPTfJ16aWAQAAgFm219BvrV2T5P5dVr8hyaVT31+a5FdGPBcAAAAwgEHfo390a+2eJJn6+vTRjQQAAAAMap/fjK+qVlbVRFVNTE5O7uvDAQAAwOPaoKH//ap6RpJMfd28pw1ba2tba+OttfGxsbEBDwcAAABMx6Ch/+Uk5059f26SPxvNOAAAAMAwpvPxepcn+dsky6pqY1W9M8lHk7ymqm5L8pqpZQAAAGCWzdvbBq21c/bw0KtGPAsAAAAwpH1+Mz4AAABg/xH6AAAA0BGhDwAAAB0R+gAAANARoQ8AAAAdEfoAAADQEaEPAAAAHRH6AAAA0BGhDwAAAB0R+gAAANARoQ8AAAAdEfoAAADQEaEPAAAAHRH6AAAA0BGhDwAAAB0R+gAAANARoQ8AAAAdEfoAAADQEaEPAAAAHRH6AAAA0BGhDwAAAB0R+gAAANCRebM9AAAAAOzsuuuue/q8efMuSnJSXKDe1fYkN23btu1dp5xyyubdbSD0AQAAOKDMmzfvomOOOeb5Y2NjP5gzZ06b7XkOJNu3b6/JyckTNm3adFGS1+9uG/8yAgAAwIHmpLGxsX8S+T9rzpw5bWxs7IfZ8WqH3W+zH+cBAACA6Zgj8vds6u9mjz0v9AEAAOAxvO9973vmBRdccPS+2v/LX/7y4+699965o9rfUO/Rr6r/mORdSVqSG5O8o7X2z6MYDAAAAJLk53/nr5Y/8PBPRnaPuace+oRtN1zw2nWj2t+wrr766ttHub+Br+hX1cIkq5KMt9ZOSjI3ydmjGgwAAACSZJSRP939feADHzhmyZIlJ734xS8+/rbbbluQJH/zN3/zxOXLlz/v+OOPP+E1r3nNcyYnJ+cmyYoVK5a9853vXDQ+Pr5s6dKlJ1599dWHvva1r33Osccee9KqVaue+eg+X/3qVz/nxBNPfP5xxx134sc+9rGjHl2/cOHCF9xzzz3zbrnllvlLly498eyzzz72uOOOO/ElL3nJcx988MGa6Z9v2Jfuz0vyxKqal+TQJP845P4AAABgVn3jG9849Etf+tKRN95443evvPLK29etW/ekJHn729/+7DVr1my89dZbv3viiSdu+cAHPvDTiJ8/f/72iYmJW97xjndMnnXWWcf90R/90Z3f+9731v/xH//xUZs2bZqbJJdddtkd69evv/mGG2747qc+9amjH12/szvvvPOQVatWbb799tvXH3744Y989rOfPWKm8w8c+q21u5N8LMmdSe5J8sPW2l8Nuj8AAAA4EHz9619/8umnn/7AYYcdtv3II4/c/trXvvaBhx56aM6PfvSjua973eseTJJ3v/vd933zm9988qPPeeMb3/hAkixfvnzLcccdt+XYY4/9yROf+MS2aNGirRs2bJifJL/7u7979LJly0445ZRTnr9p06YnrF+//pBdj71w4cKtL37xi7ckycknn/zwHXfcsWCm8w/z0v0jkrwhybOTPDPJk6rqrbvZbmVVTVTVxOTk5KCHAwAAgP2mamavmD/kkENaksyZMycLFiz46ScGzJkzJ9u2basrr7zysKuvvvqwiYmJ791yyy3fff7zn79ly5YtP9Pk8+fP/+lz586d27Zt27ZfX7r/6iR/31qbbK39JMkVSV6860attbWttfHW2vjY2NgQhwMAAIB975WvfOWDf/EXf/HUBx98sH7wgx/Mueqqq576pCc9aftTnvKUR7761a8+OUkuvvjip5166qkPTnefDzzwwNzDDz/8kcMOO2z79ddff8ijbwfYF4a5ocGdSV5UVYcm2ZLkVUkmRjIVAAAAzJKXvvSlD7/xjW+8/6STTjpx4cKFW1esWPFgknz605/++/POO+/YVatWzVm8ePHWyy+//I7p7vPMM8/84dq1a8eOP/74E57znOf88/Llyx/aV/NXa23vW+3pyVUXJnlzkm1Jrk/yrtba1j1tPz4+3iYm/FsAsG+tXbt2JPtZs2HhSPYDMIzVS+8eyX5Wrlw5kv0APJaquq61Nj7sftatW3fH8uXL7310ufeP1xvEunXrjlq+fPmS3T021F9Ua+0jST4yzD4AAADgsRzsUb6/DfvxegAAAMABROgDAABAR4Q+AAAAdEToAwAAQEeEPgAAAHRE6AMAAMAuquqUd7/73c96dPmCCy44+n3ve98zH+s569atW7BixYplz3ve805YunTpieecc86xSXLllVce9opXvOK4JLnssssOX7169TH7cvaRfQ4hAAAA7AuXXnrp8q1bt46sXxcsWLDt3HPPfcyP7Js/f377yle+csQ999yz6RnPeMa26ez3Pe95z+JVq1Z9/61vfesDSXLttdc+cddtfvVXf/WHSX440ODT5Io+AAAAB7RRRv509zd37tz2tre9bXLNmjVH7/rYrbfeOv/UU089/vjjjz/h1FNPPf62226bnySbN29+wrHHHvvjR7dbsWLFll2f+wd/8AdPe9vb3rY4Sc4888wlb3nLWxafcsopy5YsWXLS5Zdffvhwf7IdhD4AAADsxm/91m9tvuKKK46877775u68/td//dcXv+Utb7nv1ltv/e6b3/zm+84777xFSfKe97zn+6effvrxL3vZy5574YUXPv3ee++du/s9/4u77rprwbXXXnvLn//5n992/vnnH/vwww/XsHMLfQAAANiNI488cvtZZ51130c/+tGn77z++uuvf9LKlSvvT5Lzzjvv/uuuu+7JSfLe9773vhtvvHH9m970pvuvueaaw174whc+b8uWLY8Z7meeeeb9c+fOzQte8IKtixYt2nrDDTccMuzcQh8AAAD24EMf+tD3P/e5zx310EMPTauflyxZ8pPzzz//vq997Wv/b968eZmYmPiZ9+nvrKoec3kQQh8AAAD24Oijj37kl3/5l3/wuc997qhH15188skPXXTRRUckyac+9akjx8fHH0ySL3zhC0/ZunVrJcmdd94574EHHpi783v2d+eKK6444pFHHsn69esX3HXXXQuWL1/+z8PO7K77AAAA8Bg+/OEPb7r00kvHHl3+5Cc/eee555675BOf+MQxT3va07Z99rOfvSNJvvrVrz7l/e9//+IFCxZsT5ILL7xw4+LFi7d95zvf2eO+jzvuuK0rVqxYdt999z3h4x//+D8ceuihbdh5hT4AAAAHtAULFmwb9cfr7W2bhx9++PpHv1+0aNG2LVu2/HR52bJlP/7mN795667PueiiizYm2bjr+jPOOONHZ5xxxo+SZNWqVfclue/Rx1760pc+ePHFF9818z/Fngl9AAAADmh7+8x7/jWhDwAAALPgi1/84h37Yr9uxgcAAAAdEfoAAAAcaLZv3759+M+Z69TU3832PT0u9AEAADjQ3DQ5OXm42P9Z27dvr8nJycOT3LSnbbxHHwAAgAPKtm3b3rVp06aLNm3adFJcoN7V9iQ3bdu27V172kDoAwAAcEA55ZRTNid5/WzPcbDyLyMAAADQEaEPAAAAHRH6AAAA0BGhDwAAAB0R+gAAANARoQ8AAAAdEfoAAADQkaFCv6qeWlVfqKrvVdXNVXXqqAYDAAAAZm7ekM//RJKvttb+XVXNT3LoCGYCAAAABjRw6FfVU5K8LMnbk6S19uMkPx7NWAAAAMAghnnp/tIkk0k+XVXXV9VFVfWkXTeqqpVVNVFVE5OTk0McDgAAANibYUJ/XpJ/k+STrbWTkzyU5IO7btRaW9taG2+tjY+NjQ1xOAAAAGBvhgn9jUk2tta+NbX8hewIfwAAAGCWDBz6rbVNSe6qqmVTq16V5LsjmQoAAAAYyLB33f+NJJdN3XF/Q5J3DD8SAAAAMKihQr+1dkOS8RHNAgAAAAxpmPfoAwAAAAcYoQ8AAAAdEfoAAADQkWFvxgfQrdVL7x7JftZsWDiS/QAHl1H9DgGAmXJFHwAAADoi9AEAAKAjQh8AAAA6IvQBAACgI0IfAAAAOiL0AQAAoCNCHwAAADoi9AEAAKAjQh8AAAA6IvQBAACgI0IfAAAAOiL0AQAAoCNCHwAAADoi9AEAAKAjQh8AAAA6IvQBAACgI0IfAAAAOiL0AQAAoCNCHwAAADoi9AEAAKAjQh8AAAA6IvQBAACgI0IfAAAAOjJ06FfV3Kq6vqquHMVAAAAAwOBGcUX/vUluHsF+AAAAgCENFfpV9awkr0ty0WjGAQAAAIYx7BX9jyf57STb97RBVa2sqomqmpicnBzycAAAAMBjGTj0q+qMJJtba9c91nattbWttfHW2vjY2NighwMAAACmYZgr+i9J8vqquiPJ55O8sqr+50imAgAAAAYycOi31j7UWntWa21JkrOT/J/W2ltHNhkAAAAwY6O46z4AAABwgJg3ip201v46yV+PYl8AAADA4FzRBwAAgI4IfQAAAOiI0AcAAICOjOQ9+gDs2eqld49kP2s2LBzJfoA9G9XPKwDMJlf0AQAAoCNCHwAAADoi9AEAAKAjQh8AAAA6IvQBAACgI0IfAAAAOiL0AQAAoCNCHwAAADoi9AEAAKAjQh8AAAA6IvQBAACgI0IfAAAAOiL0AQAAoCNCHwAAADoi9AEAAKAjQh8AAAA6IvQBAACgI0IfAAAAOjJvtgcAYHpWL717ZPtas2HhyPYFB4JR/nwAwMHOFX0AAADoiNAHAACAjgh9AAAA6IjQBwAAgI4MHPpVtaiqvl5VN1fV+qp67ygHAwAAAGZumLvub0vym621b1fVYUmuq6qrWmvfHdFsAAAAwAwNfEW/tXZPa+3bU9//KMnNSXxeEwAAAMyikbxHv6qWJDk5ybdGsT8AAABgMEOHflU9OckXk5zfWvun3Ty+sqomqmpicnJy2MMBAAAAj2Go0K+qJ2RH5F/WWrtid9u01ta21sZba+NjY2PDHA4AAADYi2Huul9JLk5yc2vt90c3EgAAADCoYa7ovyTJv0/yyqq6Yeq/00c0FwAAADCAgT9er7X2f5PUCGcBAAAAhjSSu+4DAAAABwahDwAAAB0R+gAAANARoQ8AAAAdGfhmfAAcvFYvvXu2R/hX1mxYONsjMEMH2jkEAPwLV/QBAACgI0IfAAAAOiL0AQAAoCNCHwAAADoi9AEAAKAjQh8AAAA6IvQBAACgI0IfAAAAOiL0AQAAoCNCHwAAADoi9AEAAKAjQh8AAAA6IvQBAACgI0IfAAAAOiL0AQAAoCNCHwAAADoi9AEAAKAj82Z7AABYvfTu2R4BAKAbrugDAABAR4Q+AAAAdEToAwAAQEeEPgAAAHRkqNCvqtOq6paqur2qPjiqoQAAAIDBDBz6VTU3yR8m+aUkJyQ5p6pOGNVgAAAAwMwNc0V/RZLbW2sbWms/TvL5JG8YzVgAAADAIIYJ/YVJ7tppeePUOgAAAGCWDBP6tZt17Wc2qlpZVRNVNTE5OTnE4QAAAIC9GSb0NyZZtNPys5L8464btdbWttbGW2vjY2NjQxwOAAAA2JthQv/vkjy3qp5dVfOTnJ3ky6MZCwAAABjEvEGf2FrbVlX/Icn/SjI3ySWttfUjmwwAAACYsYFDP0laa19J8pURzQIAAAAMaZiX7gMAAAAHGKEPAAAAHRH6AAAA0BGhDwAAAB2p1tr+O1jVZJJ/2G8HZJSOSnLvbA8BI+ScpkfOa3rjnKY3j4dz+tjW2thsD/F4t19Dn4NXVU201sZnew4YFec0PXJe0xvnNL1xTrO/eOk+AAAAdEToAwAAQEeEPtO1drYHgBFzTtMj5zW9cU7TG+c0+4X36AMAAEBHXNEHAACAjgh9HlNVnVVV66tqe1WN7/LYh6rq9qq6pap+cbZmhJmqqtOmztvbq+qDsz0PzFRVXVJVm6vqpp3WHVlVV1XVbVNfj5jNGWEmqmpRVX29qm6e+v+O906td15zUKqqQ6rq2qpaN3VOXzi13jnNfiH02ZubkrwpyTU7r6yqE5KcneTEJKcl+e9VNXf/jwczM3We/mGSX0pyQpJzps5nOJh8Jjt+9+7sg0m+1lp7bpKvTS3DwWJbkt9srT0/yYuSvGfqd7PzmoPV1iSvbK0tT/LzSU6rqhfFOc1+IvR5TK21m1trt+zmoTck+XxrbWtr7e+T3J5kxf6dDgayIsntrbUNrbUfJ/l8dpzPcNBorV2T5P5dVr8hyaVT31+a5Ff261AwhNbaPa21b099/6MkNydZGOc1B6m2w4NTi0+Y+q/FOc1+IvQZ1MIkd+20vHFqHRzonLv06ujW2j3JjmhK8vRZngcGUlVLkpyc5FtxXnMQq6q5VXVDks1JrmqtOafZb+bN9gDMvqr630mO2c1DH26t/dmenrabdT7CgYOBcxfgAFVVT07yxSTnt9b+qWp3v7Lh4NBaeyTJz1fVU5N8qapOmu2ZePwQ+qS19uoBnrYxyaKdlp+V5B9HMxHsU85devX9qnpGa+2eqnpGdlxBgoNGVT0hOyL/stbaFVOrndcc9FprD1TVX2fHvVWc0+wXXrrPoL6c5OyqWlBVz07y3CTXzvJMMB1/l+S5VfXsqpqfHTeV/PIszwSj8OUk5059f26SPb0iCw44tePS/cVJbm6t/f5ODzmvOShV1djUlfxU1ROTvDrJ9+KcZj+p1rxilT2rqjcm+a9JxpI8kOSG1tovTj324SS/lh13yj2/tfaXszYozEBVnZ7k40nmJrmktfafZnkkmJGqujzJLyQ5Ksn3k3wkyZ8m+ZMki5PcmeSs1tquN+yDA1JVvTTJN5LcmGT71OrV2fE+fec1B52q+rnsuNne3Oy4uPonrbXfqaqnxTnNfiD0AQAAoCNeug8AAAAdEfoAAADQEaEPAAAAHRH6AAAA0BGhDwAAAB0R+gAAANARoQ8AAAAdEfoAAADQkf8PB4v6pTw+KfoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.boundary_handling(sc1.boundary_handling, make_slice[0.5, :, :])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the setup is complete. We can run it and look at the velocity profile in the pipe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA54AAAFlCAYAAACDRTcUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAXGUlEQVR4nO3df6zd5X0f8PfHvnjDPxIgXMAGFhsEJAhNSXSLlmVqq6YMp6tKK5UJpkZ0jcT+SLd0WtUmQ1ryz6Rq67pOW9TNa1gylhCNNGmjqiGhXSs2iWa5UFIIDoQBSYwdfBEJmB+dMX72h2866tjY3PM8Pvcev16Sde/5nnPf34/13O85fvt7flRrLQAAADDKumkPAAAAwGxTPAEAABhK8QQAAGAoxRMAAIChFE8AAACGUjwBAAAYau5U7uzcc89t27dvP5W7BGCV+MZ9j/cPPeOMvnlVffNG6P0xaC+/3DcvyWXv2NE9E4C14d577326tTZ/9PZTWjy3b9+excXFU7lLAFaJaze+t3vmuq3nd81rZ5zSh8UVqZcPdc07vO+prnlJ8sXF27pnArA2VNU3j7XdU20BAAAYSvEEAABgKMUTAACAoRRPAAAAhlI8AQAAGErxBAAAYCjFEwAAgKEUTwAAAIZSPAEAABhK8QQAAGAoxRMAAIChFE8AAACGUjwBAAAYSvEEAABgKMUTAACAoRRPAAAAhlI8AQAAGGpu2gMAMLlr1l3fPXP9lZd3zXv2urd1zUuS53b0/f/Tg29sXfNG2PBsdc17w+Nbu+Ylyc6rbuma98pDj3TNS5K7Dt/RPROA43PGEwAAgKEUTwAAAIZSPAEAABhK8QQAAGAoxRMAAIChFE8AAACGUjwBAAAY6oTFs6purar9VfXgMa775apqVXXumPEAAABY607mjOfHk+w8emNVXZzkmiTf6jwTAAAAM+SExbO1dneSZ45x1b9N8itJWu+hAAAAmB0reo1nVf1Ukidba189idveXFWLVbW4tLS0kt0BAACwhr3u4llVG5PckuRfnMztW2u7WmsLrbWF+fn517s7AAAA1riVnPG8NMmOJF+tqieSXJTkvqq6oOdgAAAAzIa51/sDrbUHkpz3/cvL5XOhtfZ0x7kAAACYESfzcSq3J7knyRVVtaeq3jd+LAAAAGbFCc94ttZuPMH127tNAwAAwMxZ0bvaAgAAwMlSPAEAABhK8QQAAGAoxRMAAIChXvfHqQAwuZ1b398176Xrru6alyTfvaLvQ8TzbznYNS9Jtm17pmve/JkvdM0bYemlTV3z9l5+Tte8JDnw5nO75p192Vld85L+x+Cd+z7aNQ9g1jjjCQAAwFCKJwAAAEMpngAAAAyleAIAADCU4gkAAMBQiicAAABDKZ4AAAAMpXgCAAAwlOIJAADAUIonAAAAQymeAAAADKV4AgAAMJTiCQAAwFCKJwAAAEMpngAAAAyleAIAADCU4gkAAMBQiicAAABDKZ4AAAAMNTftAQBWu2vWXd8986Xrru6a9+SP9P9/xDdc+kzXvPdse6JrXpK8deO+7pmnm91v2to98543be+a9+QFZ3XNS5ILs71r3oj7ibsO39E9E2BanPEEAABgKMUTAACAoRRPAAAAhlI8AQAAGErxBAAAYCjFEwAAgKEUTwAAAIY6YfGsqluran9VPfiqbf+6qr5eVX9eVZ+rqv4fsAUAAMBMOJkznh9PsvOobXcluaq19jeTPJLkQ53nAgAAYEacsHi21u5O8sxR277UWju0fPFPk1w0YDYAAABmQI/XeP5Cki8c78qqurmqFqtqcWlpqcPuAAAAWEsmKp5VdUuSQ0k+ebzbtNZ2tdYWWmsL8/Pzk+wOAACANWhupT9YVTcl+ckk726ttX4jAQAAMEtWVDyrameSX03yI621F/uOBAAAwCw5mY9TuT3JPUmuqKo9VfW+JP8hyZYkd1XV/VX1HwfPCQAAwBp1wjOerbUbj7H5YwNmAQAAYAb1eFdbAAAAOC7FEwAAgKEUTwAAAIZSPAEAABhqxZ/jCXC6WH/l5d0zv3tF37vfN1z6TNe8JPkHl36la9783IGueUlyyYb9XfM21cGueSO80DZ0zduy/qWueUmy49Klrnmfyg91zUuS737nnK55mwfcTwDMEmc8AQAAGErxBAAAYCjFEwAAgKEUTwAAAIZSPAEAABhK8QQAAGAoxRMAAIChFE8AAACGUjwBAAAYSvEEAABgKMUTAACAoRRPAAAAhlI8AQAAGErxBAAAYCjFEwAAgKEUTwAAAIZSPAEAABhK8QQAAGAoxRMAAICh5qY9AEBv1258b9e8Z697W9e8JHn+LQe75r1n2xNd85Jkfu5A17zLNnyna16SXDL3Yte8LetW/8PigcPPd81bn8Nd80Z454Df7y+8ZXPXvO9985yueUn/+7Ivvnhb1zyA18MZTwAAAIZSPAEAABhK8QQAAGAoxRMAAIChFE8AAACGUjwBAAAY6oTFs6purar9VfXgq7adU1V3VdU3lr+ePXZMAAAA1qqTOeP58SQ7j9r2wSR/1Fq7LMkfLV8GAACAH3DC4tlauzvJM0dtvi7JJ5a//0SSn+48FwAAADNipa/xPL+1ti9Jlr+e128kAAAAZsnwNxeqqpurarGqFpeWlkbvDgAAgFVmpcXzqaramiTLX/cf74attV2ttYXW2sL8/PwKdwcAAMBatdLi+fkkNy1/f1OS3+szDgAAALPmZD5O5fYk9yS5oqr2VNX7kvxakmuq6htJrlm+DAAAAD9g7kQ3aK3deJyr3t15FgAAAGbQ8DcXAgAA4PSmeAIAADCU4gkAAMBQiicAAABDKZ4AAAAMdcJ3tQVYa9ZtPb9r3nM7+v8f3bZtz3TNe+vGfV3zkuSSDfv75s292DUvSbasO/0exnr/nUesyyvp+7tzYOOZXfOS5KvbtnXNe3bHBV3zkuTszvdlANPkjCcAAABDKZ4AAAAMpXgCAAAwlOIJAADAUIonAAAAQymeAAAADKV4AgAAMJTiCQAAwFCKJwAAAEMpngAAAAyleAIAADCU4gkAAMBQiicAAABDKZ4AAAAMpXgCAAAwlOIJAADAUIonAAAAQymeAAAADKV4AgAAMNTctAcA6K2d0feu7eAbW9e8JJk/84Xumb1tqoNd87asW/0POZvX/fVpj3BCzx/+i655I9al9+/OCL2PwaUB9xO978sApskZTwAAAIZSPAEAABhK8QQAAGAoxRMAAIChFE8AAACGUjwBAAAYaqLiWVX/tKq+VlUPVtXtVbX634ceAACAU2rFxbOqLkzyT5IstNauSrI+yQ29BgMAAGA2TPpU27kkZ1bVXJKNSfZOPhIAAACzZMXFs7X2ZJJfT/KtJPuSPNta+1KvwQAAAJgNkzzV9uwk1yXZkWRbkk1V9XPHuN3NVbVYVYtLS0srnxQAAIA1aZKn2v54ksdba0uttZeTfDbJ3z76Rq21Xa21hdbawvz8/AS7AwAAYC2apHh+K8nfqqqNVVVJ3p1kd5+xAAAAmBWTvMbzy0k+k+S+JA8sZ+3qNBcAAAAzYm6SH26tfTjJhzvNAgAAwAya9ONUAAAA4DUpngAAAAyleAIAADCU4gkAAMBQiicAAABDKZ4AAAAMpXgCAAAwlOIJAADAUIonAAAAQymeAAAADKV4AgAAMJTiCQAAwFCKJwAAAEMpngAAAAyleAIAADCU4gkAAMBQiicAAABDKZ4AAAAMpXgCAAAwlOIJAADAUIonAAAAQymeAAAADKV4AgAAMJTiCQAAwFCKJwAAAEPNTXsAgN7q5UNd8zY8W13zkmTppU3dM3t7oW3omnfg8PNd85Jky7q+D2PPH/6LrnlrwYHDfY+XJHmhbe6e2VvvY3DE/UTv+zKAaXLGEwAAgKEUTwAAAIZSPAEAABhK8QQAAGAoxRMAAIChFE8AAACGmqh4VtVZVfWZqvp6Ve2uqnf2GgwAAIDZMOkHoP27JHe21n62qjYk2dhhJgAAAGbIiotnVb0hyQ8n+fkkaa0dTHKwz1gAAADMikmeantJkqUk/6Wq/qyqfruqNh19o6q6uaoWq2pxaWlpgt0BAACwFk1SPOeSvCPJb7XW3p7khSQfPPpGrbVdrbWF1trC/Pz8BLsDAABgLZqkeO5Jsqe19uXly5/JkSIKAAAAf2nFxbO19p0k366qK5Y3vTvJQ12mAgAAYGZM+q62/zjJJ5ff0faxJP9w8pEAAACYJRMVz9ba/UkWOs0CAADADJrkNZ4AAABwQoonAAAAQymeAAAADKV4AgAAMNSk72oLsOoc3vdU17w3PL61a16S7L38nK55u9/Uf8Yt61/qmrc+h7vmJcklcy92zduybvU/LB44fKhr3mOHNnbNS5LHDp7XNW/3iwOOwb19j8Gtj/f//e59XwYwTc54AgAAMJTiCQAAwFCKJwAAAEMpngAAAAyleAIAADCU4gkAAMBQiicAAABDKZ4AAAAMpXgCAAAwlOIJAADAUIonAAAAQymeAAAADKV4AgAAMJTiCQAAwFCKJwAAAEMpngAAAAyleAIAADCU4gkAAMBQiicAAABDKZ4AAAAMNTftAQB6++KLt3XN23nVLV3zkuTAm8/tmnfPm7Z3zUuSHZcudc/s7ZXs75q3qQ52zRvhhba5a95jB8/rmpckS4e2dM27Z+/2rnlJsvnrG7rmnfXA013zkuTOzvdlANPkjCcAAABDKZ4AAAAMpXgCAAAwlOIJAADAUIonAAAAQymeAAAADDVx8ayq9VX1Z1X1+z0GAgAAYLb0OOP5gSS7O+QAAAAwgyYqnlV1UZK/l+S3+4wDAADArJn0jOdvJvmVJIePd4OqurmqFqtqcWlpacLdAQAAsNasuHhW1U8m2d9au/e1btda29VaW2itLczPz690dwAAAKxRk5zxfFeSn6qqJ5J8OsmPVdV/6zIVAAAAM2PFxbO19qHW2kWtte1JbkjyP1prP9dtMgAAAGaCz/EEAABgqLkeIa21P0nyJz2yAAAAmC3OeAIAADCU4gkAAMBQiicAAABDKZ4AAAAM1eXNhQBm2SsPPdI98+zLzuqa9+QFffOS5FP5oa5579z2RNe8JDmw8czumaeb3S9u7Z55z97tXfOe+z/9f78vfPhQ17wR9xMAs8QZTwAAAIZSPAEAABhK8QQAAGAoxRMAAIChFE8AAACGUjwBAAAYSvEEAABgKMUTAACAoRRPAAAAhlI8AQAAGErxBAAAYCjFEwAAgKEUTwAAAIZSPAEAABhK8QQAAGAoxRMAAIChFE8AAACGUjwBAAAYSvEEAABgKMUTAACAoeamPQDAanfX4Tu6Z+7c+v6ueRdme9e8JPnud87pmveFt2zumpckX922rWve/JkvdM0bYemlTV3z9u7tu85JsvnrG7rmXfjwoa55SbL5T5/omnfngPsJgFnijCcAAABDKZ4AAAAMpXgCAAAwlOIJAADAUIonAAAAQ624eFbVxVX1x1W1u6q+VlUf6DkYAAAAs2GSj1M5lOSftdbuq6otSe6tqrtaaw91mg0AAIAZsOIznq21fa21+5a/P5Bkd5ILew0GAADAbOjyGs+q2p7k7Um+3CMPAACA2TFx8ayqzUl+J8kvtdaeO8b1N1fVYlUtLi0tTbo7AAAA1piJimdVnZEjpfOTrbXPHus2rbVdrbWF1trC/Pz8JLsDAABgDZrkXW0ryceS7G6t/Ua/kQAAAJglk5zxfFeS9yb5saq6f/nPT3SaCwAAgBmx4o9Taa39ryTVcRYAAABmUJd3tQUAAIDjUTwBAAAYSvEEAABgKMUTAACAoRRPAAAAhlrxu9oCsHJ37vto17xr1l3fNS9JNl95ede8733znK55SfLsjgu65i29sXXNG2HDs33fUH7r44e75iXJWQ883TXvlYce6ZqXJHcevqN7JgDH54wnAAAAQymeAAAADKV4AgAAMJTiCQAAwFCKJwAAAEMpngAAAAyleAIAADCU4gkAAMBQiicAAABDKZ4AAAAMpXgCAAAwlOIJAADAUIonAAAAQymeAAAADKV4AgAAMJTiCQAAwFCKJwAAAEMpngAAAAyleAIAADCU4gkAAMBQc9MeAIDJ3XX4jmmPcELXbnxv98yzt57fNa+dsfofFuvlQ13zDu97qmtektz54m3dMwFY25zxBAAAYCjFEwAAgKEUTwAAAIZSPAEAABhK8QQAAGCoiYpnVe2sqoer6tGq+mCvoQAAAJgdKy6eVbU+yUeTvCfJlUlurKorew0GAADAbJjkjOfVSR5trT3WWjuY5NNJruszFgAAALNikuJ5YZJvv+rynuVtAAAA8JcmKZ51jG3tB25UdXNVLVbV4tLS0gS7AwAAYC2apHjuSXLxqy5flGTv0Tdqre1qrS201hbm5+cn2B0AAABr0STF8ytJLquqHVW1IckNST7fZywAAABmxdxKf7C1dqiqfjHJF5OsT3Jra+1r3SYDAABgJqy4eCZJa+0PkvxBp1kAAACYQZM81RYAAABOSPEEAABgKMUTAACAoRRPAAAAhlI8AQAAGKpaa6duZ1VLSb55ynZ4+jo3ydPTHoIfYF1WH2uyOlmX1cm6rE7WZXWyLquTdTk13txamz964yktnpwaVbXYWluY9hz8VdZl9bEmq5N1WZ2sy+pkXVYn67I6WZfp8lRbAAAAhlI8AQAAGErxnE27pj0Ax2RdVh9rsjpZl9XJuqxO1mV1si6rk3WZIq/xBAAAYChnPAEAABhK8ZwRVXV9VX2tqg5X1cJR132oqh6tqoer6tppzXi6q6qPVNWTVXX/8p+fmPZMp7Oq2rl8TDxaVR+c9jwcUVVPVNUDy8fI4rTnOV1V1a1Vtb+qHnzVtnOq6q6q+sby17OnOePp6Djr4rFlyqrq4qr646ravfxvsQ8sb3fMTMlrrInjZYo81XZGVNVbkxxO8p+S/HJrbXF5+5VJbk9ydZJtSf4wyeWttVemNevpqqo+kuT51tqvT3uW011VrU/ySJJrkuxJ8pUkN7bWHprqYKSqnkiy0FrzOWtTVFU/nOT5JP+1tXbV8rZ/leSZ1tqvLf9nzdmttV+d5pynm+Osy0fisWWqqmprkq2ttfuqakuSe5P8dJKfj2NmKl5jTf5+HC9T44znjGit7W6tPXyMq65L8unW2v9trT2e5NEcKaFwOrs6yaOttcdaaweTfDpHjhUgSWvt7iTPHLX5uiSfWP7+EznyjzhOoeOsC1PWWtvXWrtv+fsDSXYnuTCOmal5jTVhihTP2Xdhkm+/6vKeOPCm6Rer6s+Xny7lKTfT47hYvVqSL1XVvVV187SH4a84v7W2Lznyj7ok5015Hv4/jy2rRFVtT/L2JF+OY2ZVOGpNEsfL1Ciea0hV/WFVPXiMP691pqaOsc3zqwc5wRr9VpJLk7wtyb4k/2aqw57eHBer17taa+9I8p4k719+aiFwfB5bVomq2pzkd5L8UmvtuWnPwzHXxPEyRXPTHoCT11r78RX82J4kF7/q8kVJ9vaZiKOd7BpV1X9O8vuDx+H4HBerVGtt7/LX/VX1uRx5WvTd052KZU9V1dbW2r7l10/tn/ZAJK21p77/vceW6amqM3Kk4HyytfbZ5c2OmSk61po4XqbLGc/Z9/kkN1TVX6uqHUkuS/K/pzzTaWn5Qef7fibJg8e7LcN9JcllVbWjqjYkuSFHjhWmqKo2Lb8JRKpqU5K/G8fJavL5JDctf39Tkt+b4iws89gyfVVVST6WZHdr7TdedZVjZkqOtyaOl+nyrrYzoqp+Jsm/TzKf5HtJ7m+tXbt83S1JfiHJoRx5qsEXpjboaayqbsuRp3a0JE8k+Ufff+0Hp97yW6j/ZpL1SW5trf3LKY902quqS5J8bvniXJJPWZfpqKrbk/xoknOTPJXkw0l+N8l/T/I3knwryfWtNW90cwodZ11+NB5bpqqq/k6S/5nkgRz5hIEk+ec58ppCx8wUvMaa3BjHy9QongAAAAzlqbYAAAAMpXgCAAAwlOIJAADAUIonAAAAQymeAAAADKV4AgAAMJTiCQAAwFCKJwAAAEP9P2wrlSbgvxaiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sc1.run(500)\n",
    "plt.scalar_field(sc1.velocity[domain_size[0] // 2, :, :, 0]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boundary Setup\n",
    "\n",
    "Now instead of creating a periodic, force driven channel, we want to drive the flow with a velocity bounce back boundary condition (UBB) at the inlet and a outflow boundary at the outlet. We want the inflow velocity boundary set up to already prescribe the correct, parabolic flow profile. That means we need a different velocity value at each cell. The outflow boundary condition takes the following form:\n",
    "$$ f_{\\overline{1}jkxyzt} = f_{\\overline{1}jk(x - \\Delta x)yz(t - \\Delta t)} \\left(c \\theta^{\\frac{1}{2}} -u \\right) \\frac{\\Delta t}{\\Delta x} + \\left(1 - \\left(c \\theta^{\\frac{1}{2}} -u \\right) \\frac{\\Delta t}{\\Delta x} \\right) f_{\\overline{1}jk(x - \\Delta x)yzt}$$\n",
    "More information on the outflow condition can be found in appendix F in https://doi.org/10.1016/j.camwa.2015.05.001\n",
    "\n",
    "Again, we set up a scenario, but this time without external force and without periodicity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc2 = LatticeBoltzmannStep(domain_size=domain_size, method='srt', relaxation_rate=1.9,\n",
    "                           optimization={'target': 'cpu'},\n",
    "                           #optimization={'target': 'gpu', 'gpu_indexing_params': {'block_size': (16, 8, 2)}}\n",
    "                           )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now another callback function is required to specify the velocities. This callback gets a boundaryData object, that can be used to get an array of `linkPositions`, `fluidCellPositions` and `boundaryCellPositions` of all boundary links. The object also behaves like a dict where the boundary data can be queried and set by name.\n",
    "Which data one can set, and how its called depends on the concrete boundary condition. The UBB has `vel_0` until `vel_2` data for the x,y and z components of the velocity.\n",
    "\n",
    "This example also demonstrates how to change boundary parameters during simulation. Here we introduced an activate flag: when True a parabolic inflow is set, when False no inflow velocity is prescribed. \n",
    "We'll see later how this is used to turn the inflow on and off. With this technique also time dependent boundary information can be set. Reconfiguring the boundary is costly when running GPU simulations, since the data has to be transferred from CPU to GPU."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def velocity_info_callback(boundary_data, activate=True, **_):\n",
    "    boundary_data['vel_1'] = 0\n",
    "    boundary_data['vel_2'] = 0\n",
    "    \n",
    "    if activate:\n",
    "        u_max = 0.1\n",
    "        y, z = boundary_data.link_positions(1), boundary_data.link_positions(2)\n",
    "        radius = domain_size[1] // 2\n",
    "        centered_y = y - radius\n",
    "        centered_z = z - radius\n",
    "        dist_to_center = np.sqrt(centered_y**2 + centered_z**2)\n",
    "        boundary_data['vel_0'] = u_max * (1 - dist_to_center / radius)\n",
    "    else:\n",
    "        boundary_data['vel_0'] = 0\n",
    "    \n",
    "    \n",
    "inflow = UBB(velocity_info_callback, dim=sc2.method.dim)\n",
    "\n",
    "stencil = get_stencil('D3Q27')\n",
    "outflow = ExtrapolationOutflow(stencil[4], sc2.method)\n",
    "\n",
    "sc2.boundary_handling.set_boundary(inflow, make_slice[0, :, :])\n",
    "sc2.boundary_handling.set_boundary(outflow, make_slice[-1, :, :])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly we can use the callback function from above to set the pipe geometry. This is intentionally done at the end of the setup to overwrite the in- and outflow boundaries in the outermost slices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sc2.boundary_handling.set_boundary(wall, mask_callback=pipe_geometry_callback)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the full setup, a few slices through the domain are plotted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5UAAAHiCAYAAAB4EOiDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzde3hV9Z32//uGSBSJeCCgIgcRAQGLDimtgq2H2sdai7aO47GgY8VDRRyr1Wofrf35ONqnnbHSg1C1aKvUPlZbdapTxqlQrYcJCgIqoIgcBAlQlZPBmM/vj+y0aQwkWXvt7EPer+viyl6HvdZnZef6kDvfdXBECAAAAACAJLrkuwAAAAAAQPEiVAIAAAAAEiNUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASI1Rih2z3sT3H9ibbP7D9Hdu/zHddAIqD7Tts/++dLA/bg3ew7FzbT7dlXQCQJNtftr3S9mbbh9tebvtz+a4L6AwIlZ1QO5rsJEnrJe0REd/IcVkAikimj2y33avZ/HmZADgwIi6KiP8vXzUCKGyZPx4tsL3V9lrbP7W9Zxvf29LvMt+XdGlE9IiIl9KvGMCOECqxMwMkvRIRke9CABSkNyWd2Thh+1BJu+WvHADFwvY3JN0q6SpJPSV9Wg2/d8yy3S3hZgdIWpROhQDag1DZiTWeXmb7+7b/YvtN21/ILJshaaKkb2ZOI/nYyKbt8bYX2X7X9lO2D8nMP8/2o03We932r5tMr7R9WM4PEECu/ULShCbTEyXd2zhhe4btm5pMX2V7je23bf9z0w3Z3sf2I7bft/2CpIN2tFPb5Zm+tcL2O5nTbAmzQJGwvYekGyVNjognIuLDiFgu6Z/UEAzPaaF/HG17Veb1LyT1l/Ro5neUq21vltRV0nzbb7Swz3Lbt2X6z9uZ1+WZZbNtn5p5PS5ztsWJmenP2Z6Xy+8HUAoIlfiUpMWSekn6nqS7bDsizpV0n6TvZU4j+a+mb7I9RNJMSZdLqpT0ezU0926SZks6ynYX2/tJ2kXS2Mz7BknqIenljjg4ADn1nKQ9bB9iu6uk0yW1eN217RMkXSnpeEkHS2r+h6ofS/pA0n6S/jnzb0dulTRE0mGSBkvqK+n65IcBoIMdKWlXSQ81nRkRmyU9roY+sUMR8VVJKyR9KfM7yq0R0SOzeFREtPRHqevUMBp6mKRRksZI+nZm2WxJR2def0bSMkmfbTI9u81HBnRShEq8FRE/i4iPJN2jhl/o+rThfadL+o+ImBURH6rhOobdJB0ZEcskbVJD4/6spP+UtNr2sMz0nyKiPgfHAqDjNY5WHi/pNUmrd7DeP0n6eUQsjIgtkr7TuCATSE+VdH1EbImIhWroRx9j25IukPQvEbExIjZJulnSGSkdD4Dc6yVpfUTUtbBsTWZ52s6W9N2IWBcRNWoYKf1qZtls/X2I/Ncm058VoRJoVVm+C0DerW18ERFbG35fU48dr/5X+0t6q8l7622vVMOIgfS3v/oNzrx+Vw2N+QjRnIFS8gtJcyQdqCanvrZgf0lzm0y/1eR1pRr+P1q5g+Vqtm53SXMz/UqSrIbT3gAUh/WSetkuayFY7pdZnra/+70l83r/zOtnJQ2x3UcNfxAfL+nGzI3IxqihxwHYCUYqkdTbarjuQdJfRw/66W+jFI2h8qjM68a/AvIXP6CERMRbarhhz4lqdipbM2vU0CMa9W/yukZS3U6WN7Ve0jZJIyJiz8y/nk1OfQNQ+J6VVCvpK01n2t5d0hckPSlpixr+gNRo32bbaO9NBP/u9xY19Ji3pYY/qqvhj15TJC2MiO2S/izpCklvREQuQi5QUgiVSOrXkr5o+zjbu0j6hhr+g/hzZvlsScdI2i0iVkn6k6QTJO0jidt8A6XlfEnHZk5r3ZFfSzrX9nDb3SXd0Lggc/r9Q5K+Y7u77eFquOnPx2ROnf+ZpH+33VuSbPe1/b9SOhYAORYR76nh9NOptk+wvYvtgZL+n6RVajgDYp6kE23vbXtfNdzDoal3JA1qx25nSvq27crMCOT1+vtrwGdLulR/+8P3U82mAewEoRKJRMRiSedImqqGkYMvqeGC+e2Z5UskbVZDmFREvK+GC9+fyfwCCaBERMQbEVHdyjqPS7pN0n9Lej3ztalL1XDq/VpJMyT9fCebuzqzjedsvy/pvyQNTVQ8gLyIiO9JulYN92R4X9LzajgF/riIqFVDsJwvabmkP0h6oNkm/lUNIfFd21e2YZc3SapWw40CF0h6MTOv0WxJFfrbqa7NpwHshHkEIQAAAAAgKUYqAQAAAACJESoBAAAAAIkRKgEAAAAAiREqAQAAAACJESoBAAAAAImVdeTOevXqFQMHDkz8/pqamvSKAUpEZWVlVu+fO3fu+ojIbiN5QD8B0pdNP6GXAGjUWX836cw6NFQOHDhQ1dU7fZTZTk2fPj3FaoDSMGnSpKzeb/utlErpUPQTIH3Z9BN6CYBGnfV3k86M018BAAAAAIkRKgEAAAAAiREqAQAAAACJESoBAAAAAIkRKgEAAAAAiREqAQAAAACJESoBAAAAAIkRKgEAAAAAiREqAQAAAACJleW7ALTfuGcuzHcJBeHpsdPyXQJQ9G5e1jffJRSEawetzncJQFGbWj413yUUhMm1k/NdApAXjFQCAAAAABJrNVTavtv2OtsLm82fbHux7UW2v5e7EgGUCvoJgLTQTwCgcLRlpHKGpBOazrB9jKSTJX0iIkZI+n76pQEoQTNEPwGQjhminwBAQWg1VEbEHEkbm82+WNItEVGbWWddDmoDUGLoJwDSQj8BgMKR9JrKIZKOsv287dm2P7mjFW1Psl1tu7qmpibh7gCUMPoJgLS0qZ/QSwAgXUlDZZmkvSR9WtJVkn5t2y2tGBHTI6IqIqoqKysT7g5ACaOfAEhLm/oJvQQA0pU0VK6S9FA0eEFSvaRe6ZUFoBOhnwBIC/0EAPIgaaj8raRjJcn2EEndJK1PqygAnQr9BEBa6CcAkAdlra1ge6akoyX1sr1K0g2S7pZ0d+Y23tslTYyIyGWhAIof/QRAWugnAFA4Wg2VEXHmDhadk3ItAEoc/QRAWugnAFA4Wg2VSG7cMxfmu4SSlqvv79Njp+Vku0A2bl7WN98llLRcfX+vHbQ6J9sFkppaPjXfJZS0XH1/J9dOzsl2gbQkvaYSAAAAAABCJQAAAAAgOUIlAAAAACAxQiUAAAAAIDFCJQAAAAAgMUIlAAAAACAxQiUAAAAAIDFCJQAAAAAgMUIlAAAAACAxQiUAAAAAILGyfBdQbMY9c2G+S0COteczfnrstBxWglJ387K++S4BOdaez/jaQatzWAlK2dTyqfkuATnWns94cu3kHFYCtIyRSgAAAABAYq2GStt3215ne2ELy660HbZ75aY8AKWEfgIgLfQTACgcbRmpnCHphOYzbfeTdLykFSnXBKB0zRD9BEA6Zoh+AgAFodVQGRFzJG1sYdG/S/qmpEi7KACliX4CIC30EwAoHImuqbQ9XtLqiJjfhnUn2a62XV1TU5NkdwBKGP0EQFra2k/oJQCQrnaHStvdJV0n6fq2rB8R0yOiKiKqKisr27s7ACWMfgIgLe3pJ/QSAEhXkpHKgyQdKGm+7eWSDpD0ou190ywMQKdAPwGQFvoJAORJu59TGRELJPVunM407qqIWJ9iXQA6AfoJgLTQTwAgf9rySJGZkp6VNNT2Ktvn574sAKWIfgIgLfQTACgcrY5URsSZrSwfmFo1AEoa/QRAWugnAFA42n36ayka98yF+S4BRao9PztPj52Ww0pQKG5e1jffJaBItedn59pBq3NYCQrB1PKp+S4BRao9PzuTayfnsBJ0JokeKQIAAAAAgESoBAAAAABkgVAJAAAAAEiMUAkAAAAASIxQCQAAAABIjFAJAAAAAEiMUAkAAAAASIxQCQAAAABIjFAJAAAAAEiMUAkAAAAASIxQCQAAAABIjFAJAAAAAEis1VBp+27b62wvbDLv/9p+zfbLth+2vWduywRQCugnANJCPwGAwtGWkcoZkk5oNm+WpJER8QlJSyR9K+W6AJSmGaKfAEjHDNFPAKAgtBoqI2KOpI3N5v0hIuoyk89JOiAHtQEoMfQTAGmhnwBA4Ujjmsp/lvT4jhbanmS72nZ1TU1NCrsDUMLoJwDSssN+Qi8BgHRlFSptXyepTtJ9O1onIqZHRFVEVFVWVmazOwAljH4CIC2t9RN6CQCkqyzpG21PlHSSpOMiItIrCUBnQz8BkBb6CQB0vESh0vYJkq6W9NmI2JpuSQA6E/oJgLTQTwAgP9rySJGZkp6VNNT2KtvnS/qRpApJs2zPs31HjusEUALoJwDSQj8BgMLR6khlRJzZwuy7clALgBJHPwGQFvoJABSOxNdUFrpxz1yY7xKAv9Oen8mnx07LYSVor5uX9c13CcDfac/P5LWDVuewErTH1PKp+S4B+Dvt+ZmcXDs5h5Wg2KXxSBEAAAAAQCdFqAQAAAAAJEaoBAAAAAAkRqgEAAAAACRGqAQAAAAAJEaoBAAAAAAkRqgEAAAAACRGqAQAAAAAJEaoBAAAAAAkRqgEAAAAACRGqAQAAAAAJEaoBAAAAAAk1mqotH237XW2FzaZt7ftWbaXZr7uldsyAZQC+gmAtNBPAKBwtGWkcoakE5rNu0bSkxFxsKQnM9MA0JoZop8ASMcM0U8AoCC0GiojYo6kjc1mnyzpnszreySdknJdAEoQ/QRAWugnAFA4kl5T2Sci1khS5mvvHa1oe5LtatvVNTU1CXcHoITRTwCkpU39hF4CAOnK+Y16ImJ6RFRFRFVlZWWudweghNFPAKSBXgIA6UoaKt+xvZ8kZb6uS68kAJ0M/QRAWugnAJAHSUPlI5ImZl5PlPS7dMoB0AnRTwCkhX4CAHnQlkeKzJT0rKShtlfZPl/SLZKOt71U0vGZaQDYKfoJgLTQTwCgcJS1tkJEnLmDRcelXAuAEkc/AZAW+gkAFI6c36gHAAAAAFC6CJUAAAAAgMQIlQAAAACAxAiVAAAAAIDECJUAAAAAgMQIlQAAAACAxAiVAAAAAIDECJUAAAAAgMQIlQAAAACAxAiVAAAAAIDECJUAAAAAgMQIlQAAAACAxAiVAAAAAIDEsgqVtv/F9iLbC23PtL1rWoUB6FzoJwDSQj8BgI6VOFTa7ivpMklVETFSUldJZ6RVGIDOg34CIC30EwDoeNme/lomaTfbZZK6S3o7+5IAdFL0EwBpoZ8AQAdKHCojYrWk70taIWmNpPci4g/N17M9yXa17eqamprklQIoWfQTAGlpSz+hlwBAurI5/XUvSSdLOlDS/pJ2t31O8/UiYnpEVEVEVWVlZfJKAZQs+gmAtLSln9BLACBd2Zz++jlJb0ZETUR8KOkhSUemUxaAToZ+AiAt9BMA6GDZhMoVkj5tu7ttSzpO0qvplAWgk6GfAEgL/QQAOlg211Q+L+lBSS9KWpDZ1vSU6gLQidBPAKSFfgIAHa8smzdHxA2SbkipFgCdGP0EQFroJwDQsbJ9pAgAAAAAoBMjVAIAAAAAEiNUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASyypU2t7T9oO2X7P9qu0j0ioMQOdCPwGQFvoJAHSssizf/0NJT0TEP9ruJql7CjUB6JzoJwDSQj8BgA6UOFTa3kPSZySdK0kRsV3S9nTKAtCZ0E8ApIV+AgAdL5vTXwdJqpH0c9sv2b7T9u4p1QWgc6GfAEgL/QQAOlg2obJM0j9I+mlEHC5pi6Rrmq9ke5LtatvVNTU1WewOQAmjnwBIS6v9hF4CAOnKJlSukrQqIp7PTD+ohib+dyJiekRURURVZWVlFrsDUMLoJwDS0mo/oZcAQLoSh8qIWCtppe2hmVnHSXollaoAdCr0EwBpoZ8AQMfL9u6vkyXdl7mz2jJJ52VfEoBOin4CIC30EwDoQFmFyoiYJ6kqpVoAdGL0EwBpoZ8AQMfK5ppKAAAAAEAnl+3prwXr6bHT2rzuuGcuzGElQIP2/EyisFw7aHWb1715Wd8cVgI0aM/PJArH5NrJbV53avnUHFYCNGjPzySwM4xUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASI1QCAAAAABLLOlTa7mr7JduPpVEQgM6LfgIgDfQSAOhYaYxUTpH0agrbAQD6CYA00EsAoANlFSptHyDpi5LuTKccAJ0V/QRAGuglANDxsh2pvE3SNyXVp1ALgM6NfgIgDfQSAOhgiUOl7ZMkrYuIua2sN8l2te3qmpqapLsDUMLoJwDSQC8BgPzIZqRyrKTxtpdL+pWkY23/svlKETE9IqoioqqysjKL3QEoYfQTAGmglwBAHiQOlRHxrYg4ICIGSjpD0n9HxDmpVQag06CfAEgDvQQA8oPnVAIAAAAAEitLYyMR8ZSkp9LYFoDOjX4CIA30EgDoOIxUAgAAAAASS2Wkstg9PXZam9cd98yFOawExaY9PzvoHK4dtLrN6968rG8OK0Gxac/PDkrf5NrJbV53avnUHFaCYtOenx0gLYxUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASK8t3AcXm6bHT2rzuuGcuzGElyJX2fMZANq4dtLrN6968rG8OK0GutOczBpKaXDu5zetOLZ+aw0qQK+35jIF8SDxSabuf7T/aftX2IttT0iwMQOdBPwGQFvoJAHS8bEYq6yR9IyJetF0haa7tWRHxSkq1Aeg86CcA0kI/AYAOlnikMiLWRMSLmdebJL0qifOzALQb/QRAWugnANDxUrlRj+2Bkg6X9Hwa2wPQedFPAKSFfgIAHSPrUGm7h6TfSLo8It5vYfkk29W2q2tqarLdHYASRj8BkJad9RN6CQCkK6tQaXsXNTTs+yLioZbWiYjpEVEVEVWVlZXZ7A5ACaOfAEhLa/2EXgIA6crm7q+WdJekVyPi39IrCUBnQz8BkBb6CQB0vGxGKsdK+qqkY23Py/w7MaW6AHQu9BMAaaGfAEAHS/xIkYh4WpJTrAVAJ0U/AZAW+gkAdLxU7v4KAAAAAOicEo9UonVPj52Wk+2Oe+bCnGy32OTq+wsUomsHrc7Jdm9exuP7pNx9f4FCM7l2ck62O7V8ak62W2xy9f0FCh0jlQAAAACAxAiVAAAAAIDECJUAAAAAgMQIlQAAAACAxAiVAAAAAIDECJUAAAAAgMQIlQAAAACAxAiVAAAAAIDECJUAAAAAgMQIlQAAAACAxMryXQDa7+mx0/JdAoASce2g1fkuAUAJmFw7Od8lAMijrEYqbZ9ge7Ht121fk1ZRADof+gmAtNBPAKBjJQ6VtrtK+rGkL0gaLulM28PTKgxA50E/AZAW+gkAdLxsRirHSHo9IpZFxHZJv5J0cjplAehk6CcA0kI/AYAOlk2o7CtpZZPpVZl5ANBe9BMAaaGfAEAHyyZUuoV58bGV7Em2q21X19TUZLE7ACWMfgIgLa32E3oJAKQrm1C5SlK/JtMHSHq7+UoRMT0iqiKiqrKyMovdAShh9BMAaWm1n9BLACBd2YTK/5F0sO0DbXeTdIakR9IpC0AnQz8BkBb6CQB0sMTPqYyIOtuXSvpPSV0l3R0Ri1KrDECnQT8BkBb6CQB0vMShUpIi4veSfp9SLQA6MfoJgLTQTwCgY2Vz+isAAAAAoJNzxMdusJi7ndk1kt5KebO9JK1PeZuFgmMrTsV2bAMioujuVEE/aTeOrTgV07HRS/5eMX127cWxFadiOrai7CedWYeGylywXR0RVfmuIxc4tuJUysdW6kr5s+PYilMpH1upK+XPjmMrTqV8bMg/Tn8FAAAAACRGqAQAAAAAJFYKoXJ6vgvIIY6tOJXysZW6Uv7sOLbiVMrHVupK+bPj2IpTKR8b8qzor6kEAAAAAORPKYxUAgAAAADypKhDpe0TbC+2/brta/JdT5psL7e9wPY829X5ricbtu+2vc72wibz9rY9y/bSzNe98lljUjs4tu/YXp357ObZPjGfNaJt6CfFgX5CPyl09JLiQC+hlyBdRRsqbXeV9GNJX5A0XNKZtofnt6rUHRMRh5XA7Z9nSDqh2bxrJD0ZEQdLejIzXYxm6OPHJkn/nvnsDouI33dwTWgn+klRmSH6CQoUvaSozBC9BEhN0YZKSWMkvR4RyyJiu6RfSTo5zzWhBRExR9LGZrNPlnRP5vU9kk7p0KJSsoNjQ/GhnxQJ+gkKHL2kSNBLgHQVc6jsK2llk+lVmXmlIiT9wfZc25PyXUwO9ImINZKU+do7z/Wk7VLbL2dOQSnK02c6GfpJcaOfoFDQS4obvQRIqJhDpVuYV0q3sh0bEf+ghlNovm77M/kuCG32U0kHSTpM0hpJP8hvOWgD+gkKFf2kuNBLUKjoJcipYg6VqyT1azJ9gKS381RL6iLi7czXdZIeVsMpNaXkHdv7SVLm67o815OaiHgnIj6KiHpJP1PpfXaliH5S3OgnKBT0kuJGLwESKuZQ+T+SDrZ9oO1uks6Q9Eiea0qF7d1tVzS+lvR5SQt3/q6i84ikiZnXEyX9Lo+1pKrxP6SML6v0PrtSRD8pbvQTFAp6SXGjlwAJleW7gKQios72pZL+U1JXSXdHxKI8l5WWPpIeti01fEb3R8QT+S0pOdszJR0tqZftVZJukHSLpF/bPl/SCkmn5a/C5HZwbEfbPkwNpzwtl3Rh3gpEm9BPigf9hH5SyOglxYNeQi9BuhxRSqf6AwAAAAA6UjGf/goAAAAAyDNCJQAAAAAgMUIlAAAAACAxQiUAAAAAIDFCJQAAAAAgMUIlAAAAACAxQiUAAAAAIDFCJQCgJNleZPvofNcBAECpI1QCQCdg+yzb1bY3215j+3Hb4/JdV1psz7B9U9N5ETEiIp7KU0kAAHQahEoAKHG2r5B0m6SbJfWR1F/STySdvIP1yzquOgAAUOwIlQBQwmz3lPRdSV+PiIciYktEfBgRj0bEVZl1vmP7Qdu/tP2+pHNt72/7Edsbbb9u+4Im2xyTGfV83/Y7tv8tM3/XzDY22H7X9v/Y7rODuq62vdr2JtuLbR+Xmd/F9jW238hs59e2927yvnG2/5zZ/krb59qeJOlsSd/MjMQ+mll3ue3PZV6X277N9tuZf7fZLs8sO9r2KtvfsL0uM5J7Xg4+DgAAShKhEgBK2xGSdpX0cCvrnSzpQUl7SrpP0kxJqyTtL+kfJd3cGPwk/VDSDyNiD0kHSfp1Zv5EST0l9ZO0j6SLJG1rviPbQyVdKumTEVEh6X9JWp5ZfJmkUyR9NrPvv0j6ceZ9/SU9LmmqpEpJh0maFxHTMzV/LyJ6RMSXWji+6yR9OvOeUZLGSPp2k+X7ZmrvK+l8ST+2vVcr3zMAACBCJQCUun0krY+IulbWezYifhsR9ZJ6SRon6eqI+CAi5km6U9JXM+t+KGmw7V4RsTkinmsyfx9JgyPio4iYGxHvt7CvjySVSxpue5eIWB4Rb2SWXSjpuohYFRG1kr4j6R8zp+SeLem/ImJmZrR1Q6a2tjhb0ncjYl1E1Ei6scnxNNb+3cx2fy9ps6Shbdw2AACdGqESAErbBkm92nCd5Momr/eXtDEiNjWZ95YaRvGkhpG8IZJey5zielJm/i8k/aekX2VOMf2e7V2a7ygiXpd0uRoC4zrbv7K9f2bxAEkPZ05vfVfSq2oIoX3UMAL6RvPttdH+mWNoejz7N5ne0Cx4b5XUI+G+AADoVAiVAFDanpX0gRpOKd2ZaPL6bUl7265oMq+/pNWSFBFLI+JMSb0l3SrpQdu7Z0b5boyI4ZKOlHSSpAkt7izi/ogYp4YQGZntSA3h9gsRsWeTf7tGxOrMsoPaUH9L3s7sq+nxvN3KewAAQBsQKgGghEXEe5KuV8M1gqfY7m57F9tfsP29HbxnpaQ/S/rXzM13PqGG0cn7JMn2ObYrM6fKvpt520e2j7F9qO2ukt5XwymlHzXfvu2hto/N3CjnAzVcd9m43h2S/o/tAZl1K2033qX2Pkmfs/1Ptsts72P7sMyydyQN2sm3Yqakb2e21yvzPfnlzr97AACgLQiVAFDiIuLfJF2hhhvT1KhhxO9SSb/dydvOlDRQDaN5D0u6ISJmZZadIGmR7c1quGnPGRHxgRpudvOgGgLlq5Jmq+XgVi7pFknrJa1Vw4jntZllP5T0iKQ/2N4k6TlJn8ocxwpJJ0r6hqSNkuap4aY7knSXGq7RfNd2S8d1k6RqSS9LWiDpxcw8AACQJUe0dsYQAAAAAAAtY6QSAAAAAJAYoRIAAAAAkBihEgAAAACQGKESAAAAAJAYoRIAAAAAkFhZR+6sV69eMXDgwA7bX01NTYftCwAAAChFlZWVHbq/uXPnro+Ijt0pstKhoXLgwIGqrq7usP1Nnz69w/YFAAAAlKJJkyZ16P5sv9WhO0TWOP0VAAAAAJAYoRIAAAAAkBihEgAAAACQGKESAAAAAJAYoRIAAAAAkBihEgAAAACQGKESAAAAAJAYoRIAAAAAkBihEgAAAACQWFm+C8ilcc9c2OZ1T9zvsRxWAgAAABSOikOuafO6kzQph5WgFDBSCQAAAABIrNVQaftu2+tsL2w2f7LtxbYX2f5e7koEAAAAABSqtoxUzpB0QtMZto+RdLKkT0TECEnfT780AAAAAEChazVURsQcSRubzb5Y0i0RUZtZZ10OagMAAAAAFLik11QOkXSU7edtz7b9yTSLAgAAAAAUh6R3fy2TtJekT0v6pKRf2x4UEdF8RduTpIZbRvXv3z9pnQAAAACAApR0pHKVpIeiwQuS6iX1amnFiJgeEVURUVVZWZm0TgAAAABAAUoaKn8r6VhJsj1EUjdJ69MqCgAAAABQHFo9/dX2TElHS3OR4tsAABw1SURBVOple5WkGyTdLenuzGNGtkua2NKprwAAAACA0tZqqIyIM3ew6JyUawEAAAAAFJmkp78CAAAAAECoBAAAAAAkR6gEAAAAACRGqAQAAAAAJEaoBAAAAAAkRqgEAAAAACRGqAQAAAAAJEaoBAAAAAAkRqgEAAAAACRGqAQAAAAAJEaoBAAAAAAkRqgEAAAAACRGqAQAAAAAJNZqqLR9t+11the2sOxK22G7V27KAwAAAAAUsraMVM6QdELzmbb7STpe0oqUawIAAAAAFIlWQ2VEzJG0sYVF/y7pm5Ii7aIAAAAAAMUh0TWVtsdLWh0R81OuBwAAAABQRMra+wbb3SVdJ+nzbVx/kqRJktS/f//27g4AAAAAUMCSjFQeJOlASfNtL5d0gKQXbe/b0soRMT0iqiKiqrKyMnmlAAAAAICC0+6RyohYIKl343QmWFZFxPoU6wIAAAAAFIG2PFJkpqRnJQ21vcr2+bkvCwAAAABQDFodqYyIM1tZPjC1agAAAAAARSXR3V8BAAAAAJAIlQAAAACALBAqAQAAAACJESoBAAAAAIkRKgEAAAAAiREqAQAAAACJESoBAAAAAIkRKgEAAAAAiREqAQAAAACJESoBAAAAAIkRKgEAAAAAiREqAQAAAACJESoBAAAAAIm1Gipt3217ne2FTeb9X9uv2X7Z9sO298xtmQAAAACAQtSWkcoZkk5oNm+WpJER8QlJSyR9K+W6AAAAAABFoNVQGRFzJG1sNu8PEVGXmXxO0gE5qA0AAAAAUODSuKbynyU9nsJ2AAAAAABFJqtQafs6SXWS7tvJOpNsV9uurqmpyWZ3AAAAAIACkzhU2p4o6SRJZ0dE7Gi9iJgeEVURUVVZWZl0dwAAAACAAlSW5E22T5B0taTPRsTWdEsCAAAAABSLtjxSZKakZyUNtb3K9vmSfiSpQtIs2/Ns35HjOgEAAAAABajVkcqIOLOF2XfloBYAAAAAQJFJ4+6vAAAAAIBOilAJAAAAAEiMUAkAAAAASIxQCQAAAABIjFAJAAAAAEiMUAkAAAAASIxQCQAAAABIjFAJAAAAAEiMUAkAAAAASIxQCQAAAABIjFAJAAAAAEiMUAkAAAAASIxQCQAAAABIrNVQaftu2+tsL2wyb2/bs2wvzXzdK7dlAgAAAAAKUVtGKmdIOqHZvGskPRkRB0t6MjMNAAAAAOhkWg2VETFH0sZms0+WdE/m9T2STkm5LgAAAABAEUh6TWWfiFgjSZmvvXe0ou1JtqttV9fU1CTcHQAAAACgEOX8Rj0RMT0iqiKiqrKyMte7AwAAAAB0oKSh8h3b+0lS5uu69EoCAAAAABSLpKHyEUkTM68nSvpdOuUAAAAAAIpJWx4pMlPSs5KG2l5l+3xJt0g63vZSScdnpgEAAAAAnUxZaytExJk7WHRcyrUAAAAAAIpMzm/UAwAAAAAoXYRKAAAAAEBihEoAAAAAQGKESgAAAABAYoRKAAAAAEBihEoAAAAAQGKESgAAAABAYoRKAAAAAEBihEoAAAAAQGKESgAAAABAYoRKAAAAAEBiZfkuAAAAAADSMnfu3N5lZWV3ShopBtHSUi9pYV1d3ddGjx69rvlCQiUAAACAklFWVnbnvvvue0hlZeVfunTpEvmupxTU19e7pqZm+Nq1a++UNL758qySu+1/sb3I9kLbM23vms32AAAAACBLIysrK98nUKanS5cuUVlZ+Z4aRn8/vjzphm33lXSZpKqIGCmpq6Qzkm4PAAAAAFLQhUCZvsz3tMX8mO05xmWSdrNdJqm7pLez3B4AAAAAlIwrrrhi/+uvv75Prrb/2c9+dvD69eu75mr7bZH4msqIWG37+5JWSNom6Q8R8Yfm69meJGmSJPXv3z/p7gAAAACg3Q777h9Gvbv1w9TuJbNn913q5l3/+flpbS9bs2fPfj3fNWRz+utekk6WdKCk/SXtbvuc5utFxPSIqIqIqsrKyuSVAgAAAEA7pRko27q9q6++et+BAweOPPLII4csXbq0XJL+/Oc/7zZq1KhhQ4YMGX788ccfVFNT01WSxowZM/T888/vV1VVNXTQoEEjZs+e3f3zn//8QQMGDBh52WWX7d+4zc997nMHjRgx4pDBgweP+P73v9+rcX7fvn0PXbNmTdnixYu7DRo0aMQZZ5wxYPDgwSPGjh178ObNm53mse9INqe/fk7SmxFRExEfSnpI0pHplAUAAAAAxedPf/pT94cffnjvBQsWvPLYY4+9Pn/+/N0l6dxzzz3w5ptvXrVkyZJXRowYse3qq6/+a2Ds1q1bfXV19eLzzjuv5rTTThv8s5/9bMVrr7226IEHHui1du3arpJ03333LV+0aNGr8+bNe2XatGl9Guc3tWLFil0vu+yyda+//vqinj17fnTvvffu1RHHnE2oXCHp07a727ak4yS9mk5ZAAAAAFB8/vjHP/Y48cQT362oqKjfe++96z//+c+/u2XLli6bNm3q+sUvfnGzJF1wwQUbnnvuuR6N7/nyl7/8riSNGjVq2+DBg7cNGDDgw9122y369etXu2zZsm6SdOutt/YZOnTo8NGjRx+ydu3aXRYtWvSxJ2/07du39sgjj9wmSYcffvjW5cuXl3fEMWdzTeXzth+U9KKkOkkvSZqeVmEAAAAAUIwaxtzabtdddw1J6tKli8rLy/9659ouXbqorq7Ojz32WMXs2bMrqqurX6uoqKgfM2bM0G3btn1sgLBbt25/fW/Xrl2jpXVyIaudRMQNETEsIkZGxFcjojatwgAAAACg2Bx77LGb/+M//mPPzZs3+y9/+UuXWbNm7bn77rvX77HHHh898cQTPSTprrvu2ueII47Y3NZtvvvuu1179uz5UUVFRf1LL720a+MptYUi1YtWAQAAAKAzGzdu3NYvf/nLG0eOHDmib9++tWPGjNksST//+c/fvPjiiwdcdtllXfr37187c+bM5W3d5qmnnvre9OnTK4cMGTL8oIMO+mDUqFFbclV/Eo7ouOeCVlVVRXV1dYft75WJbR92PnG/x3JYCQAAAFA4Kg65ps3rLpi4IIeVfJztuRFRlfT98+fPXz5q1Kj1jdOl/kiRjjR//vxeo0aNGth8PiOVAAAAAEpWZw2AHalDLtwEAAAAAJQmQiUAAAAAIDFCJQAAAAAgMUIlAAAAACAxQiUAAAAAIDHu/goAAAAAKVq8eHG3k0466eClS5cuapx3xRVX7N+jR4+PFi1atNtzzz1XUVFR8VFtba2/8pWvbPzBD36wRpLGjBkzdN26dbvsuuuu9du3b/cll1zyzpVXXrl+x3sqDIRKAAAAACVr8dd7jfpo84bUck/XHvvUDf3x+qweU3LTTTetOu+88/6ydetWDxkyZOQFF1ywYdiwYdsl6d577132mc98Zus777zT9eCDDz700ksv3bDrrrtGOtXnBqe/AgAAAChZaQbKtLe3devWLpJUUVFR33zZ+++/33W33XarLysrK+hAKREqAQAAAKBDffvb3z5g2LBhw/v37/+JU045ZWPfvn3rGpdNmDBh0JAhQ4YfeuihI6+88sq3y8oK/+RSQiUAAAAApMj2TuffdNNNq1577bVX1qxZM3/OnDkVs2bN2r1xnXvvvXfZkiVLXlm2bNnLP/rRj/ZdsmRJt46pOrmsQqXtPW0/aPs126/aPiKtwgAAAACgGPXp06fuvffe69p03saNG7v26tWrrum8nj171o8dO3bT7NmzezTfxv777183cuTIrXPmzNm9+bJCk+1I5Q8lPRERwySNkvRq9iUBAAAAQPHq2bNnfe/evT/83e9+VyFJ77zzTtennnqq57HHHru56Xoffvih5s6d22Pw4MG1zbexadOmLosWLeo+dOjQjy0rNIlP0LW9h6TPSDpXkiJiu6Tt6ZQFAAAAAMXrnnvuefOSSy7pf/XVV/eTpKuvvvrtESNG1EoN11Teeuut+3344YceN27c+xMmTHi38X0TJkwY1PhIkTPOOGP9UUcdtTVfx9BW2Vz1OUhSjaSf2x4laa6kKRGxpelKtidJmiRJ/fv3z2J3AAAAANA+XXvsU5f2I0Xast7o0aM/eP7555c0n/+b3/xm+Y7e88ILLyzOorS8yeabWybpHyRNjojnbf9Q0jWS/nfTlSJiuqTpklRVVVXwt8MFAAAAUDqyfaYkWpfNNZWrJK2KiOcz0w+qIWQCAAAAADqJxKEyItZKWml7aGbWcZJeSaUqAAAAAEBRyPbc4smS7rPdTdIySedlXxIAAAAAoFhkFSojYp6kqpRqAQAAAAAUmWyfUwkAAAAA6MQIlQAAAACQojfeeGOX44477qABAwaM7Nev38jzzjuv3wcffODW3nfNNdfs23T6pptu6j1o0KAR48ePP/D222/fZ8KECTl5RuMVV1yxf+/evT8xbNiw4Y3/1q9f37Wt70/teS0AAAAAUGjG/WrcqPdq30st9/Qs71n39BlP7/AxJfX19TrllFMGf+1rX1s3ZcqUN+rq6nTWWWcNmDJlSt9p06at2tm2b7/99v1uueWWtY3Td911V+Xjjz++dNiwYdtvv/32fdI6hpZcdNFF73z3u999J8l7GakEAAAAULLSDJRt2d6jjz5aUV5eXj9lypQNklRWVqY77rhj5QMPPNBr06ZNXZqPOB5zzDGDH3vssYpLLrmkb21tbZdhw4YNHz9+/IFnnXVW/1WrVpWPHz9+8I033ti76T6WLFnS7YgjjhgyZMiQ4UccccSQpUuXdqurq9MBBxxwaH19vdavX9+1S5cuox9//PEekjR69OihCxcuLE/z+9AUoRIAAAAAUrJgwYLdRo0atbXpvL333rt+v/322/7KK6/sMNj95Cc/WV1eXl7/2muvvfLII4+8ef/996/o3bv3h7Nnz15yww03rGu67kUXXdT/rLPO2rBkyZJXTj/99A0XX3xxv7KyMh144IEfvPjii7vOmjWrx/Dhw7c+9dRTPbZt2+a1a9d2GzlyZO3pp58+YM6cOd1b2v8dd9zRp/HU10996lND2nPMnP4KAAAAACmJCNmOHcxPZR8vvfTS7o8//vgbknTxxRdvvPHGGw+QpCOPPHLTk08+WfHmm2+WX3XVVWvuuuuuyjlz5mweNWrUFkl64IEH3trRNjn9FQAAAAAKwKGHHrpt3rx5uzedt3Hjxi5r167tdsghh9SWlZVFfX39X5fV1tamlsmOPvrozU8//XSPF198cffTTjvtvffff7/rk08+WTFu3LhNae2jJYRKAAAAAEjJ+PHjN33wwQddfvSjH+0jSXV1dbrkkkv6nXbaaesrKirqDzrooO2LFi3q/tFHH+n111/f5eWXX/5rAC0rK4va2tpWhzMPP/zwLXfeeedekjRt2rS9q6qqNkvS0UcfveXFF1/s0aVLl+jevXuMGDFi67333lt5zDHHbM7V8UqESgAAAABITZcuXfTb3/729YceemivAQMGjDzwwANHlpeX199+++2rJen444/f3K9fv9qhQ4eOmDJlSr/hw4f/9frLs88+u+aQQw4ZPn78+AN3to+f/vSnK37xi1/0GjJkyPCZM2fu85Of/GSlJO22226x7777bq+qqtoiSUcdddTmLVu2dBkzZsw2SWrrNZXDhg0bvnjx4m5tPWZHfOx035ypqqqK6urqDtvfKxPbfs7yifs9lsNKAAAAgMJRccg1bV53wcQFOazk42zPjYiqpO+fP3/+8lGjRq1vnO7oR4qUsvnz5/caNWrUwObzuVEPAAAAgJLVWQNgR+L0VwAAAABAYoRKAAAAAEBiWYdK211tv2SbixIBAAAAoJNJY6RyiqRXU9gOAAAAAKDIZBUqbR8g6YuS7kynHAAAAABAMcl2pPI2Sd+UVL+jFWxPsl1tu7qmpibL3QEAAABAYbM9+oILLjigcfr666/vc8UVV+y/s/fMnz+/fMyYMUOHDRs2fNCgQSPOPPPMAZL02GOPVRxzzDGDJem+++7ree211+6b2+rbL/EjRWyfJGldRMy1ffSO1ouI6ZKmSw3PqUy6PwAAAABor3vuuWdUbW1tao9SLC8vr5s4ceJOH1PSrVu3+P3vf7/XmjVr1u633351bdnu17/+9f6XXXbZO+ecc867kvTCCy/s1nyds88++z1J7yUqPIeyGakcK2m87eWSfiXpWNu/TKUqAAAAAEhBmoGyrdvr2rVrTJgwoebmm2/u03zZkiVLuh1xxBFDhgwZMvyII44YsnTp0m6StG7dul0GDBiwvXG9MWPGbGv+3ttvv32fCRMm9JekU089deBZZ53Vf/To0UMHDhw4cubMmT2zO7LkEofKiPhWRBwQEQMlnSHpvyPinNQqAwAAAIAiddVVV6176KGH9t6wYUPXpvMvuuii/medddaGJUuWvHL66advuPjii/tJ0te//vV3TjzxxCGf+cxnDr7xxht7r1+/vmvLW/6blStXlr/wwguLH3300aWXX375gK1btzpXx7MzPKcSAAAAAFK2995715922mkbbrnllt5N57/00ku7T5o0aaMkXXzxxRvnzp3bQ5KmTJmyYcGCBYu+8pWvbJwzZ07FJz/5yWHbtm3baUg89dRTN3bt2lWHHnpobb9+/WrnzZu3a+6OaMdSCZUR8VREnJTGtgAAAACgFHzrW9965/777++1ZcuWNuWugQMHfnj55ZdvePLJJ98oKytTdXX1x66rbMr2Tqc7CiOVAAAAAJADffr0+ehLX/rSX+6///5ejfMOP/zwLXfeeedekjRt2rS9q6qqNkvSgw8+uEdtba0lacWKFWXvvvtu16bXWLbkoYce2uujjz7SokWLyleuXFk+atSoD3J5PDuS6kWrAAAAAIC/ue6669bec889lY3TP/3pT1dMnDhx4A9/+MN999lnn7p77713uSQ98cQTe1x55ZX9y8vL6yXpxhtvXNW/f/+6l19+eYfbHjx4cO2YMWOGbtiwYZfbbrvtre7du+flaRuESgAAAAAlq7y8vC7tR4q0ts7WrVtfanzdr1+/um3btv11eujQodufe+65Jc3fc+edd66StKr5/JNOOmnTSSedtEmSLrvssg2SNjQuGzdu3Oa77rprZfuPIl2ESgAAAAAlq7VnSiJ7hEoAAAAAKDK/+c1vlue7hkbcqAcAAAAAkBihEgAAAEApqa+vr8/PszVKWOZ7Wt/SMkIlAAAAgFKysKampifBMj319fWuqanpKWlhS8u5phIAAABAyairq/va2rVr71y7du1IMYiWlnpJC+vq6r7W0kJCJQAAAICSMXr06HWSxue7js6E5A4AAAAASIxQCQAAAABIjFAJAAAAAEgscai03c/2H22/anuR7SlpFgYAAAAAKHzZ3KinTtI3IuJF2xWS5tqeFRGvpFQbAAAAAKDAJR6pjIg1EfFi5vUmSa9K6ptWYQAAAACAwpfKNZW2B0o6XNLzLSybZLvadnVNTU0auwMAAAAAFIisQ6XtHpJ+I+nyiHi/+fKImB4RVRFRVVlZme3uAAAAAAAFJKtQaXsXNQTK+yLioXRKAgAAAAAUi2zu/mpJd0l6NSL+Lb2SAAAAAADFIpuRyrGSvirpWNvzMv9OTKkuAAAAAEARSPxIkYh4WpJTrAUAAAAAUGRSufsrAAAAAKBzIlQCAAAAABIjVAIAAAAAEiNUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASI1QCAAAAABIjVAIAAAAAEiNUAgAAAAASI1QCAAAAABIry3cBufT02GltXvdarc5hJQAAAEABqZ2c7wpQQhipBAAAAAAkllWotH2C7cW2X7d9TVpFAQAAAACKQ+JQaburpB9L+oKk4ZLOtD08rcIAAAAAAIUvm5HKMZJej4hlEf9/e/cbqmddx3H8/eHMUVhh5ZTYWW6DUUnklDEWhtgKWSnNJ8GCQHriEwWDQmZPJMGnog8iGLYS+iNj/XGIVMOMepQ7K0VthkPUHWaeEyL+edCYfntwXeF9trPp7nO767q43y843NfvezjcX/hwn5vvuX+/69QJ4EFg52TakiRJkiQNwUqGyrXAsZH1fFtbIsnNSeaSzC0uLq7g6SRJkiRJfbOSoTLL1Oq0QtWeqtpSVVvWrFmzgqeTJEmSJPXNSobKeWDdyHoWOL6ydiRJkiRJQ7KSofIQsCnJhiSrgV3Agcm0JUmSJEkaglXj/mBVnUxyK/AHYAbYW1XPTKwzSZIkSVLvjT1UAlTVI8AjE+pFkiRJkjQwqTrt3jof3JMli8CL5+0Jz+xi4D9dN6FzYmbDZG7DZG7DZG7DZG7DZG4frMuqyjt8Dsh5HSr7IslcVW3pug+9f2Y2TOY2TOY2TOY2TOY2TOYmLbWSG/VIkiRJkqacQ6UkSZIkaWzTOlTu6boBnTMzGyZzGyZzGyZzGyZzGyZzk0ZM5ZlKSZIkSdJkTOsnlZIkSZKkCZiqoTLJjiT/SnI0ye6u+9HykuxNspDk6ZHaJ5IcTPJc+/jxLnvU6ZKsS/JYkiNJnklyW1s3u55K8qEkjyd5ss3sh23dzAYgyUySfyR5uF2bW88leSHJU0meSDLX1syt55JclGR/kmfb97gvmpu01NQMlUlmgB8BXwMuB76V5PJuu9IZ/AzYcUptN/BoVW0CHm3X6peTwPeq6nPANuCW9jVmdv31X2B7VV0BbAZ2JNmGmQ3FbcCRkbW5DcOXq2rzyL+jMLf+uw/4fVV9FriC5nVnbtKIqRkqga3A0ap6vqpOAA8COzvuScuoqr8Ar55S3gk80F4/ANx4XpvSe6qql6vq7+31GzRvumsxu96qxpvt8oL2qzCz3ksyC1wP3D9SNrdhMrceS/Ix4BrgJwBVdaKqXsPcpCWmaahcCxwbWc+3NQ3DpVX1MjTDC3BJx/3oLJKsB64E/obZ9Vq7hfIJYAE4WFVmNgz3ArcD74zUzK3/CvhjksNJbm5r5tZvG4FF4KftdvP7k1yIuUlLTNNQmWVq3vpWmrAkHwF+DXy3ql7vuh+dXVW9XVWbgVlga5LPd92Tzi7JDcBCVR3uuheds6ur6iqaozi3JLmm64b0nlYBVwE/rqorgbdwq6t0mmkaKueBdSPrWeB4R73o3L2S5FMA7eNCx/1oGUkuoBkof1FVv2nLZjcA7XauP9OcZzazfrsa+EaSF2iOcmxP8nPMrfeq6nj7uAD8luZojrn12zww3+7iANhPM2SamzRimobKQ8CmJBuSrAZ2AQc67knv3wHgpvb6JuChDnvRMpKE5szJkaq6Z+RbZtdTSdYkuai9/jDwVeBZzKzXquqOqpqtqvU072V/qqpvY269luTCJB/9/zVwHfA05tZrVfVv4FiSz7SlrwD/xNykJVI1PTtAk3yd5hzKDLC3qu7uuCUtI8mvgGuBi4FXgDuB3wH7gE8DLwHfrKpTb+ajDiX5EvBX4CnePef1A5pzlWbXQ0m+QHODiRmaPzLuq6q7knwSMxuEJNcC36+qG8yt35JspPl0Epotlb+sqrvNrf+SbKa5KdZq4HngO7S/MzE3CZiyoVKSJEmSNFnTtP1VkiRJkjRhDpWSJEmSpLE5VEqSJEmSxuZQKUmSJEkam0OlJEmSJGlsDpWSJEmSpLE5VEqSJEmSxuZQKUmSJEka2/8A4Q6vFYqfBuQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x576 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rc('figure', figsize=(14, 8))\n",
    "\n",
    "plt.subplot(231)\n",
    "plt.boundary_handling(sc2.boundary_handling, make_slice[0, :, :], show_legend=False)\n",
    "plt.title('Inflow')\n",
    "\n",
    "plt.subplot(232)\n",
    "plt.boundary_handling(sc2.boundary_handling, make_slice[0.5, :, :], show_legend=False)\n",
    "plt.title('Middle')\n",
    "\n",
    "plt.subplot(233)\n",
    "plt.boundary_handling(sc2.boundary_handling, make_slice[-1, :, :], show_legend=False)\n",
    "plt.title('Outflow')\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.boundary_handling(sc2.boundary_handling, make_slice[:, 0.5, :], show_legend=True)\n",
    "plt.title('Cross section');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run the channel for a few steps..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvUAAAHSCAYAAABsAzSeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df4xl5Xkn+O/TBQ0YnOAE7BDAA8kyyTLRGjMt7KxXkRM7GWCi9GS1GcEosccabcda2LVXGc2Q7B/OrBQpu8qPiTUWDImZGCUx403iTcuLTJxsrEz+sAN4GAzGzPQSEtoQA/6BMb+arnr2j7qdlIvqrtPue6vuPffzkY7qnnPec95z61TDc5/7nPet7g4AALC49uz2BQAAAKdGUA8AAAtOUA8AAAtOUA8AAAtOUA8AAAtOUA8AAAvutN2+gI3OO++8vuSSS3b7MgAAOAn33nvv0919/m5fx1b+wQ+e3V/68urUz3vv/S/d1d1XT/3E36S5CuovueSS3HPPPbt9GQAAnISq+svdvobj+dKXV/Pnd71+6uddueC/nDf1k56CuQrqAQBgmjrJWtZ2+zJmTk09AAAsOJl6AABGrLPaMvUAAMCck6kHAGC01mvqe7cvY+YE9QAAjJoHZQEAgLknUw8AwGh1Oqs9/vIbmXoAAFhwMvUAAIyaB2UBAGCBdZLVJQjqld8AAMCCk6kHAGDUlqH8RqYeAAAWnEw9AACj1clSDGkpqAcAYNTGP5+s8hsAAFh4MvUAAIxWpw1pCQAAzD+ZegAAxquT1fEn6mXqAQBg0cnUAwAwWp3lGP1GUA8AwIhVVlO7fREzp/wGAAAWnEw9AACj1UnWPCgLAADMO5l6AABGbRlq6gX1AACMVmc5gnrlNwAAsOBk6gEAGLW1lqkHAADmnEw9AACjtSw19YJ6AABGq1NZXYLilPG/QwAAGDlBPQAAo7bWNfVliKq6uqoerqpDVXXTFvurqt4/2X9/VV25Yd//WlUPVtUDVfXhqjrzRH0J6gEAYMqqaiXJB5Jck+TyJNdX1eWbml2T5LLJciDJzZNjL0zyvyTZ193fl2QlyXUn6k9NPQAAo7WLD8peleRQdz+SJFV1R5L9ST63oc3+JLd3dyf5VFWdW1UXTPadluSsqno5yauSPH6izgT1AACMWGW1d6U45cIkj21YP5zkTQPaXNjd91TVLyX5qyQvJPnD7v7DE3Wm/AYAAE7eeVV1z4blwKb9W3090EPaVNVrsp7FvzTJdyY5u6p+8kQXI1MPAMBodZK12eSxn+7ufSfYfzjJxRvWL8orS2iO1+btSf6iu59Kkqr6/ST/bZLfOl5nMvUAADB9dye5rKouraq9WX/Q9eCmNgeTvGMyCs6bkzzT3U9kvezmzVX1qqqqJG9L8tCJOjvloL6qLq6qP6mqhybD7rxnsv3nq+oLVXXfZLn2VPsCAICTtZqa+rKd7j6a5MYkd2U9IP9Idz9YVe+uqndPmt2Z5JEkh5L8epL/aXLsp5P8bpLPJPls1mP2W0/U3zTKb44m+Znu/kxVvTrJvVX1icm+X+3uX5pCHwAAsFC6+86sB+4bt92y4XUnueE4x74vyfuG9nXKQf3kK4InJq+fraqHsv4kLwAA7KruXRv9ZkdN9R1W1SVJ3pjk05NNN05mx7pt8hTvVsccOPbU8FNPPTXNywEAgKylpr7Mm6kF9VV1TpLfS/Le7v5a1mfE+u4kV2Q9k//LWx3X3bd2977u3nf++edP63IAAGBpTGVIy6o6PesB/W939+8nSXd/ccP+X0/ysWn0BQAAQ63PKKv8ZluTYXY+mOSh7v6VDdsv2NDsx5M8cKp9AQAArzSNTP1bkvxUks9W1X2TbT+X5PqquiLrH5AeTfLTU+gLAABOwnI8KDuN0W/+LFtPcXvnFtsAAGDHzHBG2bky/ncIAAAjN5UHZQEAYF6t9vwNQTltMvUAALDgZOoBABitTi3FkJaCegAARm1tCUa/Gf87BACAkZOpBwBgtMwoCwAALASZegAARqtThrQEAADmn0w9AACjtrYEeWxBPQAAo9WdrBrSEgAAmHcy9QAAjFhlLR6UBQAA5pxMPQAAo9VZjpp6QT0AAKNmRlkAAGDuydQDADBancqaGWUBAIB5J1MPAMCoLUNNvaAeAIDR6iRrSzD6zfjfIQAAjJxMPQAAI1ZZNaMsAAAw72TqAQAYLTX1AADAQpCpBwBg1Jahpl5QDwDAaHWX8hsAAGD+ydQDADBqqzL1AADAvJOpBwBgtDrJ2hI8KCtTDwDAiFVWe8/Ul0E9V11dVQ9X1aGqummL/VVV75/sv7+qrpxs/56qum/D8rWqeu+J+pKpBwCAKauqlSQfSPLDSQ4nubuqDnb35zY0uybJZZPlTUluTvKm7n44yRUbzvOFJB89UX+CegAARmt9RtldKb+5Ksmh7n4kSarqjiT7k2wM6vcnub27O8mnqurcqrqgu5/Y0OZtSf6/7v7LE3Wm/AYAAE7eeVV1z4blwKb9FyZ5bMP64cm2k21zXZIPb3cxMvUAAIza6mzy2E93974T7N/q64E+mTZVtTfJjyX52e0uRlAPAMBodWq3ym8OJ7l4w/pFSR4/yTbXJPlMd39xu86U3wAAwPTdneSyqrp0knG/LsnBTW0OJnnHZBScNyd5ZlM9/fUZUHqTTCFTX1UXJ7k9yXckWUtya3f/WlV9W5J/n+SSJI8m+cfd/ZVT7Q8AAE7G2i7ksbv7aFXdmOSuJCtJbuvuB6vq3ZP9tyS5M8m1SQ4leT7Ju44dX1WvyvrIOT89pL9plN8cTfIz3f2Zqnp1knur6hNJ/mmSP+7uX5yMy3lTkn85hf4AAGDudfedWQ/cN267ZcPrTnLDcY59Psm3D+3rlIP6yVcET0xeP1tVD2X9qd39Sd46afahJJ+MoB4AgB3UnazuTk39jprqdxFVdUmSNyb5dJLXHasJmvx87XGOOXBsKKCnnnpqmpcDAABLYWpBfVWdk+T3kry3u7829LjuvrW793X3vvPPP39alwMAAEnWJ5+a9jJvpjKkZVWdnvWA/re7+/cnm794bEasqrogyZPT6AsAAIZaH9Jy/AM+nvI7rKpK8sEkD3X3r2zYdTDJOyev35nkD061LwAA4JWmkal/S5KfSvLZqrpvsu3nkvxiko9U1T9L8ldJfmIKfQEAwElZ3XLi1nGZxug3f5atp7hNkred6vkBAIATm0pNPQAAzKNO5vLB1mkT1AMAMGIelAUAABaATD0AAKO2tgQPysrUAwDAgpOpBwBgtLqTVQ/KAgDAYvOgLAAAMPdk6gEAGK1OLcU49TL1AACw4GTqAQAYNUNaAgAAc0+mHgCA0epkKWrqBfUAAIyaIS0BAIC5J1MPAMB4tSEtAQCABSBTDwDAaHWWY0hLQT0AAKOm/AYAAJh7MvUAAIzWsoxTL1MPAAALTqYeAIBRW4ZMvaAeAIDR6hinHgAAWAAy9QAAjNoyjFMvUw8AAAtOph4AgPHq5XhQVqYeAAAWnEw9AACjtSyTTwnqAQAYtWUI6pXfAADAghPUAwAwWscmn5r2MkRVXV1VD1fVoaq6aYv9VVXvn+y/v6qu3LDv3Kr63ar6fFU9VFXff6K+BPUAADBlVbWS5ANJrklyeZLrq+ryTc2uSXLZZDmQ5OYN+34tyce7+3uTvCHJQyfqT009AACj1rtTU39VkkPd/UiSVNUdSfYn+dyGNvuT3N7dneRTk+z8BUmeS/IDSf5pknT3kSRHTtSZoB4AgFHbpRllL0zy2Ib1w0neNKDNhUmOJnkqyb+rqjckuTfJe7r7ueN1pvwGAABO3nlVdc+G5cCm/Vt9kuiBbU5LcmWSm7v7jVnP3L+iJn+jucrU/+eHHs+P/P2f37bd0W89Y9D5Tnvu5W3bvHjemYPOtWd18z3Y2tBvd3rP9g3X9g472dEzh302Wx14vjOeXd22zdmPfn3QueqJpwe162eHna9Xt7+2rA27VwDA+PXsZpR9urv3nWD/4SQXb1i/KMnjA9t0ksPd/enJ9t/NNkG9TD0AAEzf3Ukuq6pLq2pvkuuSHNzU5mCSd0xGwXlzkme6+4nu/uskj1XV90zavS3fWIv/ClPJ1FfVbUl+NMmT3f19k20/n+R/zHo9UJL8XHffOY3+AABgqN14ULa7j1bVjUnuSrKS5LbufrCq3j3Zf0uSO5Ncm+RQkueTvGvDKf7nJL89+UDwyKZ9rzCt8pvfTPJvkty+afuvdvcvTakPAAA4ScPHlZ+2SUL7zk3bbtnwupPccJxj70tyovKebzCV8pvu/tMkX57GuQAAgJMz65r6GyezY91WVa+ZcV8AAPAK3TX1Zd7MMqi/Ocl3J7kiyRNJfnmrRlV14NhQQC8ffX6GlwMAAOM0syEtu/uLx15X1a8n+dhx2t2a5NYk+Zazv9NYhAAATE1nZkNazpWZZeonU9we8+NJHphVXwAAsMymNaTlh5O8Neszax1O8r4kb62qK7L+AenRJD89jb4AAGCwXp+AauymEtR39/VbbP7gyZ7n5XNW8sQPfOu27b7+d9YGne/Mp7afLfal84adq14eOLvruQNmPE2S0wf0O7DPOjrsL7WODmqWVz2+/Z/FnpfOHnSus595blC7fv6FQe0yZEZZAIAN1qL8BgAAmHMze1AWAAB2W2d3ZpTdaTL1AACw4GTqAQAYsVqKIS0F9QAAjNoyjH6j/AYAABacTD0AAKPmQVkAAGDuydQDADBa3cuRqZ+roH5tb/Lc6wfMtPralwad7+WXtp9Rtl837FxrXzpjULs6c+CMp89t/6s//avDvkjZM3Dm2d4z7CmRPS8PaTToVOmVYQ2rBr6HYd0CAPyNZRj9RvkNAAAsuLnK1AMAwLQZ0hIAAJh7MvUAAIyaB2UBAGCBdWopgnrlNwAAsOBk6gEAGLUleE5Wph4AABbdfGXqT19LvmP7yaDOOuvIoNO9eNqAyaeODvtcc/YXhrV76YW9g9rt/cr2tV1nfnnY58quYe1ePG9YPdnRs7Zv8/LZK4POlZWB7faMv9YNANgFSzKjrEw9AAAsuPnK1AMAwLQtQVG9oB4AgFFTfgMAAMw9mXoAAEatl6D8RqYeAAAWnEw9AACj1VmOmnpBPQAA49VJliCoV34DAAALbq4y9a8+46X84H/1n7dtd7SHfRb55Je+d9s23/kdXxl0rqe/8LpB7V7+1tVB7fa8tP2vfu2ZQafKnpeHtRv4a0sPmAR2ddjEuekzhv2JVfl8CQDMhgdlAQCAuTdXmXoAAJi6JcjUC+oBABixWorRb5TfAADAghPUAwAwbj2DZYCqurqqHq6qQ1V10xb7q6reP9l/f1VduWHfo1X12aq6r6ru2a4v5TcAADBlVbWS5ANJfjjJ4SR3V9XB7v7chmbXJLlssrwpyc2Tn8f8YHc/PaQ/mXoAAMar12eUnfYywFVJDnX3I919JMkdSfZvarM/ye297lNJzq2qC76ZtymoBwCAk3deVd2zYTmwaf+FSR7bsH54sm1om07yh1V17xbnfgXlNwAAjNtshrR8urv3nWD/Vun8zVdyojZv6e7Hq+q1ST5RVZ/v7j89XmdzFdS/7vSv5b2v+6Nt2z3Xwy77gae3//biNWe+MOhcT68Narb1rdnCyovbtzn9+YHnOjLsL/WF1WEXt3rG9uc7esawc/XpA6anTVIrvjQCAGZlV4a0PJzk4g3rFyV5fGib7j7288mq+mjWy3mOG9SLpAAAYPruTnJZVV1aVXuTXJfk4KY2B5O8YzIKzpuTPNPdT1TV2VX16iSpqrOT/EiSB07U2VQy9VV1W5IfTfJkd3/fZNu3Jfn3SS5J8miSf9zdX5lGfwAAMNguzCjb3Uer6sYkdyVZSXJbdz9YVe+e7L8lyZ1Jrk1yKMnzSd41Ofx1ST5aVcl6vP473f3xE/U3rfKb30zyb5LcvmHbTUn+uLt/cTIu501J/uWU+gMAgLnW3XdmPXDfuO2WDa87yQ1bHPdIkjecTF9TKb+ZFO1/edPm/Uk+NHn9oST/aBp9AQDASdmlyad20iwflH1ddz+RJJPaoNfOsC8AAHilTjJsXPmFtusPylbVgWPje37ly0OHmAEAAI6ZZVD/xWMzYk1+PrlVo+6+tbv3dfe+13zbrn/GAABgZLqnv8ybWUbRB5O8c/L6nUn+YIZ9AQDA0prWkJYfTvLWrE+XezjJ+5L8YpKPVNU/S/JXSX5iGn0BAMBJmcPM+rRNJajv7uuPs+ttJ3Oes2pP/t7es7Zt9xcvf33Q+c4546Vt27xw9PRB59r77MDZWM8cNoPq3me3b3P688OeMdjz8rC/1NNeGPbFzNrp27/XWhvWZ60O/Fe06nkKAGBGPCgLAADMu1kOaQkAALuulqD8RqYeAAAWnEw9AADjNaczwE6bTD0AACw4mXoAAEaslmL0G0E9AADjpvwGAACYdzL1AACM2xJk6hcyqH+xh33B8MLL288We/beI8M6HTjh6crA0628uP1f156jQ2dtHdbn2sC7vXrm9v2unjGwzzOHdbqyd9jMvvXC9jVxvQz/cgEANljIoB4AAAZbgnyfoB4AgPHqLMXoNx6UBQCABSdTDwDAqNUSlN/I1AMAwIKTqQcAYNxk6gEAgHknqAcAgAWn/AYAgFFbhgdl5yqo/9pa8scvrGzb7qur3zHofF/66jnbtnnurGFTwO45OqhZVl4cNg7qkNlia2CftTZw5tmBs+IOmaF2z8sDz3V0YKcD3wMAAK80V0E9AABMncmnAACAeSdTDwDAeHWWYkhLQT0AAOO2BEG98hsAAFhwMvUAAIzaMgxpKVMPAAALTqYeAIBxW4JM/VwF9V8+ek5+56nvn9r51v76zG3bPPftw76s+JYBEzIlycpLw9oNmQhq6KRSQ/9QT3thWMNe2X4s173PDZtUas8LA2epOjpspq3uJfhXCQBM1xKED8pvAABgwc1Vph4AAKap2oOyAADAAhDUAwAwbl3TXwaoqqur6uGqOlRVN22xv6rq/ZP991fVlZv2r1TVf6yqj23Xl6AeAIBx6xks26iqlSQfSHJNksuTXF9Vl29qdk2SyybLgSQ3b9r/niQPDXmLgnoAAJi+q5Ic6u5HuvtIkjuS7N/UZn+S23vdp5KcW1UXJElVXZTkHyb5jSGdCeoBABi1Yw/LTnMZ4MIkj21YPzzZNrTNv07yL5IMGkdcUA8AACfvvKq6Z8NyYNP+rQrvN38c2LJNVf1okie7+96hF2NISwAAxm02Q1o+3d37TrD/cJKLN6xflOTxgW3+hyQ/VlXXJjkzybdU1W91908er7O5Cuq/fuSM/IdHv2vbdnv3Dpve9YyvbP9FxItnnD7oXENmgE2SGjjz7NrK9m16z7Anq7f8jLfV+WpYwz0DJoFdeWngv46Xh80Um9WBvzgAgMVwd5LLqurSJF9Icl2Sf7KpzcEkN1bVHUnelOSZ7n4iyc9OllTVW5P88xMF9MmcBfUAADBVuzT5VHcfraobk9yVZCXJbd39YFW9e7L/liR3Jrk2yaEkzyd51zfb38yD+qp6NMmzSVaTHN3mawoAAJiuXZpRtrvvzHrgvnHbLRted5IbtjnHJ5N8cru+dipT/4Pd/fQO9QUAAEtF+Q0AAOO2S5n6nbQTQ1p2kj+sqnu3GOoHAAA4RTuRqX9Ldz9eVa9N8omq+nx3/+mxnZNA/0CSnHbet+7A5QAAsEx240HZnTbzTH13Pz75+WSSj2Z9ytyN+2/t7n3dvW/lW86e9eUAAMDozDSor6qzq+rVx14n+ZEkD8yyTwAAWDazLr95XZKP1vqkR6cl+Z3u/viM+wQAgL+1BOU3Mw3qu/uRJG8YfMALe7Ly0DnbNnv+24ZN7/rqr23fZnXvsC8ranXYX8Oel4fN2jpkhtoeMOtskqydNnCm2IHvYeXr27c5/dlhM8XWkQHT0yZZ6yX41wYAMCOGtAQAYLx2aUbZnSaoBwBg3JYgqN+JceoBAIAZkqkHAGDcZOoBAIB5J1MPAMBoVZbjQVmZegAAWHAy9QAAjNsSZOoF9QAAjJdx6nfenpeTsw9v/1vfc2RY1dA5T2w/bevpzw2bjfW0F4b9NayeMahZTntpin9dw95C9gyb3DUrA67ttOeHzSiblwe2W10d1m5tCf5VAgCcpLkK6gEAYOqWICfoQVkAAFhwMvUAAIzbEmTqBfUAAIzaMjwoq/wGAAAWnEw9AADjJlMPAADMO5l6AADGq7MUmXpBPQAAo7YMD8rOVVB/2otrec3DL2zb7qyvDJu2deXF7WeUfelbh/0Kjrx62LStvTKs3ZCZZ/ecvTLoXHsGTsa6Z+Dkrmd8bfvf28pzRwadq48MnMZ2qD0Dfr9mnQUAlsxcBfUAADB1S5Dv86AsAAAsOJl6AABGbRlq6mXqAQBgwcnUAwAwbkuQqRfUAwAwXksyTr3yGwAAWHAy9QAAjFZNlrGbq6C+XjyS0z//2LbtTvvq+YPO9/J5r9q2zd5nh03wdOTsYX8OR08f1CxHz9r+fHv2DjvXnoHzO532wrDvnuro9u3qyMCZrI4ObLe2/YRX6+2W4PszAICTNFdBPQAATN0S5AQF9QAAjJpx6gEAgLknqAcAYNx6BssAVXV1VT1cVYeq6qYt9ldVvX+y//6qunKy/cyq+vOq+k9V9WBV/avt+hLUAwDAlFXVSpIPJLkmyeVJrq+qyzc1uybJZZPlQJKbJ9tfSvJD3f2GJFckubqq3nyi/gT1AACM2+5k6q9Kcqi7H+nuI0nuSLJ/U5v9SW7vdZ9Kcm5VXTBZ//qkzemT5YS9CuoBABivXn9QdtpLkvOq6p4Ny4FNPV+YZONY7Ycn2wa1qaqVqrovyZNJPtHdnz7R2zT6DQAAnLynu3vfCfZvNSnR5mz7cdt092qSK6rq3CQfrarv6+4HjteZTD0AAOO2O+U3h5NcvGH9oiSPn2yb7v5qkk8mufpEnc1Vpr6Prmb16S9t225lddjso3tfes32bb44bEbZXhnWbu2cYdPArp024PNUDZvFdm3vsM9me44O+72d9pUXtm1TX3120LnWBs4o2z3wMfI9A34nZp0FAHbf3Ukuq6pLk3whyXVJ/smmNgeT3FhVdyR5U5JnuvuJqjo/ycvd/dWqOivJ25P8HyfqbK6CegAAmLbdmHyqu49W1Y1J7kqykuS27n6wqt492X9LkjuTXJvkUJLnk7xrcvgFST40GUFnT5KPdPfHTtTfzIP6qro6ya9l/c38Rnf/4qz7BACA3dbdd2Y9cN+47ZYNrzvJDVscd3+SN55MXzOtqR84PicAAMzOLk0+tZNmnan/m/E5k2RSL7Q/yedm3C8AACTZnfKbnTbr0W+GjM8JAACcglln6rcdn3MyUP+BJDkzr5rx5QAAsFTmtFxm2madqR8y9uat3b2vu/ednjNmfDkAADA+sw7q/2Z8zqram/XxOQ/OuE8AAPhbHpQ9Nccbn3OWfQIAwDGV5XhQdubj1G81Pufx/N2//135xD3/14yvCACAaar6nd2+hKVnRlkAAMZtCTL1s66pBwAAZkymHgCAUasef6peUA8AwHjN6Wg106b8BgAAFpxMPQAAo7YMQ1rK1AMAwIKTqQcAYNyWIFMvqAcAYNSU3wAAAHNPph4AgHGTqQcAAOadTD0AAOPVauoBAIAFIFMPAMC4LUGmXlAPAMBoVZTfAAAAC0CmHgCAcevxp+pl6gEAYMHJ1AMAMGrLUFMvqAcAYLw6SzH6jfIbAABYcDL1AACMWq3t9hXMnkw9AAAsOJl6AADGbQlq6gX1AACM2jKMfqP8BgAAFpxMPQAA49UxoywAADD/ZOoBABg1NfUAAMDcE9QDADBuPYNlgKq6uqoerqpDVXXTFvurqt4/2X9/VV052X5xVf1JVT1UVQ9W1Xu260v5DQAAo1XZnfKbqlpJ8oEkP5zkcJK7q+pgd39uQ7Nrklw2Wd6U5ObJz6NJfqa7P1NVr05yb1V9YtOx30CmHgAApu+qJIe6+5HuPpLkjiT7N7XZn+T2XvepJOdW1QXd/UR3fyZJuvvZJA8lufBEnQnqAQAYr+7ZLNu7MMljG9YP55WB+bZtquqSJG9M8ukTdab8BgAATt55VXXPhvVbu/vWDeu1xTGbPw2csE1VnZPk95K8t7u/dqKLEdQDADBqM6qpf7q7951g/+EkF29YvyjJ40PbVNXpWQ/of7u7f3+7i1F+AwDAuO3O6Dd3J7msqi6tqr1JrktycFObg0neMRkF581JnunuJ6qqknwwyUPd/StDOpOpBwCAKevuo1V1Y5K7kqwkua27H6yqd0/235LkziTXJjmU5Pkk75oc/pYkP5Xks1V132Tbz3X3ncfrT1APAMCo7daMspMg/M5N227Z8LqT3LDFcX+Wrevtj2tm5TdV9fNV9YWqum+yXDurvgAAYJnNOlP/q939SzPuAwAAttZJ1nYpVb+DlN8AADBu44/pZz76zY1VdX9V3VZVr9mqQVUdqKp7quqep556asaXAwAA43NKQX1V/VFVPbDFsj/JzUm+O8kVSZ5I8stbnaO7b+3ufd297/zzzz+VywEAgFeonv4yb06p/Ka73z6kXVX9epKPnUpfAADA1mZWU19VF3T3E5PVH0/ywKz6AgCA4+o5TK1P2SwflP0/q+qKrD+a8GiSn55hXwAAsLRmFtR390/N6twAADDUPNbAT5shLQEAGK+OIS0BAID5J1MPAMBoVZJaggdlZeoBAGDBydQDADBua7t9AbMnqAcAYNSU3wAAAHNPph4AgPEypCUAALAIZOoBABixTpagpl5QDwDAqNX4Y3rlNwAAsOhk6gEAGLclKL+RqQcAgAUnUw8AwHh1Ukswo6xMPQAALDiZegAAxm0JauoF9QAAjNv4Y3rlNwAAsOhk6gEAGLVagvIbmXoAAFhwMvUAAIzbEmTqBfUAAIxXJzFOPQAAMO9k6gEAGK1Ke1AWAACYfzL1AACM2xJk6gX1AACM2xIE9cpvAABgwcnUAwAwXoa0BAAAFoGgHgCAUavuqS+D+q26uqoerqpDVXXTFvurqt4/2X9/VV25Yd9tVfVkVT0wpC9BPQAATFlVrST5QJJrkm/0cmcAAAjvSURBVFye5PqqunxTs2uSXDZZDiS5ecO+30xy9dD+BPUAAIxb9/SX7V2V5FB3P9LdR5LckWT/pjb7k9ze6z6V5NyqumD9kvtPk3x56FsU1AMAMGIzCOjXg/rzquqeDcuBTR1fmOSxDeuHJ9tOts0gRr8BAICT93R37zvB/tpi2+YU/5A2gwjqAQAYr85uTT51OMnFG9YvSvL4N9FmEOU3AAAwfXcnuayqLq2qvUmuS3JwU5uDSd4xGQXnzUme6e4nvpnOZOoBABi3XZh8qruPVtWNSe5KspLktu5+sKrePdl/S5I7k1yb5FCS55O869jxVfXhJG/Neu3+4STv6+4PHq8/QT0AAKM2dFz5aevuO7MeuG/cdsuG153khuMce/3J9HVK5TdV9RNV9WBVrVXVvk37fnYykP7DVfUPTqUfAADg+E41U/9Akv8+yb/duHEysP51Sf5eku9M8kdV9Xe7e/UU+wMAgJOzS5n6nXRKmfrufqi7H95i1/4kd3T3S939F1mvE7rqVPoCAAC2Nqua+guTfGrD+nEH0p8M1H8gSV7/+tfP6HIAAFhKnWRt/Jn6bYP6qvqjJN+xxa7/rbv/4HiHbbFty99md9+a5NYk2bdv3/h/4wAA7KBeivKbbYP67n77N3HeqQ2kDwAAnNisJp86mOS6qjqjqi5NclmSP59RXwAAcHzd01/mzKkOafnjk8Hwvz/J/1NVdyVJdz+Y5CNJPpfk40luMPINAADMxik9KNvdH03y0ePs+4Ukv3Aq5wcAgFM2h5n1aZtV+Q0AALBDZjWkJQAA7D5DWgIAwKLrpNd2+yJmTvkNAAAsOJl6AADGzYOyAADAvJOpBwBgvDwoCwAAI6D8BgAAmHcy9QAAjJtMPQAAMO9k6gEAGLFeiky9oB4AgPHqJGtmlAUAAOacTD0AAOO2BOU3MvUAALDgZOoBABg3mXoAAGDeydQDADBinayNP1MvqAcAYLw66TakJQAAMOdk6gEAGLclKL+RqQcAgAUnUw8AwLgtwZCWgnoAAMarO1nzoCwAADDnZOoBABi3JSi/kakHAIAFJ1MPAMCo9RLU1AvqAQAYsVZ+AwAAfHOq6uqqeriqDlXVTVvsr6p6/2T//VV15dBjNxPUAwAwXp31GWWnvWyjqlaSfCDJNUkuT3J9VV2+qdk1SS6bLAeS3HwSx34DQT0AAEzfVUkOdfcj3X0kyR1J9m9qsz/J7b3uU0nOraoLBh77DdTUAwAwbr0rD8pemOSxDeuHk7xpQJsLBx77DQT1AABw8s6rqns2rN/a3bduWK8tjtlct3O8NkOO/QaCegAARquT9IAa+G/C09297wT7Dye5eMP6RUkeH9hm74Bjv4GaegAAxqt7vfxm2sv27k5yWVVdWlV7k1yX5OCmNgeTvGMyCs6bkzzT3U8MPPYbyNQDAMCUdffRqroxyV1JVpLc1t0PVtW7J/tvSXJnkmuTHEryfJJ3nejYE/UnqAcAYNRmVH6zfb/dd2Y9cN+47ZYNrzvJDUOPPRHlNwAAsOBk6gEAGLfdGdJyR9V61n8+VNVTSf5y0+bzkjy9C5fDN3If5oP7sPvcg/ngPuw+92A+zMt9+Dvdff5uX8RWqurjWf89TdvT3X31DM77TZmroH4rVXXPNsMFsQPch/ngPuw+92A+uA+7zz2YD+4Dx6ipBwCABSeoBwCABbcIQf2t2zdhB7gP88F92H3uwXxwH3afezAf3AeSLEBNPQAAcGKLkKkHAABOYK6D+qq6uqoerqpDVXXTbl/Psqiq26rqyap6YMO2b6uqT1TVf5n8fM1uXuPYVdXFVfUnVfVQVT1YVe+ZbHcfdlBVnVlVf15V/2lyH/7VZLv7sMOqaqWq/mNVfWyy7h7ssKp6tKo+W1X3VdU9k23uww6qqnOr6ner6vOT/z98v3vAMXMb1FfVSpIPJLkmyeVJrq+qy3f3qpbGbybZPO7qTUn+uLsvS/LHk3Vm52iSn+nu/zrJm5PcMPn7dx921ktJfqi735DkiiRXV9Wb4z7shvckeWjDunuwO36wu6/YMISi+7Czfi3Jx7v7e5O8Iev/JtwDksxxUJ/kqiSHuvuR7j6S5I4k+3f5mpZCd/9pki9v2rw/yYcmrz+U5B/t6EUtme5+ors/M3n9bNb/w31h3Icd1eu+Plk9fbJ03IcdVVUXJfmHSX5jw2b3YD64Dzukqr4lyQ8k+WCSdPeR7v5q3AMm5jmovzDJYxvWD0+2sTte191PJOsBZ5LX7vL1LI2quiTJG5N8Ou7DjpuUfdyX5Mkkn+hu92Hn/esk/yLJxnne3YOd10n+sKruraoDk23uw875riRPJfl3k1K036iqs+MeMDHPQX1tsc1QPSyVqjonye8leW93f223r2cZdfdqd1+R5KIkV1XV9+32NS2TqvrRJE929727fS3kLd19ZdbLYm+oqh/Y7QtaMqcluTLJzd39xiTPRakNG8xzUH84ycUb1i9K8vguXQvJF6vqgiSZ/Hxyl69n9Krq9KwH9L/d3b8/2ew+7JLJ19yfzPrzJu7DznlLkh+rqkezXob5Q1X1W3EPdlx3Pz75+WSSj2a9TNZ92DmHkxyefFuYJL+b9SDfPSDJfAf1dye5rKouraq9Sa5LcnCXr2mZHUzyzsnrdyb5g128ltGrqsp63eRD3f0rG3a5Dzuoqs6vqnMnr89K8vYkn4/7sGO6+2e7+6LuviTr/x/4f7v7J+Me7KiqOruqXn3sdZIfSfJA3Icd091/neSxqvqeyaa3Jflc3AMm5nryqaq6Nuu1lCtJbuvuX9jlS1oKVfXhJG9Ncl6SLyZ5X5L/O8lHkrw+yV8l+Ynu3vwwLVNSVf9dkv+Q5LP52zrin8t6Xb37sEOq6r/J+oNnK1lPgnyku//3qvr2uA87rqremuSfd/ePugc7q6q+K+vZ+WS9DOR3uvsX3IedVVVXZP2B8b1JHknyrkz+2xT3YOnNdVAPAABsb57LbwAAgAEE9QAAsOAE9QAAsOAE9QAAsOAE9QAAsOAE9QAAsOAE9QAAsOAE9QAAsOD+fxc50wdEzLxtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1008x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sc2.run(20)\n",
    "plt.scalar_field(sc2.velocity[:, 0.5, :, 0])\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then switch off the inflow..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc2.boundary_handling.trigger_reinitialization_of_boundary_data(activate=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwQAAAHSCAYAAABSAwz5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df4xd5Z3n+c+nblXZxgZsYkMc2wx02j3bdGvaYSzCTHZ70k0yY5jROlkpEkhLUBTJQQtSImU17cn80elZ9QhF+bEbDYJxOp6AJhOEJslgsd4wNJtsNtKSYLppYkNoPDQNBRUb89O/q+re7/5Rx51LpXzPt7jn1r2c835JR3XvOc95znPPufdWPfV9nvN1RAgAAABAM40NuwEAAAAAhocOAQAAANBgdAgAAACABqNDAAAAADQYHQIAAACgwegQAAAAAA02PuwGdFu/fn1cccUVw24GAAAAluDxxx8/FhEbht2OxfyzP1gdr77Wrrzex588+1BE7Ki84iEYqQ7BFVdcoQMHDgy7GQAAAFgC23877Dacz6uvtfWzhy6vvN7WxmfXV17pkIxUhwAAAACoUkjqqDPsZow05hAAAAAADUaEAAAAADUWagcRgl6IEAAAAAANRoQAAAAAtTU/hyCG3YyRRocAAAAAtcak4t4YMgQAAAA0GBECAAAA1FYo1A6GDPVChAAAAABoMCIEAAAAqDUmFfdGhwAAAAC1FZLadAh6YsgQAAAA0GB0CAAAAFBrHUXlS4btHbafsX3Y9u5Fttv214vtT9q+uli/0vbPbP+V7UO2/6Rrny/afsn2E8VyQ7/nhyFDAAAAQMVstyTdKemjkqYkPWZ7X0Q81VXseklbi+WDku4qfp6V9IcRccL2hKSf2P6/IuLRYr+vRcSXq2orHQIAAADUVkjDuu3oNZIOR8RzkmT7Pkk7JXV3CHZKujciQtKjttfa3hgR05JOFGUmimVgL4IhQwAAAKi1zgCWhE2SXux6PlWsS5Wx3bL9hKSjkh6OiJ92lbu9GGK01/a6XHPOjw4BAAAAsHTrbR/oWnYt2O5F9ln4X/7zlomIdkRsk7RZ0jW2f7fYfpek90vaJmla0lfe8SsoMGQIAAAAtRWKQd129FhEbO+xfUrSlq7nmyW9vNQyEfGG7R9J2iHpYEQcObfN9jckPbj0pr8dEQIAAACgeo9J2mr7StuTkm6UtG9BmX2SPlncbehaSW9GxLTtDbbXSpLtVZI+IukXxfONXft/XNLBfhtKhAAAAAD1FVJ7CHOKI2LO9u2SHpLUkrQ3Ig7ZvrXYfrek/ZJukHRY0ilJnyp23yjpnuJORWOS7o+Ic5GAL9nepvmhRc9L+ky/baVDAAAAAAxAROzX/B/93evu7nockm5bZL8nJX3gPHXeXHEz6RAAAACgvkLpuwI1Fh0CAAAA1JjVXvRmPjiHScUAAABAgxEhAAAAQG2FpM5QEhW/exAhAAAAABqMCAEAAABqjTkEvdEhAAAAQG2F6BCUYcgQAAAA0GBECAAAAFBrnSBC0AsRAgAAAKDBiBAAAACgtphDUI4OAQAAAGorZLUZFNMTZwcAAABoMCIEAAAAqDUmFfdGhAAAAABoMCIEAAAAqC0mFZejQwAAAIAas9rBoJheODsAAABAgxEhAAAAQG2FpA7/A++JswMAAAA0WN8dAttbbP/Q9tO2D9n+bLH+i7Zfsv1EsdzQf3MBAACApWnLlS91UsWQoTlJn4+Iv7B9oaTHbT9cbPtaRHy5gmMAAAAAGIC+OwQRMS1punh83PbTkjb1Wy8AAADQrwjuMlSm0rNj+wpJH5D002LV7baftL3X9rrz7LPL9gHbB1555ZUqmwMAAACoI1e+1EllHQLbayR9V9LnIuItSXdJer+kbZqPIHxlsf0iYk9EbI+I7Rs2bKiqOQAAAAASKrntqO0JzXcGvh0R35OkiDjStf0bkh6s4lgAAABA1nymYoYM9VLFXYYs6ZuSno6Ir3at39hV7OOSDvZ7LAAAAADVqiJC8CFJN0v6ue0ninVfkHST7W2a75g9L+kzFRwLAAAAWAImFZep4i5DP5EWnVmxv9+6AQAAgH6QqbgcZwcAAABosEomFQMAAACjqh31uk1o1YgQAAAAAA1GhAAAAAC1FTK3HS1BhwAAAAC11uEuQz1xdgAAAIAGI0IAAACA2iJTcTnODgAAANBgRAgAAABQWyFz29ESRAgAAACABiNCAAAAgFrr8D/wnugQAAAAoLYipDa3He2JswMAAAAMgO0dtp+xfdj27kW22/bXi+1P2r66WL/S9s9s/5XtQ7b/pGufS2w/bPvZ4ue6fttJhwAAAAA1ZnUGsJQe1W5JulPS9ZKuknST7asWFLte0tZi2SXprmL9WUl/GBG/J2mbpB22ry227Zb0SERslfRI8bwvdAgAAACA6l0j6XBEPBcRM5Luk7RzQZmdku6NeY9KWmt7Y/H8RFFmoliia597isf3SPpYvw1lDgEAAABqKzS0OQSbJL3Y9XxK0gcTZTZJmi4iDI9L+k1Jd0bET4syl0XEtCRFxLTtS/ttKB0CAAAA1NqAMhWvt32g6/meiNjT9XyxcUWx4Pl5y0REW9I222slfd/270bEwb5afB50CAAAAIClOxYR23tsn5K0pev5ZkkvL7VMRLxh+0eSdkg6KOlIMaxo2vZGSUffYfv/DnMIAAAAUFshqxPVLwmPSdpq+0rbk5JulLRvQZl9kj5Z3G3oWklvFn/obygiA7K9StJHJP2ia59bise3SHqgvzNEhAAAAACoXETM2b5d0kOSWpL2RsQh27cW2++WtF/SDZIOSzol6VPF7hsl3VPMIxiTdH9EPFhsu0PS/bY/LekFSZ/ot610CAAAAFBrA5pDUCoi9mv+j/7udXd3PQ5Jty2y35OSPnCeOl+VdF2V7aRDAAAAgNoKSR0yFffE2QEAAAAajAgBAAAAasxqJzILNxkRAgAAAKDBiBAAAACgtphDUI6zAwAAADQYEQIAAADUGnMIeqNDAAAAgNqKMEOGSnB2AAAAgAYjQgAAAIBaaxMh6ImzAwAAADQYEQIAAADUVkjqMKm4JzoEAAAAqDEzZKgEZwcAAABoMCIEAAAAqK35TMUMGeqFCAEAAADQYEQIAAAAUGtt/gfeEx0CAAAA1FbIDBkqQXcJAAAAaLC+OwS2t9j+oe2nbR+y/dli/SW2H7b9bPFzXf/NBQAAAJamo7HKlzqp4tXMSfp8RPy2pGsl3Wb7Kkm7JT0SEVslPVI8BwAAADBC+p5DEBHTkqaLx8dtPy1pk6Sdkj5cFLtH0o8k/VG/xwMAAACyIqQ2cwh6qjTeYfsKSR+Q9FNJlxWdhXOdhkvPs88u2wdsH3jllVeqbA4AAACAEpV1CGyvkfRdSZ+LiLey+0XEnojYHhHbN2zYUFVzAAAAAEnzicmqXuqkktuO2p7QfGfg2xHxvWL1EdsbI2La9kZJR6s4FgAAAJA1f9vRek0CrloVdxmypG9Kejoivtq1aZ+kW4rHt0h6oN9jAQAAAKhWFRGCD0m6WdLPbT9RrPuCpDsk3W/705JekPSJCo4FAAAALElb9RriU7Uq7jL0E+m8Z/m6fusHAAAAMDiVzCEAAAAARlFItZsEXDU6BAAAAKgxJhWX4ewAAAAADUaEAAAAALXWYVJxT0QIAAAAgAYjQgAAAIDaipDaTCruiQ4BAAAAao1Jxb1xdgAAAIAGI0IAAACA2gqZPAQliBAAAAAADUaEAAAAALXGbUd7I0IAAAAANBgRAgAAANRWSMwhKEGEAAAAALXWibHKlwzbO2w/Y/uw7d2LbLftrxfbn7R9dbF+i+0f2n7a9iHbn+3a54u2X7L9RLHc0O/5IUIAAAAAVMx2S9Kdkj4qaUrSY7b3RcRTXcWul7S1WD4o6a7i55ykz0fEX9i+UNLjth/u2vdrEfHlqtpKhwAAAAD1FUO77eg1kg5HxHOSZPs+STsldXcIdkq6NyJC0qO219reGBHTkqYlKSKO235a0qYF+1aGIUMAAABA9TZJerHr+VSxbkllbF8h6QOSftq1+vZiiNFe2+v6bSgdAgAAANRWaP62o1UvktbbPtC17Fpw6MXCErGUMrbXSPqupM9FxFvF6rskvV/SNs1HEb6y9LPydgwZAgAAQK0NaMjQsYjY3mP7lKQtXc83S3o5W8b2hOY7A9+OiO+dKxARR849tv0NSQ++o9Z3IUIAAAAAVO8xSVttX2l7UtKNkvYtKLNP0ieLuw1dK+nNiJi2bUnflPR0RHy1ewfbG7ueflzSwX4bSoQAAAAAtTWsPAQRMWf7dkkPSWpJ2hsRh2zfWmy/W9J+STdIOizplKRPFbt/SNLNkn5u+4li3RciYr+kL9nepvmX9rykz/TbVjoEAAAAwAAUf8DvX7Du7q7HIem2Rfb7iRafX6CIuLniZtIhAAAAQL2Rqbg3OgQAAACordDQ8hC8azCpGAAAAGgwIgQAAACotc7iw/FRIEIAAAAANBgRAgAAANRXMKm4DBECAAAAoMGIEAAAAKC2hpWY7N2EDgEAAABqjQ5BbwwZAgAAABqMCAEAAABqi8Rk5YgQAAAAAA1GhAAAAAC1FkQIeqJDAAAAgFojU3FvDBkCAAAAGmykIgR//cy0rvsn/7a0XOvUTKq+sRNnS8t01qzI1fXW6VS5WL0yVc5nZ8vrarVydUWkyoWXv3ecbZuy5TJVTeTOm6o+HxW+hmxd2dfaWTmRKzeZ+B9B8mWOzXZyBZOvtbOi/LXOXJz7SjuzNnfezl6ce4/MrS4vM7sm9zrbKyt8H0lqzeRew9jZ8nJjc7ljjpV/vUmSJk7kyq14o/y9tOLN3Ptt/GQ7Va51NlfO7fLjul3tNa30u2aUjwlUJMhUXIoIAQAAANBglXQIbO+1fdT2wa51X7T9ku0niuWGKo4FAAAALEWEK1/qpKohQ9+S9O8k3btg/dci4ssVHQMAAABYIvIQlKkkQhARP5b0WhV1AQAAAFg+g55DcLvtJ4shResGfCwAAADg1zBkqLdBdgjukvR+SdskTUv6ymKFbO+yfcD2gZnZkwNsDgAAAICFBnbb0Yg4cu6x7W9IevA85fZI2iNJF124ifuaAQAAoDIhbjtaZmARAtsbu55+XNLB85UFAAAAMByVRAhsf0fShyWttz0l6Y8lfdj2Ns13zJ6X9JkqjgUAAACkBbn1ylTSIYiImxZZ/c0l1zNmza0ub9Lki6/mKkxkkBxLvkN8PDe/IS7IZT72mfJsy3HhqlRdmku+y1sVhsuSWX7TH8BOsm2ZRLpjycDXeK5cZzKX1badyKQrSZ2J8tfqir+42ityr3VuVXm59mTumJ3xZJbf5Nv8zHvK6zu9OZdK97LLX0mV237JL1Plfmv1kdIyV644mqpr0/jrqXLtZID3qTObUuX+8sTlpWUOv7UhVdfUq2tT5U5PX5Aqd8FUdRm0nfy+zGbadi6hcc6w/mLhLyU0QEcMGeqFTMUAAABAgw1sUjEAAAAwbCHV7jahVSNCAAAAADQYEQIAAADUmLntaAk6BAAAAKg15s73xpAhAAAAoMGIEAAAAKDWmFTcGxECAAAAoMGIEAAAAKC2IogQlBmpDkG0rNk1iaDFmbO5+i65uLRM54Jc2tXWXC4d5dnLVqfKTSay6c5syKVwbZ3KZWdtr8xl0h3LZPLMfq4qnsQTY+UHzmbl7axIZtJdmatvdlUye3PiMrRyb3G1ZnInuDORq+/sReWvdfaiXF3ZDMTtlbnXMLd2trTMRZedSNX1m2uPpcptXPlmqlzG8XbuhJwcO1XZMSXpVCeXPf31mfKswa+fyr2G2dO5N1zrbO4zM5b4inMusXA+CzgzEIFa4S5DvTFkCAAAAGiwkYoQAAAAAFUj6NcbEQIAAACgwYgQAAAAoNaYVNwbHQIAAADUVsh0CEowZAgAAABoMDoEAAAAqLUYwJJhe4ftZ2wftr17ke22/fVi+5O2ry7Wb7H9Q9tP2z5k+7Nd+1xi+2HbzxY/172DU/I2dAgAAACAitluSbpT0vWSrpJ0k+2rFhS7XtLWYtkl6a5i/Zykz0fEb0u6VtJtXfvulvRIRGyV9EjxvC8jNYeg05JOX1LeR7l4Ipf05sRvre23SX9n5UQuqdfp9+RO6fjx8nJnL84dc2UmkZik9qpcfTpbnuGnM17tWLxMwjFJ6kwkEpNNVpcgTJI6rWpf61h5fi1NnsxlWWqdzpWbuyDX959ZU/5eSp+PZAao8dPJ5FSz5Z+Z48msaY+dzX2H/PWaDalyk63yxIVrV55O1fWbF76SKnfxeK6+vzn1nlS5F94q/wfTiVO5JGfRyX6ek99diSSCme8GKf+5V/J9Hu3ycq46QyOApRlepuJrJB2OiOckyfZ9knZKeqqrzE5J90ZESHrU9lrbGyNiWtK0JEXEcdtPS9pU7LtT0oeL/e+R9CNJf9RPQ4kQAAAAAEu33vaBrmXXgu2bJL3Y9XyqWLekMravkPQBST8tVl1WdBhU/Ly0nxchjViEAAAAAKjcYAJ1xyJie4/ti4UlFrakZxnbayR9V9LnIuKtpTcxhw4BAAAAam1IQ4amJG3per5Z0svZMrYnNN8Z+HZEfK+rzJFzw4psb5R0tN+GMmQIAAAAqN5jkrbavtL2pKQbJe1bUGafpE8Wdxu6VtKbxR/6lvRNSU9HxFcX2eeW4vEtkh7ot6FECAAAAFBrMYS5/RExZ/t2SQ9JaknaGxGHbN9abL9b0n5JN0g6LOmUpE8Vu39I0s2Sfm77iWLdFyJiv6Q7JN1v+9OSXpD0iX7bSocAAAAAGIDiD/j9C9bd3fU4JN22yH4/0eLzCxQRr0q6rsp20iEAAABAbYWGNofgXYMOAQAAAOorJNEh6IlJxQAAAECDjVSEIMals5eU9+DizNlUfW9tKX95q47lMr1OrMqdqmwWzNapmdIyY+2VqbrGTyVS3yqfybN1aq60jFckX2hyEk8kMx870cOPVq6f21l8aN6vGUvORMpmWx5rl9c3frI8860krfzlyVS5GM+dk7G51aVlOpO5z8LZ5PutvSKZrXZleTmvLf9cSdLm97yRKnfZquOpcpn30oXjue+trauOpMpdMZnLaLx+/ESq3Jl2efbm2XbuffRmO/f90BnPlYvEYTNlJCmczECcLMf/HIF3h2FMKn43IUIAAAAANNhIRQgAAACAyhEh6IkOAQAAAGrM3GWoBEOGAAAAgAYjQgAAAIB6Y8hQT0QIAAAAgAYjQgAAAID6CjIVlyFCAAAAADQYEQIAAADUG3MIehqpDkG0pNkLy69YnDmTqu/0e8vrWvVqqiq1V2bTYOaKeTaRiTZZ19iZ8szCkuQ15ZlIJcmdxDXIRt6S5dqTufPbXlleYXsymWE0l6RarZnchRibS5ZLXPqJ47ns02NvnUqVi1UrcuVa5edu5sJUVTq7IZdtOVbmyo1NlpfbcEkus/DvrJ1Olbt0Mlff0cRJ6WRT6VasnfwQnpqbLC1zZib3HdKeSWbGbmczlCcKZb97k+lKqy43FKPcNmDZMWSoF4YMAQAAAA1WSYfA9l7bR20f7Fp3ie2HbT9b/FxXxbEAAACAJYkBLDVSVYTgW5J2LFi3W9IjEbFV0iPFcwAAAAAjpJIOQUT8WNJrC1bvlHRP8fgeSR+r4lgAAADAkhAh6GmQk4ovi4hpSYqIaduXDvBYAAAAwK8LLeFuKM009EnFtnfZPmD7QPvkyWE3BwAAAGiUQXYIjtjeKEnFz6OLFYqIPRGxPSK2t1avHmBzAAAA0EQR1S91MsgOwT5JtxSPb5H0wACPBQAAAOAdqGQOge3vSPqwpPW2pyT9saQ7JN1v+9OSXpD0iSqOBQAAACxJzf6jX7VKOgQRcdN5Nl23pHostcuTZUqzuSyuM+vLM/jOrUiegotbqWJzq3KTVjoryjN+ZjPuxkSubZ1kNmAlkt9ms/xms3h22smsoJ3yc9JJZNuVJOVOm7IJZtOJaGcylSXrOp3L2u3ZXDbr8ZMXlZZpnc2dOM8m37+t3IlLJNDWiTO5jMy/PFP+OiXpdOoLSZo+XV7fXCf3OsczqawltZT7EB6bzaWWPjlb/lpnZ5Mfmrnca/VcMlNx4pRkv5PSn61sfQDeHZhU3NPQJxUDAAAAGJ5B3nYUAAAAGDozZKgnIgQAAABAgxEhAAAAQH3VMLNw1YgQAAAAAA1GhAAAAAA1Zu4yVIIOAQAAAOqNIUM9MWQIAAAAaDAiBAAAAKg3IgQ9jV6HYCxxxSbKs/ymJYeUpTPz5pKMpjL4tmaT2XtP5zI3j5/OnbfWqUR9TmahTZ7fsWR24c6K8nJjyazH2S+H1kyyYPa1Jq7r2JlcZuFIZu128no5ce6q/iykM8ImslQnE2NrPPkiVrUyaaWlC8Zz5TI2rXg9VW7ril+myl0wdjZV7oU160rLnJzJZW5+I/nBn01mKp47Xv6rqp1LUq3ORPK7azxZLvEazF8iAEbc6HUIAAAAgCrRL++JDgEAAADqK8RdhkowqRgAAABoMCIEAAAAqDUzZKgnIgQAAABAgxEhAAAAQL0RIeiJCAEAAADQYHQIAAAAgAGwvcP2M7YP2969yHbb/nqx/UnbV3dt22v7qO2DC/b5ou2XbD9RLDf02046BAAAAKg1R/VL6THtlqQ7JV0v6SpJN9m+akGx6yVtLZZdku7q2vYtSTvOU/3XImJbsexf0slYxEjNIXBHap1OZH2czGXcnXit/OWN5RK9avJ4Lu1qZ7yVKufZ8vrGkhly3cllXfVsMiVsoj63s+lqkxlyK8x+m60rK31ngmS5TKbizPtDkjSbzGg8nvuoZ85dOht3O5sdOZkRdoT/fTFW4e0rWsk3Urpcsm0TiYvfGkt+12Qyzkvp7N4pVY8PZrwxgP5dI+lwRDwnSbbvk7RT0lNdZXZKujciQtKjttfa3hgR0xHxY9tXLEdDR/hXLAAAAFCBcPWLtN72ga5l14KjbpL0YtfzqWLdUsss5vZiiNFe2+uWeDZ+DR0CAAAAYOmORcT2rmXPgu2LxUEXxh8zZRa6S9L7JW2TNC3pK6nW9jBSQ4YAAACASoWGNQxwStKWruebJb38Dsq8TUQcOffY9jckPdhfM4kQAAAAoO5iAEu5xyRttX2l7UlJN0rat6DMPkmfLO42dK2kNyNiuleltjd2Pf24pIPnK5tFhAAAAACoWETM2b5d0kOSWpL2RsQh27cW2++WtF/SDZIOSzol6VPn9rf9HUkf1vxchSlJfxwR35T0JdvbNN8teV7SZ/ptKx0CAAAA1FqFN4NbkuKWoPsXrLu763FIuu08+950nvU3V9lGiSFDAAAAQKMRIQAAAEC9kVukp5HqELgtTb6RSEw2kUtMtnqqvK5WMlnXxMlcAqiza3OJySKRsMuRe/d2VuTOR4znAkKdleX1dZJ1pWUTFFX4gc4muuqMZ5OrJRuXyW6YTDYXyfdI9vRm3nPpJHLZBHHZa5qI945lE2IldZJvkrlOdZ+HN9urUuVemsvddvpvz67P1Xfq4tIyb5y4IFXXzInJVLnW8dz35fipRJmzuWufTvjYTpZLfgYBDBkf1Z4YMgQAAAA02EhFCAAAAIAqOYY3qfjdgggBAAAA0GBECAAAAFBvkZ1N10x0CAAAAFBvDBnqiSFDAAAAQIMRIQAAAECtMam4NyIEAAAAQIMRIQAAAEC9ESHoaaQ6BGNtaeVriSu2IpcF8+LnZ/ts0a+Mv346VW7Fmtwp9dmZ0jITJ3LZkcdmcuXGT+YCQp5tp8plRCt3zLGZXFrbKt+w2QyjY7PZcrnX0DpZ/r70WydTdUU7d608njtznYny65XN8Jz97o3xZEbYVeWv9T2rEyltJf3G6mOpcpdNvJUq957JE6lyGf9g1YupclsmXk2VayfvrHHJiveVljkyeWGqrrPJbPLRSr5Lxrg7yNuQHRlAxUaqQwAAAABUisRkpQbeIbD9vKTjktqS5iJi+6CPCQAAAPwdOgQ9LVeE4A8iIhejBwAAALBsGDIEAACAeiNC0NNy3HY0JP1X24/b3rUMxwMAAACQtBwRgg9FxMu2L5X0sO1fRMSPz20sOgm7JGly9bplaA4AAACahEnFvQ08QhARLxc/j0r6vqRrFmzfExHbI2L7+KrVg24OAAAAgC4D7RDYXm37wnOPJf1TSQcHeUwAAAAAeYMeMnSZpO/bPnes/xQRPxjwMQEAAIBfYchQTwPtEETEc5J+L1vec9LK18uzvUYyU/Hk6+XZgNsrk5mFk9l7s1lts68hVddEq7K6JClcnhXUnWR22U4y63Hy/Hqu/LVm29Zp5bKfVj3uMFVfMsOzV6zIHXQsV1/rdPl1SGUTl9LZZWdP5j6DMxeXv4aXxtem6npicnOq3OWrX0+VG3f5eXvvZC7r8dpWLkv1hc5lYr+odSZ33MnybOyrV5R/p0rSyYnc+7JdYYw6870lSapD0uPsayWjMYAkbjsKAACA+iJTcSk6BAAAAKg3OgQ9LUceAgAAAAAjiggBAAAA6o0IQU9ECAAAAIAGI0IAAACA2rKYVFyGCAEAAADQYEQIAAAAUG9ECHqiQwAAAID6Ig9BqZHqELgdmnxzrrRcTE6k6ht/M5Ohc2WqrmjlsgGPzZZnWpakzgXlr6EznhvRNbcml/U4kvW5nfjUDGmwWSSy37Ync43L1CXlv0Q8l7z2k4lsy2vX5A564epUsdl1q1LlTm0szzA7mzukIps4tvwjL0lqnSmvcOat3Gfh+ZWXpMqdmMll3L1gojyD7+kLcm27oHU2Ve7F1qlUub8+/d5UuedPlJ+T42dy52NuNvl9mbz2iUTQ6QzljfovIRmNASSNVIcAAAAAqBz93p6YVAwAAAA0GBECAAAA1BsRgp7oEAAAAKDWmFTcG0OGAAAAgAYjQgAAAIB6I0LQExECAAAAoMGIEAAAAKC+QkQIShAhAAAAQK05ql9Sx7V32H7G9mHbuxfZbttfL7Y/afvqrm17bR+1fXDBPpfYftj2s8XPdf2en5GKELjd0cRb5Vk6x06eTtUXx0+Ulhk/e1GqLs/mUmpOJLMBK5HVdqLq3uyZ3GtIZdxNZuV1J1dO7eATEm0AABPgSURBVGS5TH1jyUzFE8m3fyuX7TObCVrJl5o65opcRtjORHV9fyfbH7mmaW5N7o0+c2n5+/eyTa+n6tq+4cVUud9e/XKq3KaJ8uNeMf5qqq7N47nP6Wwyu+wlrfLvQUl6a648m/WbZ3OZ3Y+/mcuMPTab+2xlkje3ZnLnI5tNPp35OFtuGMhADAyV7ZakOyV9VNKUpMds74uIp7qKXS9pa7F8UNJdxU9J+pakfyfp3gVV75b0SETcUXQydkv6o37aSoQAAAAA9RYDWMpdI+lwRDwXETOS7pO0c0GZnZLujXmPSlpre6MkRcSPJb22SL07Jd1TPL5H0sdSremBDgEAAABQvU2SukPSU8W6pZZZ6LKImJak4uelfbZztIYMAQAAAFUbUGKy9bYPdD3fExF7ug+7yD4LW5IpM3B0CAAAAIClOxYR23tsn5K0pev5ZkkLJ6dlyix0xPbGiJguhhcdzTb4fBgyBAAAgHobzhyCxyRttX2l7UlJN0rat6DMPkmfLO42dK2kN88NB+phn6Rbise3SHog1Zoe6BAAAACgvgbRGUh0CCJiTtLtkh6S9LSk+yPikO1bbd9aFNsv6TlJhyV9Q9L/cm5/29+R9P9J+vu2p2x/uth0h6SP2n5W83cwumOpp2QhhgwBAAAAAxAR+zX/R3/3uru7Hoek286z703nWf+qpOsqbCYdAgAAANSXtfjMXfzKSHUI3AmNnZopLzjXTlZYPiLKySRn2WRXPnkmVy6R6CzOVnt5nE3+ldFOXoPktYqzicxDkjQzW17XTOI9JCnmcgmgsuXSnPhaauWyeo2tWZ0qt3Ld2lS5ibcuLK/rohWpumYvyr2G0+ty5U6emigtc0S5ZI0Hk8m/Vo/n3pfvHX+ztMzFY+XvXUla31qTKjcbuc/WpvE3UuX+3qpjpWVeviCXyPHVVbn35dzkZKpcJ5EcMMaSv+6T5SJZnTP1tUkQBmC0jVSHAAAAAKgc/fKe6BAAAACg1gaUh6A2uMsQAAAA0GBECAAAAFBvRAh6IkIAAAAANBgRAgAAANQbEYKe6BAAAACgvoJJxWUYMgQAAAA0GBECAAAA1BsRgp5Gq0NgKxIZgWNdLpOnOuXlsu8Pd5JZfjNZaCV1VuQydFYp/VojUTJTRpKy2ZFX5bLfZrItp9ovpduWzXysZLblSGSpzmbG1njyI5x8/44dL8+0veJ0LuPu5LHcZ+GCydxruPhvyj8zM4dydZ1auzFV7gdr35cqt++if1xaZvbi3DVor85+1+Te557JvZfGzpSXayUTio+fyl37VeUJniVJK18rPyeTx3OZm8dmcufXnRr89ZD8fZSS/V4F8K40Wh0CAAAAoGLMIeht4HMIbO+w/Yztw7Z3D/p4AAAAAPIG2iGw3ZJ0p6TrJV0l6SbbVw3ymAAAAMDbxACWGhn0kKFrJB2OiOckyfZ9knZKemrAxwUAAAAkMWSozKCHDG2S9GLX86liHQAAAIARMOgIwWK3OHhbH832Lkm7JGnlxEUDbg4AAAAapYZDfKo26AjBlKQtXc83S3q5u0BE7ImI7RGxfXJ89YCbAwAAAKDboDsEj0naavtK25OSbpS0b8DHBAAAAH6FScU9DXTIUETM2b5d0kOSWpL2RsShQR4TAAAAOMdiUnGZgScmi4j9kvZnym696n166MD/NuAWAQAAoEr2nwy7CegDmYoBAABQb0QIehp4pmIAAAAAo4sIAQAAAGrNQYigFzoEAAAAqK8a3hWoagwZAgAAABqMCAEAAABqjduO9kaEAAAAAGgwIgQAAACoNyIEPdEhAAAAQK0xZKg3hgwBAAAADUaEAAAAAPVGhKAnIgQAAABAgxEhAAAAQH0FcwjKECEAAAAAGowIAQAAAOqNCEFPdAgAAABQWxZDhsowZAgAAABoMCIEAAAAqLcgRNALEQIAAACgwegQAAAAoNYc1S+p49o7bD9j+7Dt3Ytst+2vF9uftH112b62v2j7JdtPFMsN/Z4fhgwBAACgvkJDucuQ7ZakOyV9VNKUpMds74uIp7qKXS9pa7F8UNJdkj6Y2PdrEfHlqtpKhAAAAACo3jWSDkfEcxExI+k+STsXlNkp6d6Y96iktbY3JvetDB0CAAAA1Jo71S8JmyS92PV8qliXKVO27+3FEKO9ttclT8N50SEAAAAAlm697QNdy64F273IPgsHL52vTK9975L0fknbJE1L+soS2rwo5hAAAACg3gYzh+BYRGzvsX1K0pau55slvZwsM3m+fSPiyLmVtr8h6cElt3wBIgQAAACotSHdZegxSVttX2l7UtKNkvYtKLNP0ieLuw1dK+nNiJjutW8xx+Ccj0s62NfJERECAAAAoHIRMWf7dkkPSWpJ2hsRh2zfWmy/W9J+STdIOizplKRP9dq3qPpLtrdpPu7xvKTP9NtWOgQAAACor9DQMhVHxH7N/9Hfve7ursch6bbsvsX6mytuJkOGAAAAgCYjQgAAAIBay2YWbioiBAAAAECDESEAAABAvREh6IkOAQAAAGrLYshQGYYMAQAAAA1GhAAAAAD1FTG0246+WxAhAAAAABqMCAEAAABqjTkEvdEhAAAAQL3RIeiJIUMAAABAgxEhAAAAQK0xZKi3gUUIbH/R9ku2nyiWGwZ1LAAAAADvzKAjBF+LiC8P+BgAAADA4kJShxBBLwwZAgAAQL3RH+hp0JOKb7f9pO29ttctVsD2LtsHbB945ZVXBtwcAAAAAN366hDY/nPbBxdZdkq6S9L7JW2TNC3pK4vVERF7ImJ7RGzfsGFDP80BAAAAfo2j+qVO+hoyFBEfyZSz/Q1JD/ZzLAAAAADVG9gcAtsbI2K6ePpxSQcHdSwAAADgvKJm/9Kv2CAnFX/J9jbNT+N4XtJnBngsAAAAAO/AwDoEEXHzoOoGAAAAsuo25r9q3HYUAAAA9RXitqMlBn3bUQAAAAAjjAgBAAAAasuSzKTinogQAAAAAA1GhAAAAAD11hl2A0YbHQIAAADUGkOGemPIEAAAANBgRAgAAABQX9x2tBQRAgAAAKDBiBAAAACgxkJiDkFPdAgAAABQa6Y/0BNDhgAAAIAGI0IAAACAemPIUE9ECAAAAIAGI0IAAACA+grJZCruiQgBAAAA0GBECAAAAFBvzCHoiQ4BAAAA6o3+QE8MGQIAAAAajAgBAAAAas0MGeqJCAEAAAAwALZ32H7G9mHbuxfZbttfL7Y/afvqsn1tX2L7YdvPFj/X9dtOOgQAAACot4jqlxK2W5LulHS9pKsk3WT7qgXFrpe0tVh2Sborse9uSY9ExFZJjxTP+0KHAAAAAPUVkjoDWMpdI+lwRDwXETOS7pO0c0GZnZLujXmPSlpre2PJvjsl3VM8vkfSx7Kn4nzoEAAAAADV2yTpxa7nU8W6TJle+14WEdOSVPy8tN+GMqkYAAAAtWXFoCYVr7d9oOv5nojY87ZD/7qFDTlfmcy+laFDAAAAACzdsYjY3mP7lKQtXc83S3o5WWayx75HbG+MiOlieNHRd9L4bgwZAgAAQL0NYVKxpMckbbV9pe1JSTdK2regzD5JnyzuNnStpDeLYUC99t0n6Zbi8S2SHujv5BAhAAAAQN0NIQ9BRMzZvl3SQ5JakvZGxCHbtxbb75a0X9INkg5LOiXpU732Laq+Q9L9tj8t6QVJn+i3rXQIAAAAgAGIiP2a/6O/e93dXY9D0m3ZfYv1r0q6rsp20iEAAABAfZ277SjOizkEAAAAQIMRIQAAAECtDei2o7VBhAAAAABoMCIEAAAAqDciBD3RIQAAAECNpfMGNBZDhgAAAIAGI0IAAACA+goRIShBhAAAAABoMCIEAAAAqDcSk/VEhwAAAAC1Rh6C3voaMmT7E7YP2e7Y3r5g27+yfdj2M7b/WX/NBAAAADAI/UYIDkr6nyT9++6Vtq+SdKOk35H0Pkl/bvu3IqLd5/EAAACApSFC0FNfEYKIeDoinllk005J90XE2Yj4G0mHJV3Tz7EAAAAAVG9Qcwg2SXq06/lUse7X2N4laZckXX755QNqDgAAABopJHWIEPRS2iGw/eeS3rvIpn8dEQ+cb7dF1i16JSJij6Q9krR9+3auFgAAACpEpuIypR2CiPjIO6h3StKWruebJb38DuoBAAAAMECDSky2T9KNtlfYvlLSVkk/G9CxAAAAgPOLqH6pkX5vO/px21OS/pGk/9P2Q5IUEYck3S/pKUk/kHQbdxgCAAAARk9fk4oj4vuSvn+ebX8q6U/7qR8AAADoW83+o1+1QQ0ZAgAAAPAuMKjbjgIAAADDx21HS9EhAAAAQI2FFJ1hN2KkMWQIAAAAaDAiBAAAAKg3JhX3RIQAAAAAaDAiBAAAAKgvJhWXokMAAACAemPIUE8MGQIAAAAajAgBAAAA6o0IQU9ECAAAAIAGI0IAAACAGgsiBCXoEAAAAKC+QlKHTMW9MGQIAAAAaDAiBAAAAKg3hgz1RIQAAAAAaDAiBAAAAKg3IgQ9ESEAAAAAGowIAQAAAGospA4Rgl7oEAAAAKC+QorgtqO9MGQIAAAAaDA6BAAAAKi3TlS/9MH2JbYftv1s8XPdecrtsP2M7cO2d5ftb/sK26dtP1Esd2faQ4cAAAAAWF67JT0SEVslPVI8fxvbLUl3Srpe0lWSbrJ9VWL//xYR24rl1kxj6BAAAACg3iKqX/qzU9I9xeN7JH1skTLXSDocEc9FxIyk+4r9svun0SEAAABAfUVInU71S38ui4jp+ebFtKRLFymzSdKLXc+ninVl+19p+y9t/z+2/4dMY7jLEAAAALB0620f6Hq+JyL2nHti+88lvXeR/f51sn4vsq4sNDEt6fKIeNX2P5T0X2z/TkS81WsnOgQAAACot8FkKj4WEdvPf8j4yPm22T5ie2NETNveKOnoIsWmJG3per5Z0svF40X3j4izks4Wjx+3/d8k/Zak7o7Lr2HIEAAAALC89km6pXh8i6QHFinzmKSttq+0PSnpxmK/8+5ve0MxGVm2f0PSVknPlTWGCAEAAABqLfof81+1OyTdb/vTkl6Q9AlJsv0+SX8WETdExJzt2yU9JKklaW9EHOq1v6Tfl/RvbM9Jaku6NSJeK2sMHQIAAADUWCV3BapURLwq6bpF1r8s6Yau5/sl7V/C/t+V9N2ltochQwAAAECDESEAAABAfYX6zixcd0QIAAAAgAYjQgAAAIB6i5GbVDxSiBAAAAAADUaEAAAAALUVkoI5BD3RIQAAAEB9RTBkqARDhgAAAIAGI0IAAACAWmPIUG9ECAAAAIAGI0IAAACAemMOQU+OGJ0Qiu1XJP3tgtXrJR0bQnPwdlyH0cB1GD6uwWjgOgwf12A0jMp1+HsRsWHYjViM7R9o/jxV7VhE7BhAvctupDoEi7F9ICK2D7sdTcd1GA1ch+HjGowGrsPwcQ1GA9cBVWAOAQAAANBgdAgAAACABns3dAj2DLsBkMR1GBVch+HjGowGrsPwcQ1GA9cBfRv5OQQAAAAABufdECEAAAAAMCAj3SGwvcP2M7YP29497PY0he29to/aPti17hLbD9t+tvi5bphtrDvbW2z/0PbTtg/Z/myxnuuwjGyvtP0z239VXIc/KdZzHZaZ7Zbtv7T9YPGca7DMbD9v++e2n7B9oFjHdVhGttfa/s+2f1H8fvhHXANUYWQ7BLZbku6UdL2kqyTdZPuq4baqMb4laeF9dXdLeiQitkp6pHiOwZmT9PmI+G1J10q6rXj/cx2W11lJfxgRvydpm6Qdtq8V12EYPivp6a7nXIPh+IOI2NZ1m0uuw/L6PyT9ICL+O0m/p/nPBNcAfRvZDoGkayQdjojnImJG0n2Sdg65TY0QET+W9NqC1Tsl3VM8vkfSx5a1UQ0TEdMR8RfF4+Oa/9LfJK7Dsop5J4qnE8US4josK9ubJf1zSX/WtZprMBq4DsvE9kWSfl/SNyUpImYi4g1xDVCBUe4QbJL0YtfzqWIdhuOyiJiW5v9YlXTpkNvTGLavkPQBST8V12HZFUNVnpB0VNLDEcF1WH7/u6R/KanTtY5rsPxC0n+1/bjtXcU6rsPy+Q1Jr0j6D8XwuT+zvVpcA1RglDsEXmQdt0RCo9heI+m7kj4XEW8Nuz1NFBHtiNgmabOka2z/7rDb1CS2/4WkoxHx+LDbAn0oIq7W/FDe22z//rAb1DDjkq6WdFdEfEDSSTE8CBUZ5Q7BlKQtXc83S3p5SG2BdMT2Rkkqfh4dcntqz/aE5jsD346I7xWruQ5DUoTmf6T5+TVch+XzIUn/o+3nNT909A9t/0dxDZZdRLxc/Dwq6fuaH9rLdVg+U5KmiiilJP1nzXcQuAbo2yh3CB6TtNX2lbYnJd0oad+Q29Rk+yTdUjy+RdIDQ2xL7dm25seJPh0RX+3axHVYRrY32F5bPF4l6SOSfiGuw7KJiH8VEZsj4grN/x74vyPifxbXYFnZXm37wnOPJf1TSQfFdVg2EfFLSS/a/vvFquskPSWuASow0onJbN+g+bGjLUl7I+JPh9ykRrD9HUkflrRe0hFJfyzpv0i6X9Llkl6Q9ImIWDjxGBWx/d9L+n8l/Vy/Gjf9Bc3PI+A6LBPb/0Dzk/Ramv8Hyv0R8W9sv0dch2Vn+8OS/teI+Bdcg+Vl+zc0HxWQ5oeu/KeI+FOuw/KyvU3zk+snJT0n6VMqvpvENUAfRrpDAAAAAGCwRnnIEAAAAIABo0MAAAAANBgdAgAAAKDB6BAAAAAADUaHAAAAAGgwOgQAAABAg9EhAAAAABqMDgEAAADQYP8/aeN8RowG9mYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sc2.run(50)\n",
    "plt.scalar_field(sc2.velocity[:, 0.5, :, 0])\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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}