overview-ns-forw.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "o4JZ84moBKMr"
   },
   "source": [
    "# Navier-Stokes Forward Simulation\n",
    "\n",
    "Now let's target a somewhat more complex example: a fluid simulation based on the Navier-Stokes equations. This is still very simple with Φ<sub>Flow</sub> (phiflow), as differentiable operators for all steps exist there. The Navier-Stokes equations (in their incompressible form) introduce an additional pressure field $p$, and a constraint for conservation of mass, as introduced in equation {eq}`model-boussinesq2d`. We're also moving a marker field, denoted by $d$ here, with the flow. It indicates regions of higher temperature, and exerts a force via a buouyancy factor $\\xi$:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\frac{\\partial \\mathbf{u}}{\\partial{t}} + \\mathbf{u} \\cdot \\nabla \\mathbf{u} &= - \\frac{1}{\\rho} \\nabla p + \\nu \\nabla\\cdot \\nabla \\mathbf{u} + (0,1)^T \\xi d\n",
    "  \\quad \\text{s.t.} \\quad \\nabla \\cdot \\mathbf{u} = 0,\n",
    "  \\\\\n",
    "  \\frac{\\partial d}{\\partial{t}} + \\mathbf{u} \\cdot \\nabla d &= 0 \n",
    "\\end{aligned}$$\n",
    "\n",
    "\n",
    "Here, we're aiming for an incompressible flow (i.e., $\\rho = \\text{const}$), and use a simple buoyancy model (the Boussinesq approximation) via the term $(0,1)^T \\xi d$. This approximates changes in density for incompressible solvers, without explicitly calculating $\\rho$. We assume a gravity force that acts along the y direction via the vector $(0,1)^T$. \n",
    "We'll solve this PDE on a closed domain with Dirichlet boundary conditions $\\mathbf{u}=0$ for the velocity, and Neumann boundaries $\\frac{\\partial p}{\\partial x}=0$ for pressure, on a domain $\\Omega$ with a physical size of $100 \\times 80$ units. \n",
    "[[run in colab]](https://colab.research.google.com/github/tum-pbs/pbdl-book/blob/main/overview-ns-forw.ipynb)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "As in the previous section, the first command with a \"!\" prefix installs the [phiflow python package from GitHub](https://github.com/tum-pbs/PhiFlow) via `pip` in your python environment. (Skip or modify this command if necessary.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "da1uZcDXdVcF",
    "outputId": "1082dc87-796c-4b57-e72e-5790fc1444c9"
   },
   "outputs": [],
   "source": [
    "#!pip install --upgrade --quiet git+https://github.com/tum-pbs/PhiFlow@develop\n",
    "\n",
    "from phi.flow import *  # The Dash GUI is not supported on Google colab, ignore the warning\n",
    "import pylab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BVV1IKVqDfLl"
   },
   "source": [
    "## Setting up the simulation\n",
    "\n",
    "The following code sets up a few constants, which are denoted by upper case names. We'll use $40 \\times 32$ cells to discretize our domain, introduce a slight viscosity via $\\nu$, and define the time step to be $\\Delta t=1.5$. \n",
    "\n",
    "We're creating a first `CenteredGrid` here, which is initialized by a `Sphere` geometry object. This will represent the inflow region `INFLOW` where hot smoke is generated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "WrA3IXDxv31P"
   },
   "outputs": [],
   "source": [
    "DT = 1.5\n",
    "NU = 0.01\n",
    "\n",
    "INFLOW = CenteredGrid(Sphere(center=tensor([30,15], channel(vector='x,y')), radius=10), extrapolation.BOUNDARY, x=32, y=40, bounds=Box(x=(0,80),y=(0,100))) * 0.2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ExA0Pi2sFVka"
   },
   "source": [
    "The inflow will be used to inject smoke into a second centered grid `smoke` that represents the marker field $d$ from above. Note that we've defined a `Box` of size $100 \\times 80$ above. This is the physical scale in terms of spatial units in our simulation, i.e., a velocity of magnitude $1$ will move the smoke density by 1 unit per 1 time unit, which may be larger or smaller than a cell in the discretized grid, depending on the settings for `x,y`. You could parametrize your simulation grid to directly resemble real-world units, or keep appropriate conversion factors in mind. \n",
    "\n",
    "The inflow sphere above is already using the \"world\" coordinates: it is located at $x=30$ along the first axis, and $y=15$ (within the $100 \\times 80$ domain box).\n",
    "\n",
    "Next, we create grids for the quantities we want to simulate. For this example, we require a velocity field and a smoke density field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "smoke = CenteredGrid(0, extrapolation.BOUNDARY, x=32, y=40, bounds=Box(x=(0,80),y=(0,100)))  # sampled at cell centers\n",
    "velocity = StaggeredGrid(0, extrapolation.ZERO, x=32, y=40, bounds=Box(x=(0,80),y=(0,100)))  # sampled in staggered form at face centers "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We sample the smoke field at the cell centers and the velocity in [staggered form](https://tum-pbs.github.io/PhiFlow/Staggered_Grids.html). The staggered grid internally contains 2 centered grids with different dimensions, and can be converted into centered grids (or simply numpy arrays) via the `unstack` function, as explained in the link above.\n",
    "\n",
    "Next we define the update step of the simulation, which calls the necessary functions to advance the state of our fluid system by `dt`. The next cell computes one such step, and plots the marker density after one simulation frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 282
    },
    "id": "WmGZdOwswOva",
    "outputId": "3ae4d68d-b586-4bbe-eca9-a223d7720949"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max. velocity and mean marker density: [\u001b[94m0.1558497\u001b[0m, \u001b[94m0.008125\u001b[0m]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x11c9fa2a0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAGdCAYAAACmQjDKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAcO0lEQVR4nO3df2xV9f3H8ddF6BVs762ltLd3/bECCiqWZQzrjcrUVkpNCAgm+GOxbAQDK2SATu3iL9yWMkycuiEuMYEtsbBhrEQTYAK2xK2w0dFUdDa069Ya2uJIem+52Atf+vn+gV69UqC3vZ/ee93zkZyk955zz31z/nh6PD331mGMMQIAxNSYeA8AAN9ExBUALCCuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwYGy8B/i6gYEBHT9+XGlpaXI4HPEeBwDCjDHq6+uT1+vVmDGXPjdNuLgeP35ceXl58R4DAC6qs7NTubm5l9wm4eKalpb2+U+OzxcASBRGkvlKpy4u4eL65aUA4gogEZkhXbLkF1oAYAFxBQALiCsAWEBcAcAC4goAFhBXALCAuAKABcQVACwgrgBgAXEFAAuIKwBYQFwBwALiCgAWEFcAsIC4AoAFxBUALCCuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsIK4AYAFxBQALiCsAWEBcAcAC4goAFhBXALCAuAKABcQVACwgrgBgAXEFAAuiiuvmzZtVVFQkl8sll8sln8+nXbt2hdfffvvtcjgcEcuKFStiPjQAJLqx0Wycm5urDRs26JprrpExRr///e+1YMECHTlyRDfccIMkafny5XruuefCr5kwYUJsJwaAJBBVXOfPnx/x+Je//KU2b96sgwcPhuM6YcIEeTye2E0IAElo2Ndcz507p+3btysYDMrn84Wff/3115WZmakZM2aoqqpKp0+fvuR+QqGQAoFAxAIAyS6qM1dJ+uCDD+Tz+dTf36/U1FTV1tbq+uuvlyQ98MADKigokNfrVXNzsx5//HG1tLTozTffvOj+qqurtX79+uH/CwAgATmMMSaaF5w5c0YdHR3y+/1644039Nprr6m+vj4c2K/av3+/SkpK1NraqilTpgy6v1AopFAoFH4cCASUl5en8yfVjqj+MQBgl5E0IL/fL5fLdckto47r15WWlmrKlCn63e9+d8G6YDCo1NRU7d69W2VlZUPaXyAQkNvtFnEFkHiGHtcR3+c6MDAQceb5VU1NTZKknJyckb4NACSVqK65VlVVqby8XPn5+err61NNTY3q6uq0Z88etbW1qaamRnfffbcmTpyo5uZmrV27VnPmzFFRUZGt+QEgIUUV1xMnTuihhx5SV1eX3G63ioqKtGfPHt11113q7OzU3r179eKLLyoYDCovL0+LFy/Wk08+aWt2AEhYI77mGmtccwWQuEbxmisA4ELEFQAsIK4AYAFxBQALiCsAWEBcAcAC4goAFhBXALCAuAKABcQVACwgrgBgAXEFAAuIKwBYQFwBwALiCgAWEFcAsIC4AoAFxBUALCCuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsIK4AYAFxBQALiCsAWEBcAcAC4goAFhBXALCAuAKABcQVACwgrgBgQVRx3bx5s4qKiuRyueRyueTz+bRr167w+v7+flVWVmrixIlKTU3V4sWL1dPTE/OhASDRRRXX3NxcbdiwQY2NjTp8+LDuvPNOLViwQB9++KEkae3atXr77be1Y8cO1dfX6/jx41q0aJGVwQEgkTmMMWYkO8jIyNDzzz+ve++9V5MmTVJNTY3uvfdeSdLHH3+s6667Tg0NDbr55puHtL9AICC3263z3XeMZDQAiDEjaUB+v18ul+uSWw77muu5c+e0fft2BYNB+Xw+NTY26uzZsyotLQ1vM336dOXn56uhoeGi+wmFQgoEAhELACS7qOP6wQcfKDU1VU6nUytWrFBtba2uv/56dXd3KyUlRenp6RHbZ2dnq7u7+6L7q66ultvtDi95eXlR/yMAINFEHddp06apqalJhw4d0sqVK1VRUaGPPvpo2ANUVVXJ7/eHl87OzmHvCwASxdhoX5CSkqKpU6dKkmbNmqW///3veumll7RkyRKdOXNGvb29EWevPT098ng8F92f0+mU0+mMfnIASGAjvs91YGBAoVBIs2bN0rhx47Rv377wupaWFnV0dMjn8430bQAgqUR15lpVVaXy8nLl5+err69PNTU1qqur0549e+R2u7Vs2TKtW7dOGRkZcrlcWr16tXw+35DvFACAb4qo4nrixAk99NBD6urqktvtVlFRkfbs2aO77rpLkvTrX/9aY8aM0eLFixUKhVRWVqZXXnnFyuAAkMhGfJ9rrHGfK4DENQr3uQIALo64AoAFxBUALCCuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsIK4AYAFxBQALiCsAWEBcAcAC4goAFhBXALCAuAKABcQVACwgrgBgAXEFAAuIKwBYQFwBwALiCgAWEFcAsIC4AoAFxBUALCCuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsiCqu1dXVmj17ttLS0pSVlaWFCxeqpaUlYpvbb79dDocjYlmxYkVMhwaARBdVXOvr61VZWamDBw/q3Xff1dmzZzV37lwFg8GI7ZYvX66urq7wsnHjxpgODQCJbmw0G+/evTvi8datW5WVlaXGxkbNmTMn/PyECRPk8XhiMyEAJKERXXP1+/2SpIyMjIjnX3/9dWVmZmrGjBmqqqrS6dOnR/I2AJB0ojpz/aqBgQGtWbNGt9xyi2bMmBF+/oEHHlBBQYG8Xq+am5v1+OOPq6WlRW+++eag+wmFQgqFQuHHgUBguCMBQMJwGGPMcF64cuVK7dq1S++//75yc3Mvut3+/ftVUlKi1tZWTZky5YL1zz77rNavXz/IK8dIcgxnNACwxEgakN/vl8vluuSWw4rrqlWrtHPnTh04cECFhYWX3DYYDCo1NVW7d+9WWVnZBesHO3PNy8sTcQWQeIYe16guCxhjtHr1atXW1qquru6yYZWkpqYmSVJOTs6g651Op5xOZzRjAEDCiyqulZWVqqmp0c6dO5WWlqbu7m5Jktvt1vjx49XW1qaamhrdfffdmjhxopqbm7V27VrNmTNHRUVFVv4BAJCIoros4HAM/r/pW7Zs0dKlS9XZ2akf/OAHOnr0qILBoPLy8nTPPffoySefvOwp9BcCgYDcbre4LAAg8Vi+5moTcQWQuIYeV75bAAAsIK4AYAFxBQALiCsAWEBcAcAC4goAFhBXALCAuAKABcQVACwgrgBgAXEFAAuIKwBYQFwBwALiCgAWEFcAsIC4AoAFxBUALCCuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsIK4AYAFxBQALiCsAWEBcAcAC4goAFhBXALCAuAKABcQVACwgrgBgAXEFAAuIKwBYEFVcq6urNXv2bKWlpSkrK0sLFy5US0tLxDb9/f2qrKzUxIkTlZqaqsWLF6unpyemQwNAoosqrvX19aqsrNTBgwf17rvv6uzZs5o7d66CwWB4m7Vr1+rtt9/Wjh07VF9fr+PHj2vRokUxHxwAEpnDGGOG++JPP/1UWVlZqq+v15w5c+T3+zVp0iTV1NTo3nvvlSR9/PHHuu6669TQ0KCbb775svsMBAJyu906333HcEcDAAuMpAH5/X65XK5Lbjmia65+v1+SlJGRIUlqbGzU2bNnVVpaGt5m+vTpys/PV0NDw6D7CIVCCgQCEQsAJLthx3VgYEBr1qzRLbfcohkzZkiSuru7lZKSovT09Ihts7Oz1d3dPeh+qqur5Xa7w0teXt5wRwKAhDHsuFZWVuro0aPavn37iAaoqqqS3+8PL52dnSPaHwAkgrHDedGqVav0zjvv6MCBA8rNzQ0/7/F4dObMGfX29kacvfb09Mjj8Qy6L6fTKafTOZwxACBhRXXmaozRqlWrVFtbq/3796uwsDBi/axZszRu3Djt27cv/FxLS4s6Ojrk8/liMzEAJIGozlwrKytVU1OjnTt3Ki0tLXwd1e12a/z48XK73Vq2bJnWrVunjIwMuVwurV69Wj6fb0h3CgDAN0VUt2I5HIPfGrVlyxYtXbpU0vkPETzyyCPatm2bQqGQysrK9Morr1z0ssDXcSsWgMQ19FuxRnSfqw3EFUDiGqX7XAEAgyOuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsIK4AYAFxBQALiCsAWEBcAcAC4goAFhBXALCAuAKABcQVACwgrgBgAXEFAAuIKwBYQFwBwALiCgAWEFcAsIC4AoAFxBUALCCuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsIK4AYAFxBQALoo7rgQMHNH/+fHm9XjkcDr311lsR65cuXSqHwxGxzJs3L1bzAkBSiDquwWBQM2fO1KZNmy66zbx589TV1RVetm3bNqIhASDZjI32BeXl5SovL7/kNk6nUx6PZ9hDAUCys3LNta6uTllZWZo2bZpWrlypkydPXnTbUCikQCAQsQBAsot5XOfNm6c//OEP2rdvn371q1+pvr5e5eXlOnfu3KDbV1dXy+12h5e8vLxYjwQAo85hjDHDfrHDodraWi1cuPCi2/zrX//SlClTtHfvXpWUlFywPhQKKRQKhR8HAoHPAztGkmO4owGABUbSgPx+v1wu1yW3tH4r1uTJk5WZmanW1tZB1zudTrlcrogFAJKd9bh+8sknOnnypHJycmy/FQAkjKjvFjh16lTEWWh7e7uampqUkZGhjIwMrV+/XosXL5bH41FbW5see+wxTZ06VWVlZTEdHAASWdTXXOvq6nTHHXdc8HxFRYU2b96shQsX6siRI+rt7ZXX69XcuXP185//XNnZ2UPafyAQkNvtFtdcASSeoV9zHdEvtGwgrgASVwL9QgsA/hcRVwCwgLgCgAXEFQAsIK4AYAFxBQALiCsAWEBcAcAC4goAFhBXALCAuAKABcQVACwgrgBgAXEFAAuIKwBYQFwBwALiCgAWEFcAsIC4AoAFxBUALCCuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsIK4AYAFxBQALiCsAWEBcAcAC4goAFhBXALCAuAKABVHH9cCBA5o/f768Xq8cDofeeuutiPXGGD399NPKycnR+PHjVVpaqmPHjsVqXgBIClHHNRgMaubMmdq0adOg6zdu3KiXX35Zr776qg4dOqSrrrpKZWVl6u/vH/GwAJAsHMYYM+wXOxyqra3VwoULJZ0/a/V6vXrkkUf06KOPSpL8fr+ys7O1detW3XfffZfdZyAQkNvt1vnuO4Y7GgBYYCQNyO/3y+VyXXLLmF5zbW9vV3d3t0pLS8PPud1uFRcXq6GhYdDXhEIhBQKBiAUAkl1M49rd3S1Jys7Ojng+Ozs7vO7rqqur5Xa7w0teXl4sRwKAuIj73QJVVVXy+/3hpbOzM94jAcCIxTSuHo9HktTT0xPxfE9PT3jd1zmdTrlcrogFAJJdTONaWFgoj8ejffv2hZ8LBAI6dOiQfD5fLN8KABLa2GhfcOrUKbW2toYft7e3q6mpSRkZGcrPz9eaNWv0i1/8Qtdcc40KCwv11FNPyev1hu8oAID/BVHH9fDhw7rjjjvCj9etWydJqqio0NatW/XYY48pGAzq4YcfVm9vr2699Vbt3r1bV155ZeymRkyd/b+98R7hAuPGll5+IyCBjeg+Vxu4z3X0EVdgqOJ0nysA4DziCgAWEFcAsIC4AoAFxBUALCCuAGABcQUAC6L+EAGSSyLewzoUQ52b+2GRqDhzBQALiCsAWEBcAcAC4goAFhBXALCAuAKABcQVACwgrgBgAR8iSGLJ+gGBWBrKMeCDBogHzlwBwALiCgAWEFcAsIC4AoAFxBUALCCuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsIK4AYAFxBQALiCsAWEBcAcAC4goAFsQ8rs8++6wcDkfEMn369Fi/DQAkNCt/Q+uGG27Q3r1f/m2jsWP5U10A/rdYqd7YsWPl8Xhs7BoAkoKVa67Hjh2T1+vV5MmT9eCDD6qjo+Oi24ZCIQUCgYgFAJJdzONaXFysrVu3avfu3dq8ebPa29t12223qa+vb9Dtq6ur5Xa7w0teXl6sRwKAUecwxhibb9Db26uCggK98MILWrZs2QXrQ6GQQqFQ+HEgEPg8sGMkOWyOlvTO/t/ey28EjRtbGu8R8I1hJA3I7/fL5XJdckvrv2lKT0/Xtddeq9bW1kHXO51OOZ1O22MAwKiyfp/rqVOn1NbWppycHNtvBQAJI+ZxffTRR1VfX69///vf+utf/6p77rlHV1xxhe6///5YvxUAJKyYXxb45JNPdP/99+vkyZOaNGmSbr31Vh08eFCTJk2K9VsBQMKKeVy3b98e610CQNLhuwUAwALiCgAWEFcAsIC4AoAFxBUALCCuAGABcQUAC/gW6yQ2lC8k+aZ/uQtfyoJExZkrAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsIK4AYAFxBQAL+BDBN1yyftCADwcg2XHmCgAWEFcAsIC4AoAFxBUALCCuAGABcQUAC4grAFhAXAHAAuIKABbwCS3waSjAAs5cAcAC4goAFhBXALCAuAKABcQVACwgrgBgAXEFAAuIKwBYYC2umzZt0re//W1deeWVKi4u1t/+9jdbbwUACcdKXP/4xz9q3bp1euaZZ/SPf/xDM2fOVFlZmU6cOGHj7QAg4TiMMSbWOy0uLtbs2bP129/+VpI0MDCgvLw8rV69Wk888cQlXxsIBOR2u3W++45YjwYAI2AkDcjv98vlcl1yy5ifuZ45c0aNjY0qLf3y8+pjxoxRaWmpGhoaYv12AJCQYv7FLf/973917tw5ZWdnRzyfnZ2tjz/++ILtQ6GQQqFQ+HEgEIj1SAAw6uJ+t0B1dbXcbnd4ycvLi/dIADBiMY9rZmamrrjiCvX09EQ839PTI4/Hc8H2VVVV8vv94aWzszPWIwHAqIt5XFNSUjRr1izt27cv/NzAwID27dsnn893wfZOp1MulytiAYBkZ+XLstetW6eKigp973vf00033aQXX3xRwWBQP/zhDy/72i9vXoj5TQwAMELnuzSUm6ysxHXJkiX69NNP9fTTT6u7u1vf+c53tHv37gt+yTWYvr6+z38yIrAAElFfX9/nt4xenJX7XEdiYGBAx48fV1pamhyO8/e5BgIB5eXlqbOzM6kuGzD36EvW2Zl7dA13bmOM+vr65PV6NWbMpa+qJtzf0BozZoxyc3MHXZes12SZe/Ql6+zMPbqGM/flzli/EPdbsQDgm4i4AoAFSRFXp9OpZ555Rk6nM96jRIW5R1+yzs7co2s05k64X2gBwDdBUpy5AkCyIa4AYAFxBQALiCsAWJDwcU3Gv8X17LPPyuFwRCzTp0+P91gXOHDggObPny+v1yuHw6G33norYr0xRk8//bRycnI0fvx4lZaW6tixY/EZ9isuN/fSpUsvOP7z5s2Lz7BfUV1drdmzZystLU1ZWVlauHChWlpaIrbp7+9XZWWlJk6cqNTUVC1evPiCb5gbbUOZ+/bbb7/gmK9YsSJOE5+3efNmFRUVhT8o4PP5tGvXrvB628c6oeOazH+L64YbblBXV1d4ef/99+M90gWCwaBmzpypTZs2Dbp+48aNevnll/Xqq6/q0KFDuuqqq1RWVqb+/v5RnjTS5eaWpHnz5kUc/23bto3ihIOrr69XZWWlDh48qHfffVdnz57V3LlzFQwGw9usXbtWb7/9tnbs2KH6+nodP35cixYtiuPUQ5tbkpYvXx5xzDdu3Binic/Lzc3Vhg0b1NjYqMOHD+vOO+/UggUL9OGHH0oahWNtEthNN91kKisrw4/PnTtnvF6vqa6ujuNUl/fMM8+YmTNnxnuMqEgytbW14ccDAwPG4/GY559/Pvxcb2+vcTqdZtu2bXGYcHBfn9sYYyoqKsyCBQviMk80Tpw4YSSZ+vp6Y8z54ztu3DizY8eO8Db//Oc/jSTT0NAQrzEv8PW5jTHm+9//vvnJT34Sv6GG6OqrrzavvfbaqBzrhD1zTfa/xXXs2DF5vV5NnjxZDz74oDo6OuI9UlTa29vV3d0dcfzdbreKi4uT4vjX1dUpKytL06ZN08qVK3Xy5Ml4j3QBv98vScrIyJAkNTY26uzZsxHHfPr06crPz0+oY/71ub/w+uuvKzMzUzNmzFBVVZVOnz4dj/EGde7cOW3fvl3BYFA+n29UjnXCfXHLF6L9W1yJpLi4WFu3btW0adPU1dWl9evX67bbbtPRo0eVlpYW7/GGpLu7W5IGPf5frEtU8+bN06JFi1RYWKi2tjb97Gc/U3l5uRoaGnTFFVfEezxJ57/9bc2aNbrllls0Y8YMSeePeUpKitLT0yO2TaRjPtjckvTAAw+ooKBAXq9Xzc3Nevzxx9XS0qI333wzjtNKH3zwgXw+n/r7+5Wamqra2lpdf/31ampqsn6sEzauyay8vDz8c1FRkYqLi1VQUKA//elPWrZsWRwn+99w3333hX++8cYbVVRUpClTpqiurk4lJSVxnOxLlZWVOnr0aEJei7+Ui8398MMPh3++8cYblZOTo5KSErW1tWnKlCmjPWbYtGnT1NTUJL/frzfeeEMVFRWqr68flfdO2MsC0f4trkSWnp6ua6+9Vq2trfEeZci+OMbfhOM/efJkZWZmJszxX7Vqld555x299957EV+v6fF4dObMGfX29kZsnyjH/GJzD6a4uFiS4n7MU1JSNHXqVM2aNUvV1dWaOXOmXnrppVE51gkb12j/FlciO3XqlNra2pSTkxPvUYassLBQHo8n4vgHAgEdOnQo6Y7/J598opMnT8b9+BtjtGrVKtXW1mr//v0qLCyMWD9r1iyNGzcu4pi3tLSoo6Mjrsf8cnMPpqmpSZLifsy/bmBgQKFQaHSOdUx+LWbJ9u3bjdPpNFu3bjUfffSRefjhh016errp7u6O92iX9Mgjj5i6ujrT3t5u/vKXv5jS0lKTmZlpTpw4Ee/RIvT19ZkjR46YI0eOGEnmhRdeMEeOHDH/+c9/jDHGbNiwwaSnp5udO3ea5uZms2DBAlNYWGg+++yzhJ27r6/PPProo6ahocG0t7ebvXv3mu9+97vmmmuuMf39/XGde+XKlcbtdpu6ujrT1dUVXk6fPh3eZsWKFSY/P9/s37/fHD582Ph8PuPz+eI49eXnbm1tNc8995w5fPiwaW9vNzt37jSTJ082c+bMievcTzzxhKmvrzft7e2mubnZPPHEE8bhcJg///nPxhj7xzqh42qMMb/5zW9Mfn6+SUlJMTfddJM5ePBgvEe6rCVLlpicnByTkpJivvWtb5klS5aY1tbWeI91gffee++LP1QWsVRUVBhjzt+O9dRTT5ns7GzjdDpNSUmJaWlpie/Q5tJznz592sydO9dMmjTJjBs3zhQUFJjly5cnxH+QB5tZktmyZUt4m88++8z8+Mc/NldffbWZMGGCueeee0xXV1f8hjaXn7ujo8PMmTPHZGRkGKfTaaZOnWp++tOfGr/fH9e5f/SjH5mCggKTkpJiJk2aZEpKSsJhNcb+seYrBwHAgoS95goAyYy4AoAFxBUALCCuAGABcQUAC4grAFhAXAHAAuIKABYQVwCwgLgCgAXEFQAsIK4AYMH/A8Qnpm1BG/wWAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def step(velocity, smoke, pressure, dt=1.0, buoyancy_factor=1.0):\n",
    "    smoke = advect.semi_lagrangian(smoke, velocity, dt) + INFLOW\n",
    "    buoyancy_force = (smoke * (0, buoyancy_factor)).at(velocity)  # resamples smoke to velocity sample points\n",
    "    velocity = advect.semi_lagrangian(velocity, velocity, dt) + dt * buoyancy_force\n",
    "    velocity = diffuse.explicit(velocity, NU, dt)\n",
    "    velocity, pressure = fluid.make_incompressible(velocity)\n",
    "    return velocity, smoke, pressure\n",
    "\n",
    "velocity, smoke, pressure = step(velocity, smoke, None, dt=DT)\n",
    "\n",
    "print(\"Max. velocity and mean marker density: \" + format( [ math.max(velocity.values) , math.mean(smoke.values) ] ))\n",
    "\n",
    "pylab.imshow(np.asarray(smoke.values.numpy('y,x')), origin='lower', cmap='magma')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A lot has happened in this `step()` call: we've advected the smoke field, added an upwards force via a Boussinesq model, advected the velocity field, and finally made it divergence free via a pressure solve.\n",
    "\n",
    "The Boussinesq model uses a multiplication by a tuple `(0, buoyancy_factor)` to turn the smoke field into a staggered, 2 component force field, sampled at the locations of the velocity components via the `at()` function. This function makes sure the individual force components are correctly interpolated for the velocity components of the staggered velocity. Note that this also directly ensure the boundary conditions of the original grid are kept. It internally also does `StaggeredGrid(..., extrapolation.ZERO,...)` for the resulting force grid. \n",
    "\n",
    "The pressure projection step in `make_incompressible` is typically the computationally most expensive step in the sequence above. It solves a Poisson equation for the boundary conditions of the domain, and updates the velocity field with the gradient of the computed pressure.\n",
    "\n",
    "Just for testing, we've also printed the mean value of the velocities, and the max density after the update. As you can see in the resulting image, we have a first round region of smoke, with a slight upwards motion (which does not show here yet). \n",
    "\n",
    "## Datatypes and dimensions\n",
    "\n",
    "The variables we created for the fields of the simulation here are instances of the class `Grid`.\n",
    "Like tensors, grids also have the `shape` attribute which lists all batch, spatial and channel dimensions.\n",
    "[Shapes in phiflow](https://tum-pbs.github.io/PhiFlow/Math.html#shapes) store not only the sizes of the dimensions but also their names and types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Smoke: (xˢ=32, yˢ=40)\n",
      "Velocity: (xˢ=32, yˢ=40, vectorᶜ=x,y)\n",
      "Inflow: (xˢ=32, yˢ=40), spatial only: (xˢ=32, yˢ=40)\n"
     ]
    }
   ],
   "source": [
    "print(f\"Smoke: {smoke.shape}\")\n",
    "print(f\"Velocity: {velocity.shape}\")\n",
    "print(f\"Inflow: {INFLOW.shape}, spatial only: {INFLOW.shape.spatial}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the phiflow output here indicates the type of a dimension, e.g., $^S$ for a spatial, and $^V$ for a vector dimension. Later on for learning, we'll also introduce batch dimensions.\n",
    "\n",
    "The actual content of a shape object can be obtained via `.sizes`, or alternatively we can query the size of a specific dimension `dim` via `.get_size('dim')`. Here are two examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape content: (32, 40, 2)\n",
      "Vector dimension: 2\n"
     ]
    }
   ],
   "source": [
    "print(f\"Shape content: {velocity.shape.sizes}\")\n",
    "print(f\"Vector dimension: {velocity.shape.get_size('vector')}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The grid values can be accessed using the `values` property. This is an important difference to a phiflow tensor object, which does not have `values`, as illustrated in the code example below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Statistics of the different simulation grids:\n",
      "\u001b[92m(xˢ=32, yˢ=40)\u001b[0m \u001b[94m0.008 ± 0.039\u001b[0m \u001b[37m(0e+00...2e-01)\u001b[0m\n",
      "\u001b[92m(~vectorᵈ=x,y, xˢ=~\u001b[94m(x=31, y=32)\u001b[92m \u001b[93mint64\u001b[92m, yˢ=~\u001b[94m(x=40, y=39)\u001b[92m \u001b[93mint64\u001b[92m)\u001b[0m \u001b[94m-8.34e-09 ± 2.8e-02\u001b[0m \u001b[37m(-1e-01...2e-01)\u001b[0m\n",
      "Reordered test tensor shape: (2, 5, 3)\n"
     ]
    }
   ],
   "source": [
    "print(\"Statistics of the different simulation grids:\")\n",
    "print(smoke.values)\n",
    "print(velocity.values)\n",
    "\n",
    "# in contrast to a simple tensor:\n",
    "test_tensor = math.tensor(numpy.zeros([3, 5, 2]), spatial('x,y'), channel(vector=\"x,y\"))\n",
    "print(\"Reordered test tensor shape: \" + format(test_tensor.numpy('vector,y,x').shape) ) # reorder to vector,y,x \n",
    "#print(test_tensor.values.numpy('y,x')) # error! tensors don't return their content via \".values\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Grids have many more properties which are documented [here](https://tum-pbs.github.io/PhiFlow/phi/field/#phi.field.Grid).\n",
    "Also note that the staggered grid has a [non-uniform shape](https://tum-pbs.github.io/PhiFlow/Math.html#non-uniform-tensors) because the number of faces is not equal to the number of cells (in this example the x component has $31 \\times 40$ cells, while y has $32 \\times 39$). The `INFLOW` grid naturally has the same dimensions as the `smoke` grid.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time evolution\n",
    "\n",
    "With this setup, we can easily advance the simulation forward in time a bit more by repeatedly calling the `step` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "0hZk5HX3w4Or",
    "outputId": "f7811af7-4b58-4ff6-a8b6-6e7bedefaa6e"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed frame 0, max velocity 0.4630013406276703\n",
      "Computed frame 1, max velocity 0.8966456651687622\n",
      "Computed frame 2, max velocity 1.4098883867263794\n",
      "Computed frame 3, max velocity 2.041128158569336\n",
      "Computed frame 4, max velocity 2.927959680557251\n",
      "Computed frame 5, max velocity 3.8394827842712402\n",
      "Computed frame 6, max velocity 4.526943683624268\n",
      "Computed frame 7, max velocity 4.86798620223999\n",
      "Computed frame 8, max velocity 5.13107442855835\n",
      "Computed frame 9, max velocity 5.483879566192627\n"
     ]
    }
   ],
   "source": [
    "for time_step in range(10):\n",
    "    velocity, smoke, pressure = step(velocity, smoke, pressure, dt=DT)\n",
    "    print('Computed frame {}, max velocity {}'.format(time_step , np.asarray(math.max(velocity.values)) ))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GMKKWQBLHIwP"
   },
   "source": [
    "Now the hot plume is starting to rise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 282
    },
    "id": "Mfl80CjZxZcL",
    "outputId": "92f3a9ba-d403-4799-a543-132ee8ed234c"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x11ce5a840>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAGdCAYAAACmQjDKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAfz0lEQVR4nO3df3BU9f3v8deCZAXJbgg/8qMkNICCiqG3FONelaJEQvxeBgRn8EdHaBkcaGAKaKvpVRHaTijOWLVFbMf5is4YaPErcHW+QBVMuNpACyU3oiVDaCqxJMFSsxsWssTs5/6Bbl0JnF3YT3ZDn4+ZM5Pd895z3jk6rxzOfj7nuIwxRgCAhOqT7AYA4HJEuAKABYQrAFhAuAKABYQrAFhAuAKABYQrAFhAuAKABVcku4GvCofDOnbsmNLT0+VyuZLdDgBEGGPU3t6u3Nxc9elz4XPTlAvXY8eOKS8vL9ltAMB5NTU1afjw4ResSblwTU9P//wn1+cLAKQKI8l8KafOL+XC9V+XAghXAKnIxHTJki+0AMACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALIgrXNetW6fCwkJ5PB55PB75fD5t27Ytsn7y5MlyuVxRy8KFCxPeNACkuiviKR4+fLhWr16tq6++WsYYvfzyy5oxY4YOHDig66+/XpK0YMECrVq1KvKZAQMGJLZjAOgF4grX6dOnR73+2c9+pnXr1mnPnj2RcB0wYICys7MT1yEA9EIXfc21q6tLGzduVDAYlM/ni7z/6quvasiQIRo3bpzKy8t16tSpC24nFAopEAhELQDQ28V15ipJ77//vnw+nzo6OjRw4EBt3rxZ1113nSTpvvvu04gRI5Sbm6u6ujo98sgjqq+v1+uvv37e7VVUVGjlypUX/xsAQApyGWNMPB84c+aMjh49Kr/fr9dee00vvviiqqurIwH7Zbt27dKUKVPU0NCgUaNGdbu9UCikUCgUeR0IBJSXl6ezJ9WuuH4ZALDLSArL7/fL4/FcsDLucP2q4uJijRo1Sr/+9a/PWRcMBjVw4EBt375dJSUlMW0vEAjI6/WKcAWQemIP10se5xoOh6POPL+strZWkpSTk3OpuwGAXiWua67l5eUqLS1Vfn6+2tvbVVlZqaqqKu3YsUNHjhxRZWWl7rzzTg0ePFh1dXVatmyZJk2apMLCQlv9A0BKiitcjx8/rgceeEDNzc3yer0qLCzUjh07dMcdd6ipqUlvv/22nnnmGQWDQeXl5Wn27Nl67LHHbPUOACnrkq+5JhrXXAGkrh685goAOBfhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYAHhCgAWEK4AYEFc4bpu3ToVFhbK4/HI4/HI5/Np27ZtkfUdHR0qKyvT4MGDNXDgQM2ePVutra0JbxoAUl1c4Tp8+HCtXr1a+/fv1759+3T77bdrxowZ+uCDDyRJy5Yt0xtvvKFNmzapurpax44d06xZs6w0DgCpzGWMMZeygczMTD311FO6++67NXToUFVWVuruu++WJB06dEjXXnutampqdNNNN8W0vUAgIK/Xq7O577qU1gAgwYyksPx+vzwezwUrL/qaa1dXlzZu3KhgMCifz6f9+/ers7NTxcXFkZqxY8cqPz9fNTU1591OKBRSIBCIWgCgt4s7XN9//30NHDhQbrdbCxcu1ObNm3XdddeppaVFaWlpysjIiKrPyspSS0vLebdXUVEhr9cbWfLy8uL+JQAg1cQdrmPGjFFtba327t2rRYsWae7cufrwww8vuoHy8nL5/f7I0tTUdNHbAoBUcUW8H0hLS9Po0aMlSRMmTNCf/vQnPfvss5ozZ47OnDmjtra2qLPX1tZWZWdnn3d7brdbbrc7/s4BIIVd8jjXcDisUCikCRMmqF+/ftq5c2dkXX19vY4ePSqfz3epuwGAXiWuM9fy8nKVlpYqPz9f7e3tqqysVFVVlXbs2CGv16v58+dr+fLlyszMlMfj0ZIlS+Tz+WIeKQAAl4u4wvX48eN64IEH1NzcLK/Xq8LCQu3YsUN33HGHJOkXv/iF+vTpo9mzZysUCqmkpETPP/+8lcYBIJVd8jjXRGOcK4DU1QPjXAEA50e4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWBBXuFZUVGjixIlKT0/XsGHDNHPmTNXX10fVTJ48WS6XK2pZuHBhQpsGgFQXV7hWV1errKxMe/bs0VtvvaXOzk5NnTpVwWAwqm7BggVqbm6OLGvWrElo0wCQ6q6Ip3j79u1Rr9evX69hw4Zp//79mjRpUuT9AQMGKDs7OzEdAkAvdEnXXP1+vyQpMzMz6v1XX31VQ4YM0bhx41ReXq5Tp05dym4AoNeJ68z1y8LhsJYuXaqbb75Z48aNi7x/3333acSIEcrNzVVdXZ0eeeQR1dfX6/XXX+92O6FQSKFQKPI6EAhcbEsAkDJcxhhzMR9ctGiRtm3bpnfffVfDhw8/b92uXbs0ZcoUNTQ0aNSoUeesf/LJJ7Vy5cpuPtlHkutiWgMAS4yksPx+vzwezwUrLypcFy9erK1bt2r37t0qKCi4YG0wGNTAgQO1fft2lZSUnLO+uzPXvLw8Ea4AUk/s4RrXZQFjjJYsWaLNmzerqqrKMVglqba2VpKUk5PT7Xq32y232x1PGwCQ8uIK17KyMlVWVmrr1q1KT09XS0uLJMnr9ap///46cuSIKisrdeedd2rw4MGqq6vTsmXLNGnSJBUWFlr5BQAgFcV1WcDl6v6f6S+99JLmzZunpqYmfec739HBgwcVDAaVl5enu+66S4899pjjKfQXAoGAvF6vuCwAIPVYvuZqE+EKIHXFHq7cWwAALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcACwhUALCBcAcCCuMK1oqJCEydOVHp6uoYNG6aZM2eqvr4+qqajo0NlZWUaPHiwBg4cqNmzZ6u1tTWhTQNAqosrXKurq1VWVqY9e/borbfeUmdnp6ZOnapgMBipWbZsmd544w1t2rRJ1dXVOnbsmGbNmpXwxgEglbmMMeZiP/zJJ59o2LBhqq6u1qRJk+T3+zV06FBVVlbq7rvvliQdOnRI1157rWpqanTTTTc5bjMQCMjr9eps7rsutjUAsMBICsvv98vj8Vyw8pKuufr9fklSZmamJGn//v3q7OxUcXFxpGbs2LHKz89XTU1Nt9sIhUIKBAJRCwD0dhcdruFwWEuXLtXNN9+scePGSZJaWlqUlpamjIyMqNqsrCy1tLR0u52Kigp5vd7IkpeXd7EtAUDKuOhwLSsr08GDB7Vx48ZLaqC8vFx+vz+yNDU1XdL2ACAVXHExH1q8eLHefPNN7d69W8OHD4+8n52drTNnzqitrS3q7LW1tVXZ2dndbsvtdsvtdl9MGwCQsuI6czXGaPHixdq8ebN27dqlgoKCqPUTJkxQv379tHPnzsh79fX1Onr0qHw+X2I6BoBeIK4z17KyMlVWVmrr1q1KT0+PXEf1er3q37+/vF6v5s+fr+XLlyszM1Mej0dLliyRz+eLaaQAAFwu4hqK5XJ1PzTqpZde0rx58ySdnUTw0EMPacOGDQqFQiopKdHzzz9/3ssCX8VQLACpK/ahWJc0ztUGwhVA6uqhca4AgO4RrgBgAeEKABYQrgBgAeEKABYQrgBgAeEKABYQrgBgAeEKABYQrgBgAeEKABYQrgBgAeEKABYQrgBgAeEKABYQrgBgwUU9oBC9SV/HCpfLuUaS3P2GXGozkqRQ5z+ci0w4pm0ZdcVUBfQ0zlwBwALCFQAsIFwBwALCFQAsIFwBwALCFQAsIFwBwALCFQAsYBJBinLF8J9m5KD/cKyZO2S8Y83/+tqnMfWUN6wtpjonrSfSHWt2tQyOaVvr/97sWFMbqHSsCYdPxbQ/IFacuQKABYQrAFhAuAKABYQrAFhAuAKABYQrAFhAuAKABYQrAFjgMsak1G3aA4GAvF6vzua+K9ntWBDbXf/HDprtWPP/vu880eCKud92rDE5OTH1ZPoPiKnOieu084B91/HjMW3LvPGeY03xCq9jTXXgBed9mY6YesLlzEgKy+/3y+PxXLAy7jPX3bt3a/r06crNzZXL5dKWLVui1s+bN08ulytqmTZtWry7AYBeLe5wDQaDGj9+vNauXXvemmnTpqm5uTmybNiw4ZKaBIDeJu57C5SWlqq0tPSCNW63W9nZ2RfdFAD0dla+0KqqqtKwYcM0ZswYLVq0SCdOnDhvbSgUUiAQiFoAoLdLeLhOmzZNr7zyinbu3Kmf//znqq6uVmlpqbq6un8EckVFhbxeb2TJy8tLdEsA0OMSfsvBe+65J/LzDTfcoMLCQo0aNUpVVVWaMmXKOfXl5eVavnx55HUgECBgAfR61se5jhw5UkOGDFFDQ0O3691utzweT9QCAL2d9XD9+OOPdeLECeXEOJYSAC4HcV8WOHnyZNRZaGNjo2pra5WZmanMzEytXLlSs2fPVnZ2to4cOaIf/ehHGj16tEpKShLaeG/lveqamOrem5bmWNP3obucN9TU5FjS509/jqUlqfOz2Oqc9HP+384Mch74L0n6zp2OJds++i/HmoL//JZjTYv/3ZhaAqSLCNd9+/bptttui7z+4nrp3LlztW7dOtXV1enll19WW1ubcnNzNXXqVP3kJz+R2+1OXNcAkOLiDtfJkyfrQjNmd+zYcUkNAcDlgBu3AIAFhCsAWEC4AoAFhCsAWEC4AoAFhCsAWEC4AoAFCb9xCy4sv5/zTCBJyvjfRc5F//V7x5JTf3S+haO/9cpYWtKpUL+Y6pxc1f+MY403O7ZHqvSfcMixJm3BrY41o/+z3rGmRczQQuw4cwUACwhXALCAcAUACwhXALCAcAUACwhXALCAcAUACwhXALCASQQJ5IrhcN581ddj29j+vziWHN/mPND+4LGvOdacCcf2N9YVU5Wzrn86b+nKlu4fxf5V3/hnq2NNZp/3HWsmeIc51rwXcH70jiQZ4zxJApc/zlwBwALCFQAsIFwBwALCFQAsIFwBwALCFQAsIFwBwALCFQAsYBJBD+sMm5jq/L9vc6x596N8x5qTnzn//fT2i23Afixi+fVMDNMRmjtie+pB8K/OkyT+5/a/O9a87293rDEmcccJlz/OXAHAAsIVACwgXAHAAsIVACwgXAHAAsIVACwgXAHAAsIVACxgEkECGX3mWPNG8K2YtnX/X29xrNl7wnmgfXZ/51H97TFMNJCkK/s4b6srhkkEwS7nSQSnPovtuQc1nzj/L+zvdJ5ssa/ztRj2xiQCxC7uM9fdu3dr+vTpys3Nlcvl0pYtW6LWG2P0xBNPKCcnR/3791dxcbEOHz6cqH4BoFeIO1yDwaDGjx+vtWvXdrt+zZo1eu655/TCCy9o7969uuqqq1RSUqKODufnPQHA5SLuywKlpaUqLS3tdp0xRs8884wee+wxzZgxQ5L0yiuvKCsrS1u2bNE999xzad0CQC+R0C+0Ghsb1dLSouLi4sh7Xq9XRUVFqqmp6fYzoVBIgUAgagGA3i6h4drS0iJJysrKino/Kysrsu6rKioq5PV6I0teXl4iWwKApEj6UKzy8nL5/f7I0tTUlOyWAOCSJTRcs7OzJUmtra1R77e2tkbWfZXb7ZbH44laAKC3S2i4FhQUKDs7Wzt37oy8FwgEtHfvXvl8vkTuCgBSWtyjBU6ePKmGhobI68bGRtXW1iozM1P5+flaunSpfvrTn+rqq69WQUGBHn/8ceXm5mrmzJmJ7BsAUprLGBPbc0c+V1VVpdtuu+2c9+fOnav169fLGKMVK1boN7/5jdra2nTLLbfo+eef1zXXXBPT9gOBgLxer86eVMc2S6c36e8eHlPd7f3vc6xpM6cca76RPsix5o7szph6ykg741jTGXb+x9AnoTTHmv/zcWz/7Xd3/l/HmtOftTnWtJ380LEmlhl4uNwZSWH5/X7HS5hxn7lOnjxZF8pjl8ulVatWadWqVfFuGgAuG0kfLQAAlyPCFQAsIFwBwALCFQAsIFwBwALCFQAsIFwBwAIe89LDToc+jqlu+2cvONYMcGc51kztf69jzT/P9I2pJ8l58P/pLue/1/+IYX+D3TE1pGOfvBdboSMe4YLE4swVACwgXAHAAsIVACwgXAHAAsIVACwgXAHAAsIVACwgXAHAAiYRpKiurnbHmqHubzvWtJx2ftBE/76xTSLoF8Of4tNdzk8Q+PSMc82JUGyD+r8+qMSx5m+f/ndM2wISiTNXALCAcAUACwhXALCAcAUACwhXALCAcAUACwhXALCAcAUAC5hEkLKcB/+HY7h7/snOsGPNiVBskwiu7Ov8t7gjhkkEx2OY2NAZdu4bSGWcuQKABYQrAFhAuAKABYQrAFhAuAKABYQrAFhAuAKABYQrAFhAuAKABQmfofXkk09q5cqVUe+NGTNGhw4dSvSu/u0ZOc9iCsYwQyv4WWx/YwOdzrOvOmJ4Osupz5xnaBnnkrN1MRwDIBmsTH+9/vrr9fbbb/9rJ1cwyxbAvxcrqXfFFVcoOzvbxqYBoFewcs318OHDys3N1ciRI3X//ffr6NGj560NhUIKBAJRCwD0dgkP16KiIq1fv17bt2/XunXr1NjYqFtvvVXt7d0/KrqiokJerzey5OXlJbolAOhxLmNi/erg4rS1tWnEiBF6+umnNX/+/HPWh0IhhUKhyOtAIPB5wPaR5PwFyr+zEYOmOdZMcP0Px5rcq2K7OpTdPzFfaP096PwllP9MDBuStN8ccKz56NPtMW0LcGYkheX3++XxeC5Yaf2bpoyMDF1zzTVqaGjodr3b7Zbb7bbdBgD0KOvjXE+ePKkjR44oJyfH9q4AIGUkPFwffvhhVVdX629/+5v+8Ic/6K677lLfvn117733JnpXAJCyEn5Z4OOPP9a9996rEydOaOjQobrlllu0Z88eDR06NNG7Qgz+Gm51rOkfiu1fFf36OD8OpjPsfAk/luup9eb8I0yA3iDh4bpx48ZEbxIAeh3uLQAAFhCuAGAB4QoAFhCuAGAB4QoAFhCuAGAB4QoAFnAX614slhuSpA2a5VjTEsqIaX9dpn9MdU6auv7pWBPs82lM2+KmLEhVnLkCgAWEKwBYQLgCgAWEKwBYQLgCgAWEKwBYQLgCgAWEKwBYwCSCy9zhT193rEkbdH9M2zrR6fwgSSPnpwyc7hN0rPnbp/8dU09AquLMFQAsIFwBwALCFQAsIFwBwALCFQAsIFwBwALCFQAsIFwBwALCFQAsYIYW9MGnrya7BeCyw5krAFhAuAKABYQrAFhAuAKABYQrAFhAuAKABYQrAFhAuAKABdbCde3atfr617+uK6+8UkVFRfrjH/9oa1cAkHKshOtvf/tbLV++XCtWrNCf//xnjR8/XiUlJTp+/LiN3QFAynEZY0yiN1pUVKSJEyfqV7/6lSQpHA4rLy9PS5Ys0aOPPnrBzwYCAXm9Xp3NfVeiWwOAS2AkheX3++XxeC5YmfAz1zNnzmj//v0qLi7+10769FFxcbFqamoSvTsASEkJv3HLP/7xD3V1dSkrKyvq/aysLB06dOic+lAopFAoFHkdCAQS3RIA9LikjxaoqKiQ1+uNLHl5ecluCQAuWcLDdciQIerbt69aW1uj3m9tbVV2dvY59eXl5fL7/ZGlqakp0S0BQI9LeLimpaVpwoQJ2rlzZ+S9cDisnTt3yufznVPvdrvl8XiiFgDo7azcLHv58uWaO3euvvWtb+nGG2/UM888o2AwqO9+97uOn/3X4IWED2IAgEt0NpdiGWRlJVznzJmjTz75RE888YRaWlr0jW98Q9u3bz/nS67utLe3f/6TEQELIBW1t7d/PmT0/KyMc70U4XBYx44dU3p6ulyus+NcA4GA8vLy1NTU1KsuG9B3z+utvdN3z7rYvo0xam9vV25urvr0ufBV1ZR7hlafPn00fPjwbtf11muy9N3zemvv9N2zLqZvpzPWLyR9KBYAXI4IVwCwoFeEq9vt1ooVK+R2u5PdSlzou+f11t7pu2f1RN8p94UWAFwOesWZKwD0NoQrAFhAuAKABYQrAFiQ8uHaG5/F9eSTT8rlckUtY8eOTXZb59i9e7emT5+u3NxcuVwubdmyJWq9MUZPPPGEcnJy1L9/fxUXF+vw4cPJafZLnPqeN2/eOcd/2rRpyWn2SyoqKjRx4kSlp6dr2LBhmjlzpurr66NqOjo6VFZWpsGDB2vgwIGaPXv2OXeY62mx9D158uRzjvnChQuT1PFZ69atU2FhYWSigM/n07Zt2yLrbR/rlA7X3vwsruuvv17Nzc2R5d133012S+cIBoMaP3681q5d2+36NWvW6LnnntMLL7ygvXv36qqrrlJJSYk6Ojp6uNNoTn1L0rRp06KO/4YNG3qww+5VV1errKxMe/bs0VtvvaXOzk5NnTpVwWAwUrNs2TK98cYb2rRpk6qrq3Xs2DHNmjUriV3H1rckLViwIOqYr1mzJkkdnzV8+HCtXr1a+/fv1759+3T77bdrxowZ+uCDDyT1wLE2KezGG280ZWVlkdddXV0mNzfXVFRUJLErZytWrDDjx49PdhtxkWQ2b94ceR0Oh012drZ56qmnIu+1tbUZt9ttNmzYkIQOu/fVvo0xZu7cuWbGjBlJ6Scex48fN5JMdXW1Mebs8e3Xr5/ZtGlTpOYvf/mLkWRqamqS1eY5vtq3McZ8+9vfNj/4wQ+S11SMBg0aZF588cUeOdYpe+ba25/FdfjwYeXm5mrkyJG6//77dfTo0WS3FJfGxka1tLREHX+v16uioqJecfyrqqo0bNgwjRkzRosWLdKJEyeS3dI5/H6/JCkzM1OStH//fnV2dkYd87Fjxyo/Pz+ljvlX+/7Cq6++qiFDhmjcuHEqLy/XqVOnktFet7q6urRx40YFg0H5fL4eOdYpd+OWL8T7LK5UUlRUpPXr12vMmDFqbm7WypUrdeutt+rgwYNKT09PdnsxaWlpkaRuj/8X61LVtGnTNGvWLBUUFOjIkSP68Y9/rNLSUtXU1Khv377Jbk/S2bu/LV26VDfffLPGjRsn6ewxT0tLU0ZGRlRtKh3z7vqWpPvuu08jRoxQbm6u6urq9Mgjj6i+vl6vv/56EruV3n//ffl8PnV0dGjgwIHavHmzrrvuOtXW1lo/1ikbrr1ZaWlp5OfCwkIVFRVpxIgR+t3vfqf58+cnsbN/D/fcc0/k5xtuuEGFhYUaNWqUqqqqNGXKlCR29i9lZWU6ePBgSl6Lv5Dz9f3ggw9Gfr7hhhuUk5OjKVOm6MiRIxo1alRPtxkxZswY1dbWyu/367XXXtPcuXNVXV3dI/tO2csC8T6LK5VlZGTommuuUUNDQ7JbidkXx/hyOP4jR47UkCFDUub4L168WG+++abeeeedqNtrZmdn68yZM2pra4uqT5Vjfr6+u1NUVCRJST/maWlpGj16tCZMmKCKigqNHz9ezz77bI8c65QN13ifxZXKTp48qSNHjignJyfZrcSsoKBA2dnZUcc/EAho7969ve74f/zxxzpx4kTSj78xRosXL9bmzZu1a9cuFRQURK2fMGGC+vXrF3XM6+vrdfTo0aQec6e+u1NbWytJST/mXxUOhxUKhXrmWCfkazFLNm7caNxut1m/fr358MMPzYMPPmgyMjJMS0tLslu7oIceeshUVVWZxsZG895775ni4mIzZMgQc/z48WS3FqW9vd0cOHDAHDhwwEgyTz/9tDlw4ID56KOPjDHGrF692mRkZJitW7eauro6M2PGDFNQUGBOnz6dsn23t7ebhx9+2NTU1JjGxkbz9ttvm29+85vm6quvNh0dHUnte9GiRcbr9ZqqqirT3NwcWU6dOhWpWbhwocnPzze7du0y+/btMz6fz/h8viR27dx3Q0ODWbVqldm3b59pbGw0W7duNSNHjjSTJk1Kat+PPvqoqa6uNo2Njaaurs48+uijxuVymd///vfGGPvHOqXD1RhjfvnLX5r8/HyTlpZmbrzxRrNnz55kt+Rozpw5Jicnx6SlpZmvfe1rZs6cOaahoSHZbZ3jnXfe+eJBZVHL3LlzjTFnh2M9/vjjJisry7jdbjNlyhRTX1+f3KbNhfs+deqUmTp1qhk6dKjp16+fGTFihFmwYEFK/EHurmdJ5qWXXorUnD592nz/+983gwYNMgMGDDB33XWXaW5uTl7Txrnvo0ePmkmTJpnMzEzjdrvN6NGjzQ9/+EPj9/uT2vf3vvc9M2LECJOWlmaGDh1qpkyZEglWY+wfa245CAAWpOw1VwDozQhXALCAcAUACwhXALCAcAUACwhXALCAcAUACwhXALCAcAUACwhXALCAcAUACwhXALDg/wNfFI826xLOtQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pylab.imshow(smoke.values.numpy('y,x'), origin='lower', cmap='magma')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wnbQJvA-HPSL"
   },
   "source": [
    "Let's compute and show a few more steps of the simulation. Because of the inflow being located off-center to the left (with x position 30), the plume will curve towards the right when it hits the top wall of the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 489
    },
    "id": "tkhCOzc0ITsj",
    "outputId": "f6366c12-1eb5-4ff6-e0d7-94b806bfd8e4"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing time step 0\n",
      "Computing time step 1\n",
      "Computing time step 2\n",
      "Computing time step 10\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABQcAAAFNCAYAAAC9sdkvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQbklEQVR4nO3de3iU1bn//8/kNDknhEMCEhAFoYJgS5UGq7WKIm3duout2lqx22/bbdH9Veuusmu12u6itd+ttkXswao9UFqtVO3PQhUlHgoqKOIRBUFOCcdkcp4kM8/vD7bBlHWPCUxOs96v65rrgnXP/czKkPueZxbPzAoFQRAIAAAAAAAAgHfS+noCAAAAAAAAAPoGi4MAAAAAAACAp1gcBAAAAAAAADzF4iAAAAAAAADgKRYHAQAAAAAAAE+xOAgAAAAAAAB4isVBAAAAAAAAwFMsDgIAAAAAAACeYnEQAAAAAAAA8BSLgynoe9/7nkKhUF9PA0AfoQcAfqDWAX9R/4Df6AFINhYH8aEee+wxfe973+vy/e+66y7dd999SZ1DNBrVtddeqxEjRignJ0fTpk3T448/ntTHAODW1z1gxYoVCoVCztuqVauS9jiA7/q61hsaGnTjjTfqrLPOUklJiUKhUMLjv/nmmzrrrLOUn5+vkpISfeUrX9Hu3buTNh/AJwOp/i+55BLnOcGECROSNh/AN33dA1588UVdfvnlmjhxovLy8jRq1Ch98Ytf1Ntvv+28P+cAyZfR1xNA//fYY49pwYIFXW4Wd911l4YMGaJLLrkkaXO45JJL9OCDD+rKK6/UuHHjdN999+kzn/mMnnrqKX3yk59M2uMAOFh/6AGS9B//8R864YQTOo2NHTs2qY8B+Kyva33Pnj26+eabNWrUKE2ZMkUrVqww77tt2zadcsopKioq0g9/+EM1NDToxz/+sV599VW98MILysrKSsqcAF8MpPqXpHA4rF/96ledxoqKipIyF8BHfd0Dbr31Vj333HP6whe+oMmTJ6u6ulo/+9nP9LGPfUyrVq3SpEmTOu7LOUDPYHEQ/d4LL7ygxYsX67bbbtM111wjSbr44os1adIkffvb39Y//vGPPp4hgN5w8skn67zzzuvraQDoIcOHD1dVVZXKysq0evXqg/4z4IN++MMfqrGxUWvWrNGoUaMkSSeeeKLOOOMM3Xffffr617/eW9MGkATdqX9JysjI0EUXXdRLswPQ066++motWrSo08Le+eefr+OOO0633HKLfve733WMcw7QM/hY8QD37LPP6oQTTlB2draOPvpo/fznP+9y7jPPPKMvfOELGjVqlMLhsMrLy3XVVVepubm54z6XXHKJFixYIEmdLtu3HHnkkXr99ddVWVnZcd9TTz31kH8+SXrwwQeVnp7eqcizs7N16aWXauXKldq6dethHR8YyHzoAR9UX1+v9vb2pB0PGCh8qPVwOKyysrIu3ffPf/6zPve5z3W8KZCkGTNm6JhjjtGf/vSnw5oH0N9Q/26xWEx1dXWH9bjAQOBDD5g+ffpBV/yNGzdOEydO1JtvvtlpnHOAnsGVgwPYq6++qjPPPFNDhw7V9773PbW3t+vGG29UaWlpl/IfeOABNTU16bLLLtPgwYP1wgsv6Kc//am2bdumBx54QJL0jW98Qzt27NDjjz+u3/72tx96zDvuuENXXHGF8vPz9Z3vfEeSOuYTj8e1b9++Ls2tqKhImZmZkqSXX35ZxxxzjAoLCzvd58QTT5QkrV27VuXl5V06LpBKfOkB7/vqV7+qhoYGpaen6+STT9Ztt92mj3/84106HjCQ+VbrH2b79u3atWuXs/5PPPFEPfbYY906HtCfUf9uTU1NKiwsVFNTkwYNGqQLL7xQt956q/Lz8w/peEB/5XMPCIJAO3fu1MSJEzvGOAfoQQEGrHPPPTfIzs4O3nvvvY6xN954I0hPTw+68k/b1NR00Nj8+fODUCjU6Zhz587t0vHeN3HixOBTn/rUQeObNm0KJHXp9tRTT3U63mmnnXbQ8V5//fVAUnD33Xd3eW5AKvGlBzz33HPB7Nmzg3vuuSd4+OGHg/nz5weDBw8OsrOzg5deeqnL8wIGKl9q/YNefPHFQFJw7733mrHf/OY3B8X+8z//M5AUtLS0dPnnAPoz6v9g1113XXDttdcGf/zjH4M//OEPwZw5cwJJwUknnRS0tbV1+WcABgIfe8D7fvvb3waSgnvuuadjjHOAnsOVgwNULBbTsmXLdO6553a6nPYjH/mIZs6c2aUV85ycnI4/NzY2qrm5WdOnT1cQBHr55Zc7HTcZysrKurzD8JQpUzr+3NzcrHA4fNB9srOzO+KAb3zqAdOnT9f06dM7/v4v//IvOu+88zR58mTNmzdPS5cuTeo8gf7Ep1rvqvdf9z/s3MAVBwYS6t9t/vz5nf5+wQUX6JhjjtF3vvMdPfjgg7rgggsO6bhAf+NzD3jrrbc0d+5cVVRUaM6cOR3jnAP0HBYHB6jdu3erublZ48aNOyg2fvz4LjWKLVu26IYbbtAjjzyimpqaTrFIJJK0ub4vOztbM2bM6HZeTk6OotHoQeMtLS0dccA3PvUAl7Fjx+qcc87RQw89pFgspvT09KQcF+hvfK91l/df9zk3QKqj/rvuqquu0ne/+1098cQTLA4iZfjaA6qrq/XZz35WRUVFHfsPvI9zgJ7D4qCnYrGYzjjjDO3bt0/XXnutJkyYoLy8PG3fvl2XXHKJ4vF4jzzm7t27u3TfkpKSji8kHT58uLZv337QfaqqqiRJI0aMSN4kAU8MpB5gKS8vV2trqxobGw/6TlIA+6VCrf+z4cOHSzpwHvBBVVVVKikp4YoBQKlZ/5acnBwNHjy4y991BvhgIPaASCSiWbNmqba2Vs8888xB7/U5B+g5LA4OUEOHDlVOTo7eeeedg2Lr16//0PxXX31Vb7/9tu6//35dfPHFHeOuS4AT7VTkYt1/69atGjNmTJeO8dRTT3XseHT88cfrqaeeUl1dXacFgOeff74jDvjGpx5geffdd5Wdnc2XjyOlUesHO+KIIzR06FCtXr36oNgLL7zAeQFSBvXfdfX19dqzZ4+GDh2alOMB/YFvPaClpUVnn3223n77bT3xxBM69thjD8rhHKDnsDg4QKWnp2vmzJn6y1/+oi1btnR8V8Cbb76pZcuWdSlf2r8D0PuCINCdd9550H3z8vIkSbW1tSouLv7QY+fl5am2tvag8UP9/oHzzjtPP/7xj/WLX/xC11xzjaT9lxHfe++9mjZtGjsVw0s+9YDdu3cfdLL/yiuv6JFHHtGsWbOUlpbWpWMCA5FPtd4ds2fP1v3336+tW7d2nAcsX75cb7/9tq666qpDOibQ31D/B2tpaVFbW5sKCgo6jX//+99XEAQ666yzun1MoL/yqQfEYjGdf/75WrlypR5++GFVVFSYeZwD9AwWBwewm266SUuXLtXJJ5+sb37zm2pvb9dPf/pTTZw4UevWrUuYO2HCBB199NG65pprtH37dhUWFurPf/7zQd9DIElTp06VJP3Hf/yHZs6cqfT09ITf5TF16lQtXLhQP/jBDzR27FgNGzZMp5122iF//8C0adP0hS98QfPmzdOuXbs0duxY3X///dq8ebPuueeebh8PSBW+9IDzzz9fOTk5mj59uoYNG6Y33nhDv/jFL5Sbm6tbbrml28cDBhpfal2Sfvazn6m2tlY7duyQJD366KPatm2bJOmKK65QUVGRJOm//uu/9MADD+jTn/60/u///b9qaGjQbbfdpuOOO05f/epXD+mxgf6I+u9c/9XV1froRz+qCy+8UBMmTJAkLVu2TI899pjOOussnXPOOYf02EB/5UsP+Na3vqVHHnlEZ599tvbt26ff/e53neIXXXRRx585B+ghvb9BMpKpsrIymDp1apCVlRUcddRRwd133x3ceOONXdqG/I033ghmzJgR5OfnB0OGDAm+9rWvBa+88kogKbj33ns77tfe3h5cccUVwdChQ4NQKPShx66urg4++9nPBgUFBYEk5xbn3dXc3Bxcc801QVlZWRAOh4MTTjghWLp06WEfFxjofOgBd955Z3DiiScGJSUlQUZGRjB8+PDgoosuCt55553DOi4wkPhQ60EQBKNHjw4kOW+bNm3qdN/XXnstOPPMM4Pc3NyguLg4+PKXvxxUV1cf9hyA/ob6P1D/NTU1wUUXXRSMHTs2yM3NDcLhcDBx4sTghz/8YdDa2nrYcwD6Ix96wKc+9Smz/l1z4Rwg+UJB8IFrTAEAAAAAAAB4gy9qAgAAAAAAADzF4iAAAAAAAADgKRYHAQAAAAAAAE+xOAgAAAAAAAB4isVBAAAAAAAAwFMsDgIAAAAAAACeyujrCfyzeDyuHTt2qKCgQKFQqK+nA6SEIAhUX1+vESNGKC2tf/+fAD0ASL6B0gOofyD5Bkr9S/QAoCfQAwB/daf++93i4I4dO1ReXt7X0wBS0tatWzVy5Mi+nkZC9ACg5/T3HkD9Az2nv9e/RA8AehI9APBXV+q/3y0OFhQU/O+fQv97A3D4AknBB+qr/6IHAD1hYPQA6h/oCQOj/iV6ANAz6AGAv7pe//1ucfDA5cM0BCC5ggFxeT49AOgp/b8HUP9AT+n/9S/RA4CeQw8A/NW1+u/fXzoAAAAAAAAAoMewOAgAAAAAAAB4isVBAAAAAAAAwFMsDgIAAAAAAACeYnEQAAAAAAAA8BSLgwAAAAAAAICnWBwEAAAAAAAAPMXiIAAAAAAAAOApFgcBAAAAAAAAT7E4CAAAAAAAAHiKxUEAAAAAAADAUywOAgAAAAAAAJ5icRAAAAAAAADwFIuDAAAAAAAAgKdYHAQAAAAAAAA8xeIgAAAAAAAA4KmMvp4AAAAAAABAagqZkcyMwWbs6ILTzNjpeeOd48Nz7cdKt0OqbXWPr48YAUmvB2+bsW1NLzjHW6LVZk6gdjOGnseVgwAAAAAAAICnWBwEAAAAAAAAPMXiIAAAAAAAAOApFgcBAAAAAAAAT7E4CAAAAAAAAHiKxUEAAAAAAADAUxndufPChQu1cOFCbd68WZI0ceJE3XDDDZo1a5Yk6dRTT1VlZWWnnG984xu6++67kzNbAH2KHgD4i/oH/EYPAPxF/XdNdtYI5/jlI/6PmXP9qRvMWN4Ee7kmlLfDPZ4WMnOUKGYI2mJmLL43x4ytevyLzvE714fNnCU1d9qPFW8yY0iObi0Ojhw5UrfccovGjRunIAh0//3365xzztHLL7+siRMnSpK+9rWv6eabb+7Iyc3NTe6MAfQZegDgL+of8Bs9APAX9Q+kvm4tDp599tmd/v7f//3fWrhwoVatWtXRFHJzc1VWVpa8GQLoN+gBgL+of8Bv9ADAX9Q/kPoO+TsHY7GYFi9erMbGRlVUVHSM//73v9eQIUM0adIkzZs3T01NiS//jEajqqur63QD0P/RAwB/Uf+A3+gBgL+SVf8SPQDoT7p15aAkvfrqq6qoqFBLS4vy8/O1ZMkSHXvssZKkL33pSxo9erRGjBihdevW6dprr9X69ev10EMPmcebP3++brrppkP/CQD0KnoA4C/qH/AbPQDwV7LrX6IHAP1JKAiCoDsJra2t2rJliyKRiB588EH96le/UmVlZUdj+KAnn3xSp59+ujZs2KCjjz7aebxoNKpoNNrx97q6OpWXl2v/RY3d/8JMAC6BpLgikYgKCwsP60j0AGAgSk4PoP6BgYhzAMBv/fMcQEq9HtC7G5Jkusd7dUOSqBlb9Xipc5wNSXpb1+u/21cOZmVlaezYsZKkqVOn6sUXX9Sdd96pn//85wfdd9q0aZKUsCmEw2GFw/YvCID+hR4A+Iv6B/xGDwD8lez6l+gBQH/S7cXBfxaPxzut9n/Q2rVrJUnDhw8/3IcB0E/RAwB/Uf+A3+gBgL9Su/7tK+yGFHzUjN1x9FnO8S+cbV8dmD4oy57GIVzpp4wE20pkpduxdHcs0ZWI6SV5ZuykUe4r/U58u8rMeWDZVWbsPzdVOserI8+ZOfuvmkNXdWtxcN68eZo1a5ZGjRql+vp6LVq0SCtWrNCyZcu0ceNGLVq0SJ/5zGc0ePBgrVu3TldddZVOOeUUTZ48uafmD6AX0QMAf1H/gN/oAYC/qH8g9XVrcXDXrl26+OKLVVVVpaKiIk2ePFnLli3TGWecoa1bt+qJJ57QHXfcocbGRpWXl2v27Nm6/vrre2ruAHoZPQDwF/UP+I0eAPiL+gdSX7c3JOlpdXV1Kioq0kD9ElKgf0rel5H3NHoA0BMGRg+g/oGeMDDqX6IHAD2DHtB9yf5Y8WYz55A/Vpzu/vhwKNFHhw/hY8WH9NFmSYq2OYfb3o6YKQ8sO8qM8bHiQ9X1+k/wgXQAAAAAAAAAqYzFQQAAAAAAAMBTh71bMQAAAAAAQCrIzrJ3Wf7eqM+asS98dqNzPJRpX5MVr3N//FaSQgl2Hg7lGh8DTrRbcVqimPvjw6GsBEtGiT5ybORlHpdpplw4fIcZa77nVOf41Rvs3Y8bmt3/HnDjykEAAAAAAADAUywOAgAAAAAAAJ5icRAAAAAAAADwFIuDAAAAAAAAgKdYHAQAAAAAAAA8xeIgAAAAAAAA4KkE+1IDAAAAAACkopBz9MScL5gZ/+czG+zDxd3DbdXt3ZlUh4wi44CS0jLccw/lucclSen2tWGh7Ex3IB7Yx0skzZhfVoIlqCH5ZmjO5zY6x5/99YVmzm9bfmTGgqDVnoenuHIQAAAAAAAA8BSLgwAAAAAAAICnWBwEAAAAAAAAPMXiIAAAAAAAAOApFgcBAAAAAAAAT7E4CAAAAAAAAHgqwT7SAAAAAID+JnQIb+MCtffATICBKytzqHP8xmMzzZwgbh+veXvIPV6fZeZk5yaqy5gZycwxJpLmnoMkKdpmhoKY+3ihvLCdU99sxkI57rwgwRyUbl+7lnFMsXP8lulbzZxlT041YzsjK+15eIorBwEAAAAAAABPsTgIAAAAAAAAeIrFQQAAAAAAAMBTLA4CAAAAAAAAnmJxEAAAAAAAAPBUt7a5WrhwoRYuXKjNmzdLkiZOnKgbbrhBs2bNkiS1tLToW9/6lhYvXqxoNKqZM2fqrrvuUmlpadInDqD30QMAf1H/gN/oAYcu0c7CgwommbEvFp5txs4c7t7xM9GVH9tb7B1Yn98dOMf/3rLCzNlVv9aMxeMNCWaCgSZV639I3gTn+PFHVps5zbvs3YC3bh/kHG9stWuvrNCuldIiO5aW7+4roYIcM0fhBLsw76t3jretqbJzEuzcHDKaUdoge+fmtEHZ9gEz3AccPMne/bii8lNm7GG96Bz3eVf3bl05OHLkSN1yyy1as2aNVq9erdNOO03nnHOOXn/9dUnSVVddpUcffVQPPPCAKisrtWPHDn3+85/vkYkD6H30AMBf1D/gN3oA4C/qH0h9oSAI3P9N1EUlJSW67bbbdN5552no0KFatGiRzjvvPEnSW2+9pY985CNauXKlPvGJT3TpeHV1dSoqKtL+dUt7VR5AdwSS4opEIiosLEzqkekBwEDQMz2A+gcGAs4B+hpXDqJvDYxzAKn3e8CI4lOc46/OGmXmhDLsy+WSfuXgkXYsPNZ9lV3aiGIz51CuHGzfELFz+sGVg+3v2vM7f4H97/jwvh87x1PvysGu1/8hf+dgLBbT4sWL1djYqIqKCq1Zs0ZtbW2aMWNGx30mTJigUaNGaeXKlYf6MAD6KXoA4C/qH/AbPQDwF/UPpKZufeegJL366quqqKhQS0uL8vPztWTJEh177LFau3atsrKyVFxc3On+paWlqq62P7MfjUYVjUY7/l5XV9fdKQHoRfQAwF/UP+A3egDgr2TXv0QPAPqTbl85OH78eK1du1bPP/+8LrvsMs2ZM0dvvPHGIU9g/vz5Kioq6riVl5cf8rEA9Dx6AOAv6h/wGz0A8Fey61+iBwD9SbcXB7OysjR27FhNnTpV8+fP15QpU3TnnXeqrKxMra2tqq2t7XT/nTt3qqyszDzevHnzFIlEOm5bt27t9g8BoPfQAwB/Uf+A3+gBgL+SXf8SPQDoT7r9seJ/Fo/HFY1GNXXqVGVmZmr58uWaPXu2JGn9+vXasmWLKioqzPxwOKxwOHy40wDQR+gBgL+of8BvPvaAUMj+Mv2PFl3sHJ8/3r1JgSR98gR7MSSrbLM9j3zjbVxagk0c2u3dA77eEnOOxyL2lVzLn5puxm5e797c4PnI3WbO/i/Ox0BxuPUv9X0PmJnt3pAkPfyumbOrKt+MrdtX7Bzf0pRu5pyRZv/eF+yJmrGsEe6aDSJNZk4QtTfbaN/e7Bx/pvIIM2dfgo1WRue55zH2iL1mTtEke+5pg92blaQPtTcxOWGI3a8fjeQ6x2Mxfz/a3q3FwXnz5mnWrFkaNWqU6uvrtWjRIq1YsULLli1TUVGRLr30Ul199dUqKSlRYWGhrrjiClVUVHRrhyIA/Rc9APAX9Q/4jR4A+Iv6B1JftxYHd+3apYsvvlhVVVUqKirS5MmTtWzZMp1xxhmSpNtvv11paWmaPXu2otGoZs6cqbvuuqtHJg6g99EDAH9R/4Df6AGAv6h/IPWFgiDoV9dv19XVqaioSPu/DjHBZfEAuiGQFFckElFhYWFfTyYhegDQEwZGD6D+gZ4wMOpfGhg9oHc/Vmx/ZC/ZHysOzI8V2x9DXP6U/ZFjPlbcn9ADLF8tu945fvtZh/ax4ud3DnWOJ/xYcVnEjI0p22fGiie7x9OG5pg5/f9jxfa/ufWx4qDN3bsk6dafjTRjN2y8wzmeeh8r7nr9d3tDEgAAAAAAAACpgcVBAAAAAAAAwFMsDgIAAAAAAACe6taGJAAAAACQykIJ3iJ9qvCbZmzxDPf3gw2Zvtt+sDT7+8FCOfb3eSnL+A6ztEO79iPU7v7errRB9vd5nXneTjN24mr39xte98x3zJzf7vmFGWtt22XGgETS0nLN2MRB7u+4i0XtOtpWV2DGntvtzlvZ+J6ZMybP/k6/saPt3pE+dphzPGi1v1ewZrn7ewUl6edrxznHf737ZTMnS3b/Gq+jnOMXNdnfVXpGtv085WW1OsdDOXa/njrI/nnzs93Pe6Qx1b5zsOu4chAAAAAAAADwFIuDAAAAAAAAgKdYHAQAAAAAAAA8xeIgAAAAAAAA4CkWBwEAAAAAAABPsVsxAAAAAA+5dyqdUnyxmfGXc+1dc/MmGG+t3Bv37p9Bgp02lRe287KMvDT3z/Sh2t2TDDLazJS0DPs6k0GfiDrHby/YaOa0/vVrZuz3u293jsfjTWYOIEmZ6YVmbES2ezfutjb7d3tfq72L+ObGFuf4ptZVZk6k7QtmLFyUoHnkZbvHmyNmyq/XjTFjN2+62zk+JG+CmZMdyjdj1XH3PF6LDDVzPrK12IyNG17jHE/PNnZul3RUsf1cxOLu3Y99xpWDAAAAAAAAgKdYHAQAAAAAAAA8xeIgAAAAAAAA4CkWBwEAAAAAAABPsTgIAAAAAAAAeIrFQQAAAAAAAMBTGX09AQAAAADobRkZg5zj90zJM3Nyj4zaB2wPnMNBmntckkLxBLHsTPuxwll2zNIeM0NBS1v3j5dAKOx+m5kz3v55r91cY8aeWlPhHN9W+2SCWdiPBX8U544xY6XZrc7xhobsQ3qsI3LCzvGRoY+ZObVt9vVaoUSrNVH33GPv1Zkpd+3cZMba2t15RwXHmznH5ReZscy0kHN8aNiuy6Y2u+e17nGPZxfYvSstwROYneGee4OZkfq4chAAAAAAAADwFIuDAAAAAAAAgKdYHAQAAAAAAAA8xeIgAAAAAAAA4CkWBwEAAAAAAABPsTgIAAAAAAAAeCrR5tgHmT9/vh566CG99dZbysnJ0fTp03Xrrbdq/PjxHfc59dRTVVlZ2SnvG9/4hu6+++7kzBhAn6EHAP6i/gG/pWIPmJR/rnN84if2mDlB1D5evDHuHA9lhcycUJodU0ubnZeebucZgvpmOxiLuXNqEuRkJLjOJB44hxP9vOOm15qxb1VNd47Pa3rTzGlp3WHG0D0Duf4zQ7lmbFDYXdAtbfYySTjN/bstSSPz3DUxqGmwmdPYboaUbhxPklljT/5thJmyve7/M2NB4O43O9I2mTknZ33UjB1X5P7BijLtHzgnw+557VF374g3unuXJO1pzDFj6WlhI5KgJ8v+t08F3bpysLKyUnPnztWqVav0+OOPq62tTWeeeaYaGxs73e9rX/uaqqqqOm4/+tGPkjppAH2DHgD4i/oH/EYPAPxF/QOpr1tXDi5durTT3++77z4NGzZMa9as0SmnnNIxnpubq7KysuTMEEC/QQ8A/EX9A36jBwD+ov6B1HdY3zkYiUQkSSUlJZ3Gf//732vIkCGaNGmS5s2bp6amJvMY0WhUdXV1nW4ABgZ6AOAv6h/wGz0A8Fcy6l+iBwD9SbeuHPygeDyuK6+8UieddJImTZrUMf6lL31Jo0eP1ogRI7Ru3Tpde+21Wr9+vR566CHncebPn6+bbrrpUKcBoI/QAwB/Uf+A3+gBgL+SVf8SPQDoT0JBEBzStypedtll+tvf/qZnn31WI0eONO/35JNP6vTTT9eGDRt09NFHHxSPRqOKRg98EWhdXZ3Ky8u1/6LGRF8GCaDrAklxRSIRFRYWJuWI9ABgIEluD6D+gYGEcwDL8cX/5hxf9Q17E45Qgs9dBcb37CfakCR9UJb9WGUFdiwv256IoT9sSKJ2+61nLNJqxu7688G/P5I0b8M9Zg4bkryvf54DSL3TA0YWn27GHjn+WOd4Vrq9ycWmOrsuV9e4N7l4Yvc+M+fTQ0rM2Hcv2GDGMqeUOsf//v/s678+t/qXZqw9Vu8cP2rQZ8yc8wcld0OSUfkNZmzkyFrneG65maKXnx5qxs5Z94JzfGdklX3AAbkhSdfr/5CuHLz88sv117/+VU8//XTChiBJ06ZNkySzKYTDYYXD1k4xAPojegDgL+of8Bs9APBXMutfogcA/Um3FgeDINAVV1yhJUuWaMWKFRozZsyH5qxdu1aSNHz48EOaIID+gx4A+Iv6Pxz21Q8FuWPN2BeLzjdjXxnj/h6nSFummfP8XvsN2OKal53jmyPLzZx4PPF3SSG1DNwekG5Grh3j3jghaN9k5rQl+Eq0aL37sdKz4mZOfp59hVJamX1FkbKNKw5b7KvvtM99ZZAkNf4j4hxvb7X7V2a2fRVNpjH19AL77Wda2P63OmP4Xuf4bTvsHrqDKweTpv/Xv/27kxMqMmPZGe6r2Fpj9vEORU2a+/dXkoZlDzJjGWPtuQctbc7xv2zLN3PaY7VmzLKl/jkztiI+yowdU1DsHC/MdM9bksIJriq0Pu8aq7f763uNuWasuc26mnMgXh2YHN1aHJw7d64WLVqkhx9+WAUFBaqurpYkFRUVKScnRxs3btSiRYv0mc98RoMHD9a6det01VVX6ZRTTtHkyZN75AcA0HvoAYC/qH/Ab/QAwF/UP5D6urU4uHDhQknSqaee2mn83nvv1SWXXKKsrCw98cQTuuOOO9TY2Kjy8nLNnj1b119/fdImDKDv0AMAf1H/gN/oAYC/qH8g9XX7Y8WJlJeXq7Ky8rAmBKD/ogcA/qL+Ab/RAwB/Uf9A6kuwxRQAAAAAAACAVMbiIAAAAAAAAOCpbn2sGAAAwG/2DoYnFV3mHH/0czVmTvFX7Z02Fdg78KnG+IhX2P7o19kF9v8Jf3/Qac7xUPVxZs7tV9vPxQ2bfuccb2yxd4EFekJGeoEZO2lktXO8vcGuo13b7OPtbnDvjDkop8XMySm1tz9O3+PeQViSlOWuv6DKzmlY3WjGnnu93Dm+odHe5Xx0rr0z8pRhe5zjw8bZc8gosnvK2I+5++jgNaPNHPYq9kcoZP/uDI+7f7clKR64d85tabePV9WSacZe3ec+XlXb62bO+PyRZiw0rNiMxd9x9681DfZuwIeivd3uKeuij5mxx3d8yTk+aqz93KanJfj4ety9c3pbvb2j+jsN9r9VffN79mN5iisHAQAAAAAAAE+xOAgAAAAAAAB4isVBAAAAAAAAwFMsDgIAAAAAAACeYnEQAAAAAAAA8BSLgwAAAAAAAICnMvp6AgAAAP1Luhn5ZNFcM/b0khLneDD6ePuh1rxuhuI7au289sCOWRL8l3CoKNsdGF1q5ly5pMyMnfz1f3GOn/VypZmzr/4VMyYdws8LSMrPGWnGMjJizvGG3WEz54097jqXpA2NWc7x8pwcM2f4vjozlnX0EWZM7e65x17baaYsXn20GfvNllrn+NbQWjNnXDDJjF3U6u4PM7O2mTmDjmo1Y0HM3QOOD48yc14Luf89JCkI7MdCaonJXSuSFBgvLbWt9u/OM7vs16OljYuc48U5o82cI/IbzZjSiu2YcR6wM7TJzjkkcTPS0rbXjP2j/WXn+AmRqWbOiPyGrk/rf9XsyjVjq/dQ593BlYMAAAAAAACAp1gcBAAAAAAAADzF4iAAAAAAAADgKRYHAQAAAAAAAE+xOAgAAAAAAAB4isVBAAAAAAAAwFMZfT0BAACA/uTYQV80Y08/MMhO3FfnHA5eetZMCRJNJMv+P9xQQaY70B63H6stZsbie5vdgd2bzJy0d3aYsak/negcf+iydDNn5ovVZizaZseARI7MONGMZWS0Ocd37i0wc16JhM3Y6j2tzvFji7PMnLMy7S4QpCXoAc2NzvGG9fbxlmyNmrGV9b92jhfkjDZzNmbZbyU3N37cOb6r1n5uC+r2mrE0o+VlhEJmTihkzy8I3P9WGJiCwH59y0yw5BEL2p3jK/fmmDmPNDzQ9Yn9r2Eaa8YGFTTZie2FZiiU4/654rKfi0Nj95R43K6jmlb3+cMr+z5m5nx8kP2859W5H+uvm0eYOZWtD5uxIHD3f59x5SAAAAAAAADgKRYHAQAAAAAAAE+xOAgAAAAAAAB4isVBAAAAAAAAwFMsDgIAAAAAAACe6tZuxfPnz9dDDz2kt956Szk5OZo+fbpuvfVWjR8/vuM+LS0t+ta3vqXFixcrGo1q5syZuuuuu1RaWpr0yQPoXfQAwF+pWf/unXN/ObHMzAh22Ltptr+zzzkeSrd300zLN7bglBTKtHf2VYbx/7tpCXbutI+moMm9Y2PQau9+3N7o3i1VktKXr3WOn/zbT5o5x00924ytrv2VEUm43zOSqH/3APu3+wgNMWNtbbuc43tbss2cjXV2TbwQ/4dzPFZTYebU7bR3P86t3mPGFHHX36p3jjBT/l53lxkLAvfPNTzrODPnX4onmLHyXPfx4oH9b9XWkqBLGbH3jF2bJSlIsJMquqd/139iuSF7t/C6qPu19Ind7tdzSWpprTFj2VmDnOPFgb3rcOHQnWbMfK2X1LbZ/bvfGm+wj5dsxm7PktTW7t6FeXOzPb+XE+xm/g9jF/mFO1ebOTUNr5sxJX1X54GvW1cOVlZWau7cuVq1apUef/xxtbW16cwzz1TjB04Or7rqKj366KN64IEHVFlZqR07dujzn/980icOoPfRAwB/Uf+A3+gBgL+ofyD1devKwaVLl3b6+3333adhw4ZpzZo1OuWUUxSJRHTPPfdo0aJFOu200yRJ9957rz7ykY9o1apV+sQnPpG8mQPodfQAwF/UP+A3egDgL+ofSH2H9Z2DkUhEklRSUiJJWrNmjdra2jRjxoyO+0yYMEGjRo3SypUrnceIRqOqq6vrdAMwMNADAH9R/4Df6AGAv5JR/xI9AOhPDnlxMB6P68orr9RJJ52kSZMmSZKqq6uVlZWl4uLiTvctLS1VdXW18zjz589XUVFRx628vPxQpwSgF9EDAH9R/4Df6AGAv5JV/xI9AOhPDnlxcO7cuXrttde0ePHiw5rAvHnzFIlEOm5bt249rOMB6B30AMBf1D/gN3oA4K9k1b9EDwD6k2595+D7Lr/8cv31r3/V008/rZEjR3aMl5WVqbW1VbW1tZ3+12Dnzp0qK3Pv/BcOhxUO27t0Aeh/6AGAv6h/wG/0AMBfyax/iR4A9CfdWhwMgkBXXHGFlixZohUrVmjMmDGd4lOnTlVmZqaWL1+u2bNnS5LWr1+vLVu2qKKiInmzBtAn6AGAv1Kx/suK3POquK3UzIk++Ip9wHjgHE7LDpkpQVbMjIWK7DdMoZIC9/Hqmsyc+N5mO9bsnkesPm7Pwf6xFLREneOZT68xc24+ZpAZO3uN++eNxfh+qt4yUHvAFlWZsba2LOd4U7v9Fikjzf7F39v4tnP8lbxsM6ep+UgzprwcMxTscf/uP7PHfqwg3mI/lqEusD8SGtIEMzY61/1YOZlt9vES9JS6fe7nYmPam2ZOoHb7gOiW/l7/Idm/PJHAfl2sb891ju9K22zmxBLUUSze6hzPVLqZk1WW4Bd/iP26WLNxu3M8kP26nWyB3Oc9ktTaHnGOv9y+zMyp3nGcGdsXe8893rjezFEvPhepoFuLg3PnztWiRYv08MMPq6CgoOP7A4qKipSTk6OioiJdeumluvrqq1VSUqLCwkJdccUVqqioYIciIAXQAwB/Uf+A3+gBgL+ofyD1dWtxcOHChZKkU089tdP4vffeq0suuUSSdPvttystLU2zZ89WNBrVzJkzdddddyVlsgD6Fj0A8Bf1D/iNHgD4i/oHUl+3P1b8YbKzs7VgwQItWLDgkCcFoH+iBwD+ov4Bv9EDAH9R/0DqO+TdigEAAAAAAAAMbCwOAgAAAAAAAJ5icRAAAAAAAADwVLe+cxAAcKhCCSLpZiwvZ7QZGxv+VLdnsaltpRmra9roHA+CtgRH/PDvoAH6q2N1gnO87ZFXzJz2Ovt4rQ3u06r0cNzMySuwjxcqyLGDsZgds8Ttem3c7B6vr8kzc0Ih+3jFw5qd48FrtWbOjM/Vm7FgdbsZAxK9FrWEGs1YtC23m0eThmYneD0Pua+72NX0mpmTkz3CjMWPOdqMpdU3uR+r2e43QcKfzN1TdjXYc1+qMWZsfMFI5/gR+fa/R3N9phm7f/0RzvGa6LNmTqJzL85fUkui3+3atJoEme4e0KYWMyMe2K9H0baIezzLzkkL2+8D4mWlZizW/oZzPDMtwblD0tnnIu0x92t6XZP93NY1bbIfynjegwRzoM67hysHAQAAAAAAAE+xOAgAAAAAAAB4isVBAAAAAAAAwFMsDgIAAAAAAACeYnEQAAAAAAAA8BS7FfcL9g5FoZAdC2cO6fYjRdv22MHAvbsZOwABXZeZ4a7Li4b8u5nzg09sM2Nl/5pgp9AJxk7GcXunQm3KN0P1jx/vHL/t6XFmzh07fm/GGprdux8D/cVppe6tgttqqsyc3dvtGqqqcx9vZLG9xXFWsb1zZ6jA3r037cihzvH4u3vNnPo37dfzDe+5j7dyT6GZU5Jl95qKtt3O8eHhBNs9JzjfyMoodo63tLp3bAXelxcUm7FY3H2dRKI9bjfWJdqp1H2eHQplmTmJdv0ONTTYE8lyv43LTk80++5rb3fvvipJ65seN2O/ee+LzvGGmP3+Zdcme+4/qfqjc7yuab2ZA5/Yrx974++ZsUib+1w6W/ZrvbUruSS1x9w1WxPaZ+YkulwrVG+fBwwa2ewcz1C2fcBe5f43CYw1hw/HukNP48pBAAAAAAAAwFMsDgIAAAAAAACeYnEQAAAAAAAA8BSLgwAAAAAAAICnWBwEAAAAAAAAPMXiIAAAAAAAAOCpjL6eQKoJGU/pUYM+a+bMGTLFjH3uiBozVj6stsvzet/OvQVm7Mnqwc7x+7ZXmTlr6xaZsXi8qesTAwaQjPRiM/ZfR37TOX79G580c0LRFjMW1EXsiTQ3Gwe0/98nGF1uxvLOdveHG3NyzZzjTmwzY1993d0fmqLvmTlAsoVCWWYsJz1wjkcj9unR67vdr5WS9Eok7Byf1OIel6TTj2g0Y2kTRpgxFRe6c2oazJTqv8fM2M/fcdf/My1vmjkfTT/GjJVlu+dX2lZv5jz54igzFou3mjEgkR3x183YW7Wfdo7HgpCZ80psg/1gQdw9HjLGJe2N5Jmxsr32+wA1us8BcjLs44Vk/1zubihJdt9oad1pxlbpT87xd7aPN3NqWzabsaao/X4ESGRP4xtm7Kmd7h4wOD7MzNmY4LwiFrhf0/dpm5kTRO3X+lCt/T4gfEyOczw3VGzm9A92t0Hf4spBAAAAAAAAwFMsDgIAAAAAAACeYnEQAAAAAAAA8BSLgwAAAAAAAICnWBwEAAAAAAAAPMXiIAAAAAAAAOCpjO4mPP3007rtttu0Zs0aVVVVacmSJTr33HM74pdcconuv//+TjkzZ87U0qVLD3uy/Ue6GRk/6PPO8Ve+aT/VGXOOMGPB8I/bsZxcM2YZ1Nxkxj6ya5dz/JuP2luvz7jx62assu5u53gQtJg56N/8qv+QGflE/sVm7IaHhriP9vRzZk77MxvMWN16ex4tTZnux0oLzJycvDYzlj/OPZ5x0lFmzuyHP2LG/nbCV5zj91ffYuYEajdj6HsDsQcEQcyMvRVxj0db7Nf5dxuzzNhTOxuc4+815Jk5Hy2zY3m7jQlKUoP7tTS+o85MeWnPSDP2+72/do63tteaOR8fPN6MFWe1Osez3C1SkrSs2n5u29r32InoFQOx/iVpb+ObZuzx6jOd4xOL7OPti71nxgIZr79B3MzJyrB7lPbZPSDY6+43mWl2T0k+e+4trTud4zta3e83JM4B+ruB2gPa22vM2B9qFjnHS8JHmzlpoW4voWhXw6t2MN1eC1BNfYI893Veo+OjzJS3Ery/kdW/4I1uXznY2NioKVOmaMGCBeZ9zjrrLFVVVXXc/vCHPxzWJAH0D9Q/4Dd6AOAv6h/wGz0ASG3dXvaeNWuWZs2alfA+4XBYZWVlhzwpAP0T9Q/4jR4A+Iv6B/xGDwBSW4985+CKFSs0bNgwjR8/Xpdddpn27t1r3jcajaqurq7TDcDA1Z36l+gBQKrhHADwF+cAgN/oAcDAlfTFwbPOOku/+c1vtHz5ct16662qrKzUrFmzFIu5v5Ni/vz5Kioq6riVl5cne0oAekl361+iBwCphHMAwF+cAwB+owcAA1v3v03zQ1xwwQUdfz7uuOM0efJkHX300VqxYoVOP/30g+4/b948XX311R1/r6uroykAA1R361+iBwCphHMAwF+cAwB+owcAA1uPfKz4g4466igNGTJEGza4d+MMh8MqLCzsdAOQGj6s/iV6AJDKOAcA/MU5AOA3egAwsCT9ysF/tm3bNu3du1fDhw/v6YfqNUV5x5ix587Kco6nf+tf7QNu3WqG0l58yc5ra7djlkz7nzwYVOQOXPQZM+dv7/3ZjI359ced49WRZ80cpJaBXP8hpZuxC0YatSIptHO3c3zXz98zczZWDzNjbfGQGUu2rPfizvGjXrPnPmxOmxn7yhh3j/rNLneflKQgfgh9Df1W/+gB9seZ/tH8rnM8nG3XXXZ6YMbWBc84x99qyjdzLq453oyN2LbPjKUNzXWOt2y0a3Jzk30O0NK6w/04ae7HkaQX2980Y8/vnegcf2BxsZnzq52/MmMYePpH/UuxWKMZ+2PkUef4+PqTzJzG1l0JHs39OhoEUTMjO8t+3QtG2OcHoYj75xoctntUKK03X3/dvdeeHVJNf+kBiX7rGlu2O8db2+vNnHjQ/Vppj9nfpRg02+cpQdkQO/bKNuf4mPwcMycUsd/fBOIc3HfdXhxsaGjotPq/adMmrV27ViUlJSopKdFNN92k2bNnq6ysTBs3btS3v/1tjR07VjNnzkzqxAH0Puof8Bs9APAX9Q/4jR4ApLZuLw6uXr1an/70pzv+/v53BMyZM0cLFy7UunXrdP/996u2tlYjRozQmWeeqe9///sKh8PJmzWAPkH9A36jBwD+ov4Bv9EDgNTW7cXBU089VUFgX5q7bNmyw5oQgP6L+gf8Rg8A/EX9A36jBwCprcc3JAEAAAAAAADQP7E4CAAAAAAAAHiqx3crTkWjMt278EpS8XemuQN//ruZ0/SCvXtRZGe2nRfNNGOWvJxWM1ZU1uIcz5n6lpmT9bWTzdjYX693jleL3YrR/6VnFJqxfzvjHTMW+a37/1zWbj3CzGmK2TuHHYp4krcCrN9aasamP2LvZHzSKe7xrDXFZk5La1NXpwUcts3N/3CO//j5r5g5eQnOnBqiVc7x5qh7N0RJerdhuhmbfsYoMxYMc+9gmNvytJmz+s/2OYAlHm82Y+/UPGTG/iOy1DkeCtlPYCzBbo7AobN3Aq1pcJ/jvpRRbea0xxO9TlkvwAl2CE3wmh2MPdqODRnsHD9j+PNmTtrb9m7FcfH6C/8Egfv9b1tbgp17Q4mur3LvWJ6oB7TX200gPUEPSG9397ZJv7d3Wk6rzjVjvAaDKwcBAAAAAAAAT7E4CAAAAAAAAHiKxUEAAAAAAADAUywOAgAAAAAAAJ5icRAAAAAAAADwFIuDAAAAAAAAgKcy+noC/VUowVNzUt6RduKaN53Du/7m3iZdkl7bcYQZa43b67chexam2D47K7vavR368ft2mjklaa+asalFw5zjz9VlmTlB0GrGgN6UHrJ/T999tcSMRaJh5/jbDdlmTnZaYMZG5UbNWEnYXS/RmN03dhvzk6StTZnu48Xt52LQenedS9LOZvfPHItT5+gfGls2Ocdv2/IzM6ckd5wZa2nd0+05bG1ON2OhvbV2Ypr79Txoc7+WS1JL3I7Z7P6USDzedEh5QG8KAvdrbLRtd6/NYXd9nhkr37XLTmx2v7c49qwGMyV3VakZq2uqtR8L8Eyg9gTBQ3kXbmveY687ZO/bZ08jw33+8KWPbTRzrnt3qBlraK4zY/ADVw4CAAAAAAAAnmJxEAAAAAAAAPAUi4MAAAAAAACAp1gcBAAAAAAAADzF4iAAAAAAAADgKRYHAQAAAAAAAE/Z+2bD1BYPzFjk77XO8WffG2XmNLTba7RFmbEuz+t9CaanQPbW61Utmc7xxnePMHOmL91uxl6N1LvnEHT/ZwJ6W7St2oz9ecsQMzY8O+4cf3mvXZgfH2LXZTRu94eaaJZzvD1IcLyYfbzdUXfeS3tbzJwR2Xlm7K/b3cdra99r5gD9QXt7xIztqnvhEI5o139mgv+mDar3mbFQmru+YrWtZs6UQTlmbGmtPQ8gdVm1mdxz1VCC8+/B+U123h67BygadQ4HTe1mSml4ohmra1pvPxaAD0jwZttkv9jnlNg1G9qT4Jw50uAczj/K7jcjsj9qxt5u3mg/FrzAlYMAAAAAAACAp1gcBAAAAAAAADzF4iAAAAAAAADgKRYHAQAAAAAAAE+xOAgAAAAAAAB4qtu7FT/99NO67bbbtGbNGlVVVWnJkiU699xzO+JBEOjGG2/UL3/5S9XW1uqkk07SwoULNW7cuGTOu8cFsncNerTxcTP25Xc/6Rx/fq97J2BJKsuxdzyqT7CTcXaaOy+WYAOlxpi9e1FTuzu2crf9axJps3dhXt32oBFht+KBypf638+ulSW7t5mxcellzvFX4xvMnKNbJ5ixjJDdO4qN3cxbEu1I3GrHdjW7m0dlywNmzrrNR5mx6vo1ZgwDkz89INmvU3Y/yTJeyyUpqK63j5jt7g1Bsz33E0rsnYy1yQ4Bkk/1n3yhtCwzlpuboC53JdituN1d60E0bqb827BJZmxezRIjcig7syIV0QMOXXpathkL2e1B2lNjxyLunc7jTXYPmFVwrBl7uybdiPDe3RfdvnKwsbFRU6ZM0YIFC5zxH/3oR/rJT36iu+++W88//7zy8vI0c+ZMtbS0HPZkAfQt6h/wGz0A8Bf1D/iNHgCktm5fOThr1izNmjXLGQuCQHfccYeuv/56nXPOOZKk3/zmNyotLdVf/vIXXXDBBYc3WwB9ivoH/EYPAPxF/QN+owcAqS2p3zm4adMmVVdXa8aMGR1jRUVFmjZtmlauXOnMiUajqqur63QDMPAcSv1L9AAgVXAOAPiLcwDAb/QAYOBL6uJgdXW1JKm0tLTTeGlpaUfsn82fP19FRUUdt/Ly8mROCUAvOZT6l+gBQKrgHADwF+cAgN/oAcDA1+e7Fc+bN0+RSKTjtnXr1r6eEoBeRA8A/EX9A36jBwB+owcA/UdSFwfLyvbv0rlz585O4zt37uyI/bNwOKzCwsJONwADz6HUv0QPAFIF5wCAvzgHAPxGDwAGvm5vSJLImDFjVFZWpuXLl+v444+XJNXV1en555/XZZddlsyH6lP10e1m7P+96d4CvDbYZ+Yc3z7IjJ1R1mbGirNaneNtcXvNd3fU3iv9kW0h5/jTbc+YOX/ZWGvG6ps2mjGkntSr/8CMrK39tRl7OzzaOZ6Rnm3mbKo/xoxVDLZ3eCvNbXKON7RmmjnZjblm7JEd7uM1tdj/i9vYssWMSbEEMaSa1OsByWS/LhdlxM1Y6NzpZix25FHO8Yxxr5k5syKrzFjGKyXO8fZ2+/wFeB/1n1hWRrEZi7Xb/SE2a4YZs2QcZfeAK2rsHnDDxsHO8bb2Pd2eA/xDD0gsO8tdX5IUa3S/B5ek2BmnmbFQY4NzPGvM22bOTbvXmrFf/mGkc7wp+p6Zg9TS7cXBhoYGbdiwoePvmzZt0tq1a1VSUqJRo0bpyiuv1A9+8AONGzdOY8aM0Xe/+12NGDFC5557bjLnDaAPUP+A3+gBgL+of8Bv9AAgtXV7cXD16tX69Kc/3fH3q6++WpI0Z84c3Xffffr2t7+txsZGff3rX1dtba0++clPaunSpcrOtq+aATAwUP+A3+gBgL+of8Bv9AAgtXV7cfDUU09VENgftwuFQrr55pt18803H9bEAPQ/1D/gN3oA4C/qH/AbPQBIbX2+WzEAAAAAAACAvsHiIAAAAAAAAOApFgcBAAAAAAAAT3X7OwchNUe3mbGl7Xc7x3PDpWbOmTkXmrF9rekJZpLlHG2O2Wu+exIcb3DYPb5j93MJ5pBI7BDzgIGrKfqeczwtLd/MGTnCrtmaVrtN56S7e0BdW6aZU99u94BxuQXO8acj7WYOgA+XEx5uxs6esMWMhXbYX+IeynfXa2jX3q5P7APu+sg3neNff/UHh3Q8AAeMyJtqxkpPt1+XQ9vt9xxBhjsvtLvGzomFzNh/j3X3gG+/xffHAYdrVPaJZizvdHudIF61w4yFGhvdgX21dk7I/s7I/xo9xzl+/dvfN3Mk+3gYeLhyEAAAAAAAAPAUi4MAAAAAAACAp1gcBAAAAAAAADzF4iAAAAAAAADgKRYHAQAAAAAAAE+xW3GSxWL1zvGh4U+ZOdXN9i4/Oen2DmaZxtJuc4KdyGpa7djeqHt34SMHzTRzNtc8ZsYAHDCm6HQzVtNq5+1LsFtxfoZ7V+L6djunvt3uAc3tcef40YPOMXM21jxsxgDs96WSS8xYe/tWO7F6nxkKFVQ7x4OdtWZO4C5xSVKa0RpOK7rKzHkycrt9QAAdLh36UTMWq95gxjLfs/tDKBx2jgc77B3L4+1mSEWZ7gZxatGVZs6KyB32AQF0uOKIo81YfIe9K3lagh6gFvcbiGCb3QNirfa1YcPC7h5wWoIewHlAauHKQQAAAAAAAMBTLA4CAAAAAAAAnmJxEAAAAAAAAPAUi4MAAAAAAACAp1gcBAAAAAAAADzF4iAAAAAAAADgqYy+nkDqCZyjccXMjIY297bhkrQ3mm7GstPda7stsZCZs6vZPT9Jaovb8wDQc1ra7bqsabP/D6eozd3C69vtvrGrxe4PAHpGesiuu+aWTDMW27zPPmamu/7j1fVmTtNe+7HSjPMXAIcvmuAUu2mrHSvcVG3GQmF3Pce3RMyc5nq7B1ji9AbgsDUmeH8efbfFjGUPqbIP2tbuHG7fbJ8HNDdmmbHWuHuO7WKNwBdcOQgAAAAAAAB4isVBAAAAAAAAwFMsDgIAAAAAAACeYnEQAAAAAAAA8BSLgwAAAAAAAICnkr44+L3vfU+hUKjTbcKECcl+GAD9EPUP+I0eAPiNHgD4i/oHBraMnjjoxIkT9cQTTxx4kIweeZgBJUiwBXhjW4JYu71+W9fm3m68JWbPo6k9MGOBEUo0d+CfUf/d1xyza6y+Ld2M1RqxhnZ3b9h/PLsHpHMtOZKAHtA9u+rzzNjQd3aZsZy03c7xlg0tZs623UPNWFx23wC6gx5wsETn5tu3F5uxnNf2mbH0Ivfz2vC2/WDbdw82Y+0BPQCHj/p3a0pwbr7z3XwzdkTB3m4/Vt1G+znfvq/IjFlvEWJK0MCQUnqkWjMyMlRWVtYThwbQz1H/gN/oAYDf6AGAv6h/YODqketE3nnnHY0YMUJHHXWUvvzlL2vLli098TAA+iHqH/AbPQDwGz0A8Bf1DwxcSb9ycNq0abrvvvs0fvx4VVVV6aabbtLJJ5+s1157TQUFBQfdPxqNKhqNdvy9rq4u2VMC0Eu6W/8SPQBIJZwDAH6jBwD+4n0AMLAlfXFw1qxZHX+ePHmypk2bptGjR+tPf/qTLr300oPuP3/+fN10003JngaAPtDd+pfoAUAq4RwA8Bs9APAX7wOAga3Hv36+uLhYxxxzjDZs2OCMz5s3T5FIpOO2devWnp4SgF7yYfUv0QOAVMY5AOA3egDgL94HAANLjy8ONjQ0aOPGjRo+fLgzHg6HVVhY2OkGIDV8WP1L9AAglXEOAPiNHgD4i/cBwMCS9I8VX3PNNTr77LM1evRo7dixQzfeeKPS09N14YUXJvuhUsa78Z1mLCdqN9PMtHTneFvc2IdcUqTV3op8fcAXxuLwUP+2NLnrVZI2tdSbsfzMIjOWleb+/53WuD2PugTB16PuXhSEEhwQ+AB6gFtju11Dz+2x3wgVbY6asWENDc7x2lr39zpJ0uaGPDO2rrbH/78YHqAHuNVE7XPz53aVmLH8N+wekJvb6hzfta/YzHm7zu4P6+tCZgzoCurfVuMuV0nSM9tLzdipadVmLD3NfW5RleA8YH1dvhnbaL8dgSeSvji4bds2XXjhhdq7d6+GDh2qT37yk1q1apWGDh2a7IcC0M9Q/4Df6AGA3+gBgL+of2BgS/ri4OLFi5N9SAADBPUP+I0eAPiNHgD4i/oHBjY+QwIAAAAAAAB4isVBAAAAAAAAwFMsDgIAAAAAAACeSvp3DsLtvZqlZixr0OfNWHW02IzFgpxuz2NrbJ8Za0yrcY4nmjuArnmn5iEzljXoy2asqinXjMXimc5xe09EaXOCnZHr09z94d2aRxMcEcCH+f2u/zZjaaHvmLFh4cFmrLzR3Rua2u1Tu+d2Z5mx9ZEm5/iTkdvNHABd88sd3zdjbfHrzdjgLHsjh6LMdud4fYIesHKvHXst4j4/eDpyp5kDoGtuf+9mM3Zpq90D8jLsHpCT7t6tONJm1/k/9tjXhr1eX+ccfy7yMzMHqYUrBwEAAAAAAABPsTgIAAAAAAAAeIrFQQAAAAAAAMBTLA4CAAAAAAAAnmJxEAAAAAAAAPAUi4MAAAAAAACAp+x9rtFr3ql5yIxlDfqyGdvbFnaOB4qZOc1pjWZsc81jZgxAz3m95vdmLFz8f8zYzkZ3D4grbuY0pNWZsUS9CEDP+O3O/zZj2enfNWM5GfnO8bZ4YOZUNbWasWWRH5sxAD3nvuofmLHMNLsHZKRlOsdj9imAtjW1mLEnI7fbiQB6zD1Vdg+QrjcjaaGQc7w9wXnAtuYmM/Z05M4E84APuHIQAAAAAAAA8BSLgwAAAAAAAICnWBwEAAAAAAAAPMXiIAAAAAAAAOApFgcBAAAAAAAAT7E4CAAAAAAAAHgqFASBvdd1H6irq1NRUZH2r1u6t+cG0F2BpLgikYgKCwv7ejIJ0QOAnjAwegD1D/SEgVH/Ej0A6Bn0AMBfXa9/rhwEAAAAAAAAPMXiIAAAAAAAAOApFgcBAAAAAAAAT7E4CAAAAAAAAHiKxUEAAAAAAADAUz22OLhgwQIdeeSRys7O1rRp0/TCCy/01EMB6Geof8Bv9ADAX9Q/4Dd6ADAw9cji4B//+EddffXVuvHGG/XSSy9pypQpmjlzpnbt2tUTDwegH6H+Ab/RAwB/Uf+A3+gBwMAVCoIgSPZBp02bphNOOEE/+9nPJEnxeFzl5eW64oordN111yXMraurU1FRkfavW4aSPTXAU4GkuCKRiAoLC3v0kQ6n/iV6ANAzBkYPoP6BnjAw6l+iBwA9gx4A+Kvr9Z/0KwdbW1u1Zs0azZgx48CDpKVpxowZWrly5UH3j0ajqqur63QDMDB1t/4legCQSjgHAPzFOQDgN3oAMLAlfXFwz549isViKi0t7TReWlqq6urqg+4/f/58FRUVddzKy8uTPSUAvaS79S/RA4BUwjkA4C/OAQC/0QOAga3PdyueN2+eIpFIx23r1q19PSUAvYgeAPiL+gf8Rg8A/EYPAPqPjGQfcMiQIUpPT9fOnTs7je/cuVNlZWUH3T8cDiscDnf8/cBXICb9qxABj+2vpx74itFOulv/Ej0A6B39swdQ/0Bv6J/1L9EDgN5BDwD81fX6T/qVg1lZWZo6daqWL1/eMRaPx7V8+XJVVFR8aH59ff3//mn/Fydy48YtGbf9zeBAffWMw61/iR7AjVvP3AZGD6D+uXHridvAqH+JHsCNW8/c6AHcuPl763r9J/3KQUm6+uqrNWfOHH384x/XiSeeqDvuuEONjY366le/+qG5I0aM0NatW1VQUKBQKKS6ujqVl5dr69atPb67Un/G83AAz8UBXX0ugiBQfX29RowY0eNzOpz6lzr3gPr6ev6t/xe/9wfwXByQaj2AcwAbz8V+PA8HpFr9S5wDWPi9P4Dn4gB6gD/4vT+A52K/nqj/HlkcPP/887V7927dcMMNqq6u1vHHH6+lS5ce9OWkLmlpaRo5cuRB44WFhV7/47+P5+EAnosDuvJcFBUV9cpcDqf+pc49IBQKSeLf+oN4Lg7guTggVXoA5wAfjudiP56HA1Kl/iXOAT4Mz8UBPBcH0AP8wXNxAM/Ffsms/x5ZHJSkyy+/XJdffnlPHR5AP0b9A36jBwD+ov4Bv9EDgIGpz3crBgAAAAAAANA3+v3iYDgc1o033thpFyMf8TwcwHNxQKo/F6n+83UHz8UBPBcHpPpzkeo/X3fwXOzH83BAqj8Xqf7zdQfPxQE8Fwek+nOR6j9fd/BcHMBzsV9PPA+hoKf3NAcAAAAAAADQL/X7KwcBAAAAAAAA9AwWBwEAAAAAAABPsTgIAAAAAAAAeIrFQQAAAAAAAMBT/XpxcMGCBTryyCOVnZ2tadOm6YUXXujrKfW4p59+WmeffbZGjBihUCikv/zlL53iQRDohhtu0PDhw5WTk6MZM2bonXfe6ZvJ9rD58+frhBNOUEFBgYYNG6Zzzz1X69ev73SflpYWzZ07V4MHD1Z+fr5mz56tnTt39tGMe8bChQs1efJkFRYWqrCwUBUVFfrb3/7WEU/l54Ae4G8PoP4P8LUHUP/+1r9ED/ggegA94H2+9ADq/wDqn/p/ny/1L9EDPqg3e0C/XRz84x//qKuvvlo33nijXnrpJU2ZMkUzZ87Url27+npqPaqxsVFTpkzRggULnPEf/ehH+slPfqK7775bzz//vPLy8jRz5ky1tLT08kx7XmVlpebOnatVq1bp8ccfV1tbm84880w1NjZ23Oeqq67So48+qgceeECVlZXasWOHPv/5z/fhrJNv5MiRuuWWW7RmzRqtXr1ap512ms455xy9/vrrklL3OaAH+N0DqP8DfOwB1L/f9S/RAz6IHkAPeJ8vPYD6P4D6p/7f50v9S/SAD+rVHhD0UyeeeGIwd+7cjr/HYrFgxIgRwfz58/twVr1LUrBkyZKOv8fj8aCsrCy47bbbOsZqa2uDcDgc/OEPf+iDGfauXbt2BZKCysrKIAj2/+yZmZnBAw880HGfN998M5AUrFy5sq+m2SsGDRoU/OpXv0rp54AeQA/4IOq/s1TvAdQ/9f/P6AGd0QNSHz3gAOq/M+o/9VH/ndEDOuupHtAvrxxsbW3VmjVrNGPGjI6xtLQ0zZgxQytXruzDmfWtTZs2qbq6utPzUlRUpGnTpnnxvEQiEUlSSUmJJGnNmjVqa2vr9HxMmDBBo0aNStnnIxaLafHixWpsbFRFRUXKPgf0ADefewD1v58PPYD6d/O5/iV6wPvoAQPzZ0oGn3sA9b8f9T8wf6Zk8Ln+JXrA+3q6B2Qkc7LJsmfPHsViMZWWlnYaLy0t1VtvvdVHs+p71dXVkuR8Xt6Ppap4PK4rr7xSJ510kiZNmiRp//ORlZWl4uLiTvdNxefj1VdfVUVFhVpaWpSfn68lS5bo2GOP1dq1a1PyOaAHuPnaA3yvf8mvHkD9u/la/xI9QKIHSPQAX3sA9U/9S9S/r/Uv0QOk3usB/XJxEPhnc+fO1WuvvaZnn322r6fSJ8aPH6+1a9cqEonowQcf1Jw5c1RZWdnX0wJ6he/1L9ED4Dd6AD0A/qL+qX/4jR7Qez2gX36seMiQIUpPTz9ol5WdO3eqrKysj2bV997/2X17Xi6//HL99a9/1VNPPaWRI0d2jJeVlam1tVW1tbWd7p+Kz0dWVpbGjh2rqVOnav78+ZoyZYruvPPOlH0O6AFuPvYA6n8/n3oA9e/mY/1L9ID30QMG9s+UDD72AOp/P+p/YP9MyeBj/Uv0gPf1Vg/ol4uDWVlZmjp1qpYvX94xFo/HtXz5clVUVPThzPrWmDFjVFZW1ul5qaur0/PPP5+Sz0sQBLr88su1ZMkSPfnkkxozZkyn+NSpU5WZmdnp+Vi/fr22bNmSks/HB8XjcUWj0ZR9DugBbj71AOo/sVTuAdS/m0/1L9EDPgw9wD8+9QDqPzHq3z8+1b9ED/gwPdYDkrRhStItXrw4CIfDwX333Re88cYbwde//vWguLg4qK6u7uup9aj6+vrg5ZdfDl5++eVAUvA///M/wcsvvxy89957QRAEwS233BIUFxcHDz/8cLBu3brgnHPOCcaMGRM0Nzf38cyT77LLLguKioqCFStWBFVVVR23pqamjvv8+7//ezBq1KjgySefDFavXh1UVFQEFRUVfTjr5LvuuuuCysrKYNOmTcG6deuC6667LgiFQsHf//73IAhS9zmgB/jdA6j/A3zsAdS/3/UfBPSAD6IH0AN86wHU/wHUP/XvW/0HAT3gg3qzB/TbxcEgCIKf/vSnwahRo4KsrKzgxBNPDFatWtXXU+pxTz31VCDpoNucOXOCINi/jfl3v/vdoLS0NAiHw8Hpp58erF+/vm8n3UNcz4Ok4N577+24T3Nzc/DNb34zGDRoUJCbmxv867/+a1BVVdV3k+4B//Zv/xaMHj06yMrKCoYOHRqcfvrpHc0gCFL7OaAH+NsDqP8DfO0B1L+/9R8E9IAPogfQA3zrAdT/AdQ/9e9b/QcBPeCDerMHhIIgCLp/vSEAAAAAAACAga5ffucgAAAAAAAAgJ7H4iAAAAAAAADgKRYHAQAAAAAAAE+xOAgAAAAAAAB4isVBAAAAAAAAwFMsDgIAAAAAAACeYnEQAAAAAAAA8BSLgwAAAAAAAICnWBwEAAAAAAAAPMXiIAAAAAAAAOApFgcBAAAAAAAAT7E4CAAAAAAAAHjq/wdKoe3ehSxwBAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1600x500 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "steps = [[ smoke.values, velocity.values.vector[0], velocity.values.vector[1] ]]\n",
    "for time_step in range(20):\n",
    "  if time_step<3 or time_step%10==0: \n",
    "    print('Computing time step %d' % time_step)\n",
    "  velocity, smoke, pressure = step(velocity, smoke, pressure, dt=DT)\n",
    "  if time_step%5==0:\n",
    "    steps.append( [smoke.values, velocity.values.vector[0], velocity.values.vector[1]] )\n",
    "\n",
    "fig, axes = pylab.subplots(1, len(steps), figsize=(16, 5))\n",
    "for i in range(len(steps)):\n",
    "    axes[i].imshow(steps[i][0].numpy('y,x'), origin='lower', cmap='magma')\n",
    "    axes[i].set_title(f\"d at t={i*5}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also take a look at the velocities. The `steps` list above already stores `vector[0]` and `vector[1]` components of the velocities as numpy arrays, which we can show next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig, axes = pylab.subplots(1, len(steps), figsize=(16, 5))\n",
    "# for i in range(len(steps)):\n",
    "#     axes[i].imshow(steps[i][1].numpy('y,x'), origin='lower', cmap='magma')\n",
    "#     axes[i].set_title(f\"u_x at t={i*5}\")\n",
    "    \n",
    "# fig, axes = pylab.subplots(1, len(steps), figsize=(16, 5))\n",
    "# for i in range(len(steps)):\n",
    "#     axes[i].imshow(steps[i][2].numpy('y,x'), origin='lower', cmap='magma')\n",
    "#     axes[i].set_title(f\"u_y at t={i*5}\")\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ooqVxCPM8PXl"
   },
   "source": [
    "It looks simple here, but this simulation setup is a powerful tool. The simulation could easily be extended to more complex cases or 3D, and it is already fully compatible with backpropagation pipelines of deep learning frameworks. \n",
    "\n",
    "In the next chapters we'll show how to use these simulations for training NNs, and how to steer and modify them via trained NNs. This will illustrate how much we can improve the training process by having a solver in the loop, especially when the solver is _differentiable_. Before moving to these more complex training processes, we will cover a simpler supervised approach in the next chapter. This is very fundamental: even when aiming for advanced physics-based learning setups, a working supervised training is always the first step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "\n",
    "You could create a variety of nice fluid simulations based on this setup. E.g., try changing the spatial resolution, the buoyancy factors, and the overall length of the simulation run."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "overview-ns-forw.ipynb",
   "provenance": []
  },
  "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}