{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 0,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-2ad2c23f9a0a820a",
     "locked": true,
     "schema_version": 1,
     "solution": false
    },
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Set-up-the-boundary-values\" data-toc-modified-id=\"Set-up-the-boundary-values-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Set up the boundary values</a></span></li><li><span><a href=\"#Set-up-the-diffusivity-matrix\" data-toc-modified-id=\"Set-up-the-diffusivity-matrix-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Set up the diffusivity matrix</a></span></li><li><span><a href=\"#Steady-state-the-old-fashioned-way\" data-toc-modified-id=\"Steady-state-the-old-fashioned-way-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Steady state the old-fashioned way</a></span></li><li><span><a href=\"#Repeat-the-steady-state-calculation-using-a-sparse-matrix\" data-toc-modified-id=\"Repeat-the-steady-state-calculation-using-a-sparse-matrix-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Repeat the steady state calculation using a sparse matrix</a></span></li><li><span><a href=\"#Now-do-the-2d-problem\" data-toc-modified-id=\"Now-do-the-2d-problem-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Now do the 2d problem</a></span><ul class=\"toc-item\"><li><span><a href=\"#Set-up-the-problem-domain-and-build-A-and-b\" data-toc-modified-id=\"Set-up-the-problem-domain-and-build-A-and-b-5.1\"><span class=\"toc-item-num\">5.1&nbsp;&nbsp;</span>Set up the problem domain and build A and b</a></span></li></ul></li><li><span><a href=\"#Compare-dense-and-sparse\" data-toc-modified-id=\"Compare-dense-and-sparse-6\"><span class=\"toc-item-num\">6&nbsp;&nbsp;</span>Compare dense and sparse</a></span></li><li><span><a href=\"#Plot-your-saved-snapshots-on-a-grid\" data-toc-modified-id=\"Plot-your-saved-snapshots-on-a-grid-7\"><span class=\"toc-item-num\">7&nbsp;&nbsp;</span>Plot your saved snapshots on a grid</a></span></li><li><span><a href=\"#Conclusions\" data-toc-modified-id=\"Conclusions-8\"><span class=\"toc-item-num\">8&nbsp;&nbsp;</span>Conclusions</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introducing sparse matrices\n",
    "\n",
    "This is a rerun of the 2D_Transient_Assignment with two differences:\n",
    "\n",
    "1. It imports common functions from a new python module called module_2D.py\n",
    "\n",
    "1. It uses the scipy.sparse compressed storage format for rows (csr) with\n",
    "   the sparse solver spsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "lines_to_next_cell": 2,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-e58a84ffce31f93a",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "It looks like I've been imported by another python script\n",
      "My file location is: \n",
      "/Users/phil/repos/eosc213_students/notebooks/2D_Assignment_sparse/module_2D.py\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "import matplotlib.cm as cmap\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from module_2D import Boundary_Def\n",
    "from module_2D import build_2D_matrix\n",
    "from module_2D import mat2vec\n",
    "from module_2D import Problem_Def\n",
    "from module_2D import vec2mat\n",
    "from mpl_toolkits.axes_grid1 import AxesGrid\n",
    "from scipy.sparse import csr_matrix\n",
    "from scipy.sparse.linalg import spsolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the boundary values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "c0 = 1  # mg/L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-f97ab996f12bd2f2",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# Here we create 4 boundaries, west has a constant concentration at c0, east has a constant boundary at 0;\n",
    "west = Boundary_Def(\"const\", val=c0)\n",
    "east = Boundary_Def(\"const\", val=0)\n",
    "\n",
    "# For 1D problem, the used boundaries are west and east.\n",
    "\n",
    "# The other south and north boundaries have a zero flux (impermeable)\n",
    "\n",
    "north = Boundary_Def(\"flux\", val=0)\n",
    "south = Boundary_Def(\"flux\", val=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-ea536fa5b54285e2",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "bc_dict = {\"west\": west, \"north\": north, \"east\": east, \"south\": south}\n",
    "# The latter array bc_dict will be sent to the different functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-7bfb0e282680103b",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Set up the diffusivity matrix\n",
    "\n",
    "Going from 1 to 2 dimensions changes nothing conceptually. There are, however a couple of changes required for the coding perspective. Indeed, whether the problem is 1D or 2D or 3D, the stucture of the system of equation Ac = b is the same. Matrix $A$ will always be a $n \\times n$ matrix, while $c$ and $b$ will always be column vector of size $n$. In 2D, $n = n_x \\times n_y$, while in 3D, it will be $n = n_x \\times n_y \\times n_z$.  The individual equation for a cell still\n",
    "produces a single row in the A and b matrices, but in 2D that cell has 4 neighbours instead\n",
    "of 2, and 3D it has 6 neighbors intead of 4.\n",
    "\n",
    "However, the fact is that in every case, the solution is stored in one vector, representing either a 1/2/3D solution. For these higher dimension problems, a two-way conversion between vector and matrix is required. To plot the 2D result, for example, we will use colourmap plots, which require the solution to be plotted to be represented as a 2D array (matrix).\n",
    "\n",
    "The function vec2mat(...) (specifically $vector\\ to\\ matrix$) does this: it converts a vector into the relevant 2D matrix, using n_x and n_y.\n",
    "\n",
    "The reverse function is usually required to initialize the initial condition. It is mac2vec(...). These two functions are defined here below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "decreasing_factor = 0.1  # Feel free to change if you want to see the impact\n",
    "# (you can go higher than 1 ... But be careful, if diffusion speeds up significantly,\n",
    "# the accuracy with respect to the chosen timestep might not be so good if you speed things up! )\n",
    "# Initial value is 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "lines_to_next_cell": 2,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-42d47ff484f9a502",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "Diff = 2e-9 * 100 * 24 * 3600  # dm²/day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "lines_to_next_cell": 2,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-b9bd7eaf90619a64",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x115c46ac8>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAD8CAYAAABTjp5OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAG3tJREFUeJzt3X+QXWWd5/H3h07CL+VXopQm2UrcZHSaOIikAsqUhcaVsDo2tRtqQ6mbcVOVKgv8OZYbZgqoZaRWqtxBZkCnuiCaRYuAGafsdbPGKYNlsbUbaBAn6cSMvSGShmgIhMgPQxL47B/nJNzc3F/p26T7NJ9XVRfnPOc5535vdfLJw3PPuY9sExERE98p411ARER0JoEdEVERCeyIiIpIYEdEVEQCOyKiIhLYEREV0TawJa2WtEfSlpq28yT9k6Rfl/899/UtMyJifElaImm7pGFJqxocP1XSfeXxTZLmlO3TJT0g6QVJd9Sdc42kzZL+WdKPJc1oVUMnI+zvAEvq2lYBP7U9H/hpuR8RMSlJ6gHuBK4EeoFrJPXWdVsB7LM9D7gNuLVsPwDcAHy57ppTgNuBD9r+E+Cfgeta1dE2sG3/HHi2rrkPWFNurwGuanediIgKWwQM295h+yCwliIHa9Xm4jpgsSTZftH2gxTBXUvlz5mSBJwFPNWqiCmjLP5827sBbO+W9NZmHSWtBFYCnHnGqRe/61+/fZQvGRFvJI9sfnyv7beM9vwrLr/Qzzz7fKevNcSxgdpvu79mfyawq2Z/BLik7jJH+9g+LGk/MB3Y2+g1bR+S9BlgM/Ai8Gvg2lZ1jjawO1a+6X6AhX/yDj+0/quv90tGxCTQM/sTv+nm/GeefZ5O86Zn9icO2F7YoosatNV/r0cnfV7rLE0FPgNcBOwA/g64Hmha9GjvEvmdpLeVL/o2YM8orxMRUQUjwOya/VkcP31xtE85P302x08n13oPgO3/5+JLne4H3t+qiNEG9gCwvNxeDvxwlNeJiKiCh4H5kuZKmgYso8jBWrW5uBTY6Nbfrvck0CvpyLTPvwG2tSqi7ZSIpHuBy4EZkkaAm4CvAfdLWgE8AVzd7joREVVVzklfB2wAeoDVtock3QwM2h4A7gbukTRMMbJeduR8STspPlScJukq4CO2t0r6L8DPJR0CfgP8eas62ga27WuaHFrc7tyIiMnC9npgfV3bjTXbB2gyeLU9p0n73wN/32kNedIxIqIiEtgRERWRwI6IqIgEdkRERSSwIyIqIoEdEVERCeyIiIpIYEdEVEQCOyKiIhLYEREVkcCOiKiIBHZEREUksCMiKiKBHRFREQnsiIiKSGBHRHRA0hJJ2yUNS1rV4Pipku4rj2+SNKdsny7pAUkvSLqj7pxpkvol/YukX0n6961qeN0X4Y2IqDpJPcCdFMt4jQAPSxqwvbWm2wpgn+15kpYBtwL/gWI19huABeVPrb8C9tj+I0mnAOe1qiMj7IiI9hYBw7Z32D4IrAX66vr0AWvK7XXAYkmy/aLtBymCu95/Av4rgO1Xbe9tVUQCOyKivZnArpr9kbKtYR/bh4H9wPRmF5R0Trn515IelfR9See3KiJTIhExOb30B/TI5k57z5A0WLPfb7u/Zl8NzqlfEb2TPrWmALOA/237S5K+BHwd+FSrEyIi3uj22l7Y4vgIMLtmfxbwVJM+I5KmAGdTrJ7ezDPAS8A/lvvfp5gHbypTIhER7T0MzJc0V9I0YBkwUNdnAFhebi8FNtpuOsIuj/0P4PKyaTGwtVl/yAg7IqIt24clXQdsAHqA1baHJN0MDNoeAO4G7pE0TDGyXnbkfEk7gbOAaZKuAj5S3mHyn8tzvgE8DXy6VR0J7IiIDtheD6yva7uxZvsAcHWTc+c0af8N8IFOa8iUSERERSSwIyIqIoEdEVERCeyIiIpIYEdEVEQCOyKiIhLYEREVkcCOiKiIBHZEREUksCMiKiKBHRFREV0FtqQvShqStEXSvZJOG6vCIiLiWKMObEkzgc8BC20voPgGq2Wtz4qIiNHqdkpkCnB6+WXdZ3D8F3pHRMQYGXVg236SYjmbJ4DdwH7bP6nvJ2mlpEFJg08/+/zoK42IeIPrZkrkXIpVgucCbwfOlPTJ+n62+20vtL3wLee9efSVRkS8wXUzJfJh4HHbT9s+BPwAeP/YlBURMbFIWiJpu6RhSasaHD9V0n3l8U2S5pTt0yU9IOkFSXc0ufaApC3taugmsJ8ALpV0hiRRrEe2rYvrRURMSJJ6gDuBK4Fe4BpJvXXdVgD7bM8DbgNuLdsPADcAX25y7X8HvNBJHd3MYW8C1gGPApvLa/W3PCkiopoWAcO2d9g+CKylmBKu1QesKbfXAYslyfaLth+kCO5jSHoT8CXgq50U0dWajrZvAm7q5hoREa8Hv3iQQ5ue7LT7DEmDNfv9tmsHoDOBXTX7I8Alddc42qdctHc/MB3Y2+J1/xr4b8BLnRSZRXgjImCv7YUtjqtBm0fR57XO0nuAeba/eGS+u508mh4R0d4IMLtmfxbHP3dytE/5bMrZwLMtrvk+4GJJO4EHgT+S9LNWRSSwIyLaexiYL2mupGkUT3UP1PUZAJaX20uBjbabjrBtf8v2223PAf4U+Bfbl7cqIlMiERFtlHPS1wEbKL6GY7XtIUk3A4O2B4C7gXskDVOMrI9+VUc5ij4LmCbpKuAjtreeaB0J7IiIDtheD6yva7uxZvsAcHWTc+e0ufZOYEG7GjIlEhFRERlhV4ge2TzeJUSHfPG7x7uEmIQywq6IhHW16JHN+Z3FmMsIuwL0yOYTeQAgJoipl8xEj2zOaDvGTAJ7AjsyQjsS1r/bcvp4lhMn6HyeTGjHmEpgT1CNwnpoz/TxLClOwAVvfYbfbTmd8yl+f1PJvHZ0L4E9AdWG9ZFRdcK6Wob2TH8ttBf8gUObnmRqeSzBHaOVwJ5gauer68P6l89ljeNqKf+R3fLMa6GdKZLoQgJ7gmg2Xz20Z/rRoN6y79XxKS5O2IJzT+GXz53GheccKP7B3fIMkHnt6E4CewJoNV+dsK6m135fr4V25rWjWwnscdZqvro2rB87tHNc6ovRe8/UOWVwH5nKajyvndCOTiWwx1H9yLpWo5H18EsPnJzComvzzvggjx3ayXumzjmmvf7DSCDTI9GxPOkYEVERCeyIiIpIYEdEVEQCOyKiA5KWSNouaVjSqgbHT5V0X3l805F1GiVNl/SApBck3VHT/wxJ/1PSryQNSfpauxoS2BERbUjqAe4ErgR6gWsk9dZ1WwHssz0PuA24tWw/ANwAfLnBpb9u+13ARcBlkq5sVUcCOyKivUXAsO0dtg8Ca4G+uj59wJpyex2wWJJsv2j7QYrgPsr2S7YfKLcPAo9SLO7bVG7ri4hJ6dAfTjmRb7icIWmwZr/fdn/N/kxgV83+CHBJ3TWO9inXgNxP8f0Ee9u9uKRzgD8Dbm/VL4EdEQF7bS9scVwN2upXRO+kz/EXlqYA9wJ/a3tHq76ZEomIaG8EmF2zPwt4qlmfMoTPplg9vZ1+4Ne2v9GuYwI7IqK9h4H5kuZKmgYsAwbq+gwAy8vtpcBG2y1H2JK+ShHsX+ikiEyJRES0Uc5JXwdsAHqA1baHJN0MDNoeAO4G7pE0TDGyXnbkfEk7gbOAaZKuAj4C/B74K+BXwKOSAO6wfVezOhLYEREdsL0eWF/XdmPN9gHg6ibnzmly2Ubz3k1lSiQioiIS2BERFZHAjoioiAR2RERFJLAjIiqiq8CWdI6kdeW3TW2T9L6xKiwiIo7V7W19twM/tr20vJn8jDGoKSIiGhh1YEs6C/gA8Odw9NumDo5NWRERUa+bKZF3AE8D35b0C0l3STqzvpOklZIGJQ0+/ezzXbxcRMQbWzeBPQV4L/At2xcBLwLHrcJgu9/2QtsL33Lem7t4uYiIN7ZuAnsEGLG9qdxfRxHgERHxOhh1YNv+LbBL0jvLpsXA1jGpKiIijtPtXSKfBb5X3iGyA/h09yVFREQjXQW27ceAVqs0RETEGMmTjhERFZHAjojogKQlkrZLGpZ03B1xkk6VdF95fJOkOWX7dEkPSHpB0h1151wsaXN5zt+qXMWgmQR2REQbknqAO4ErgV7gGkm9dd1WAPtszwNuA24t2w8ANwBfbnDpbwErgfnlz5JWdSSwIyLaWwQM295RPtW9Fuir69MHrCm31wGLJcn2i7YfpAjuoyS9DTjL9v8p137878BVrYrIEmERMSkdODyFoT3TO+0+Q9JgzX6/7f6a/ZnArpr9EeCSumsc7VOuAbkfmA7sbfKaM8vr1F5zZqsiE9gREbDXdqs73hrNLdeviN5Jn276Z0okIqIDI8Dsmv1ZwFPN+kiaApxNsXp6q2vOanPNYySwIyLaexiYL2lu+aDgMmCgrs8AsLzcXgpsLOemG7K9G3he0qXl3SH/EfhhqyIyJRIR0UY5J30dsAHoAVbbHpJ0MzBoewC4G7hH0jDFyHrZkfMl7QTOAqZJugr4iO2twGeA7wCnA/+r/GkqgR0R0QHb64H1dW031mwfAK5ucu6cJu2DwIJOa8iUSERERSSwIyIqIoEdEVERCeyIiIpIYEdEVEQCOyKiIhLYEREVkcCOiKiIBHZEREUksCMiKiKBHRFREQnsiIiKSGBHRFREAjsioiIS2BERFZHAjojogKQlkrZLGpa0qsHxUyXdVx7fJGlOzbHry/btkq6oaf+ipCFJWyTdK+m0VjUksCMi2pDUA9wJXAn0AtdI6q3rtgLYZ3secBtwa3luL8XqMxcAS4BvSuqRNBP4HLDQ9gKKlWyW0UICOyKivUXAsO0dtg8Ca4G+uj59wJpyex2wuFyrsQ9Ya/tl248Dw+X1oFj16/Ry0d4zyCK8ERFtzZA0WPOzsu74TGBXzf5I2dawj+3DwH5gerNzbT8JfB14AtgN7Lf9k1ZFZk3HiJiUXjp8Cr98ruWUcK29the2OK4GbfUrojfr07Bd0rkUo++5wHPA9yV90vZ3mxWREXZERHsjwOya/VkcP31xtE85xXE2xerpzc79MPC47adtHwJ+ALy/VREJ7IiI9h4G5kuaK2kaxYeDA3V9BoDl5fZSYKNtl+3LyrtI5gLzgYcopkIulXRGOde9GNjWqohMiUREtGH7sKTrgA0Ud3Ostj0k6WZg0PYAcDdwj6RhipH1svLcIUn3A1uBw8C1tl8BNklaBzxatv8C6G9VRwI7IqIDttcD6+vabqzZPgBc3eTcW4BbGrTfBNzUaQ2ZEomIqIiuA7u8AfwXkn40FgVFRERjYzHC/jxtJsojIqJ7XQW2pFnAR4G7xqaciIhoptsR9jeArwCvNusgaeWRp4eefvb5Ll8uIuKNa9SBLeljwB7bj7TqZ7vf9kLbC99y3ptH+3IREW943YywLwM+LmknxRehfEhS00cqIyKiO6MObNvX255lew7FDeIbbX9yzCqLiIhj5D7siIiKGJMnHW3/DPjZWFwrIiIaywg7IqIiEtgRERWRwI6IqIgEdkRERSSwIyIqIoEdEdEBSUskbZc0LGlVg+OnSrqvPL5J0pyaY9eX7dslXVHTfo6kdZJ+JWmbpPe1qiGBHRHRhqQe4E7gSqAXuEZSb123FcA+2/OA24Bby3N7KR4uvABYAnyzvB7A7cCPbb8LuJA233yawI6IaG8RMGx7h+2DFF/H0VfXpw9YU26vAxaXazX2AWttv2z7cWAYWCTpLOADFEuLYfug7edaFZHAjoiAGUe+VbT8WVl3fCawq2Z/pGxr2Mf2YWA/ML3Fue8Anga+XS4Cc5ekM1sVmTUdI2JS+sMrZsu+pt/8XG+v7YUtjqtBmzvs06x9CvBe4LO2N0m6HVgF3NCsiIywIyLaGwFm1+zPAp5q1kfSFOBsitXTm507AozY3lS2r6MI8KYS2BER7T0MzJc0V9I0ig8RB+r6DADLy+2lFN9g6rJ9WXkXyVxgPvCQ7d8CuyS9szxnMbC1VRGZEomIaMP2YUnXARuAHmC17SFJNwODtgcoPjy8R9Iwxch6WXnukKT7KcL4MHCt7VfKS38W+F75j8AO4NOt6khgR0R0wPZ6YH1d24012weAq5ucewtwS4P2x4BWc+fHyJRIRERFZIQ9jnzxu9Ejm5l6yUwObXrymGMXnnOAXz532jFt88744MksL7r0nqlzjmu74K3PAHD+gj+c5GpiMsgIe5z54ncDMPWSmUf/Eh/5S33hOQdYcO4pLDj3lIZ/+WPiOvL7WnBu8VfswnMOHBfWUy+ZiS9+99E/AxHtZIQ9AdSOtM+nHGlveYahPdOPjrQXnHsK7JszrnXGiVlw7ilceM4B4NiR9dRLiuctEtRxohLYE8TRkTZwaNOTxShsyzPl0emvhXZURsI6xloCe4Lxxe8+JrR/t+X08i/89PEuLUbhgrc+c8wUCCSsY/QS2BPQkdAGOJ8nj4b20J6EdlU0mq+GhHV0J4E9QTWb147qqP9wMaJbCewJrP62v9wKVi0ZVcdYS2BPcPUfRkY1JKzj9ZDArojaee2ohoR1jLUEdoUkACLe2HJjb0RERSSwIyIqIoEdEVERCeyIiA5IWiJpu6RhSasaHD9V0n3l8U2S5tQcu75s3y7pirrzespFeH/UroYEdkREG5J6gDuBK4Fe4BpJvXXdVgD7bM8DbgNuLc/tpVh95gJgCfDN8npHfB7Y1kkdCeyIiPYWAcO2d9g+CKwF+ur69AFryu11wGJJKtvX2n7Z9uPAcHk9JM0CPgrc1UkRCeyICJghabDmZ2Xd8ZnArpr9kbKtYR/bh4H9FN/a1urcbwBfAV7tpMjchx0Rk9JLPshjh3Z22n2v7VZrK6pBmzvs07Bd0seAPbYfkXR5J0WOeoQtabakByRtkzQk6fOjvVZExAQ3Asyu2Z8FPNWsj6QpwNkUq6c3O/cy4OOSdlJMsXxI0ndbFdHNlMhh4C9s/zFwKXBtg0n4iIjJ4GFgvqS5kqZRfIg4UNdnAFhebi8FNtp22b6svItkLjAfeMj29bZn2Z5TXm+j7U+2KmLUUyK2dwO7y+3nJW2jmJfZOtprRkRMRLYPS7oO2AD0AKttD0m6GRi0PQDcDdwjaZhiZL2sPHdI0v0U2XgYuNb2K6OpY0zmsMv7DS8CNjU4thJYCfCvZs4Yi5eLiDjpbK8H1te13VizfQC4usm5twC3tLj2z4Cftauh67tEJL0J+AfgC7Z/36CQftsLbS98y3lv7vblIiLesLoKbElTKcL6e7Z/MDYlRUREI93cJSKKOZtttv9m7EqKiIhGuhlhXwZ8iuJWlMfKn387RnVFRESdbu4SeZDGN4RHRMTrII+mR0RURAI7IqIiEtgRERWRwI6IqIgEdkRERSSwIyIqIoEdEVERCeyIiIpIYEdEVEQCOyKiIhLYEREdkLRE0nZJw5JWNTh+qqT7yuObynUCjhy7vmzfLumKsu2El1lMYEdEtCGpB7gTuBLoBa5psCTiCmCf7XnAbcCt5bm9FKvPXAAsAb5ZXu+El1lMYEdEtLcIGLa9w/ZBikVz++r69AFryu11wOLya6j7gLW2X7b9ODAMLLK92/ajUCyzCBxZZrGpBHZERHszgV01+yMcH65H+9g+DOwHpndybqtlFmuNyZqOERETzcuvPs/wSw902n2GpMGa/X7b/TX7jb5K2nX7zfq0PLfdMou1EtgREbDX9sIWx0eA2TX7s4CnmvQZkTQFOJti9fSm557oMouZEomIaO9hYL6kuZKmUXyIOFDXZwBYXm4vBTbadtm+rLyLZC4wH3hoNMssZoQdEdGG7cOSrgM2AD3AattDkm4GBm0PUITvPZKGKUbWy8pzhyTdD2yluDPkWtuvSPpTimUWN0t6rHypv7S9vlkdCeyIiA6UQbq+ru3Gmu0DwNVNzr0FuKWu7YSXWcyUSERERSSwIyIqIoEdEVERCeyIiIpIYEdEVEQCOyKiIhLYEREVkcCOiKiIBHZEREUksCMiKiKBHRFREQnsiIiKSGBHRFREAjsioiIS2BERFZHAjoioiK4CW9ISSdslDUtaNVZFRURMNO3yrlwC7L7y+KZyJfQjx64v27dLuqLTa9YbdWBL6gHuBK4EeoFrJPWO9noRERNVh3m3Athnex5wG3BreW4vxXJhFwBLgG9K6hlNhnYzwl4EDNveYfsgsBbo6+J6ERETVSd51wesKbfXAYvLhXb7gLW2X7b9ODBcXu+EM7SbNR1nArtq9keAS+o7SVoJrCx3X+6Z/YktXbzmRDQD2DveRbwO8r6qYzK+J4B3dnOyfXDDgZefmNFh99MkDdbs99vur9nvJO+O9ikX7d0PTC/b/2/duTPL7bYZWqubwG60eKSPayjedD+ApEHbC7t4zQlnMr4nyPuqksn4nqB4X92cb3vJWNVCZ3nXrE+z9kYzHMdlaK1upkRGgNk1+7OAp7q4XkTERNVJ3h3tI2kKcDbwbItzTzhDuwnsh4H5kuZKmkYxqT7QxfUiIiaqTvJuAFhebi8FNtp22b6svItkLjAfeKjDax5j1FMi5RzNdcAGoAdYbXuozWn9bY5X0WR8T5D3VSWT8T3BBHpfzfJO0s3AoO0B4G7gHknDFCPrZeW5Q5LuB7YCh4Frbb8CcKIZquIfgIiImOjypGNEREUksCMiKuKkBPZkfIRd0mxJD0jaJmlI0ufHu6axUj6F9QtJPxrvWsaKpHMkrZP0q/J39r7xrmksSPpi+edvi6R7JZ023jWdKEmrJe2RtKWm7TxJ/yTp1+V/zx3PGieK1z2wJ/Ej7IeBv7D9x8ClwLWT5H0BfB7YNt5FjLHbgR/bfhdwIZPg/UmaCXwOWGh7AcUHV8vGt6pR+Q7FI9u1VgE/tT0f+Gm5/4Z3MkbYk/IRdtu7bT9abj9PEQAzW5818UmaBXwUuGu8axkrks4CPkDxKT62D9p+bnyrGjNTgNPL+37PoILPQtj+OcVdFbVqH/NeA1x1UouaoE5GYDd6pLPywVar/Faui4BN41vJmPgG8BXg1fEuZAy9A3ga+HY51XOXpDPHu6hu2X4S+DrwBLAb2G/7J+Nb1Zg53/ZuKAZHwFvHuZ4J4WQEdkePsFeVpDcB/wB8wfbvx7uebkj6GLDH9iPjXcsYmwK8F/iW7YuAF5kE/4tdzuv2AXOBtwNnSvrk+FYVr6eTEdiT9hF2SVMpwvp7tn8w3vWMgcuAj0vaSTF19SFJ3x3fksbECDBi+8j/Aa2jCPCq+zDwuO2nbR8CfgC8f5xrGiu/k/Q2gPK/e8a5ngnhZAT2pHyEvfzaxLuBbbb/ZrzrGQu2r7c9y/Ycit/TRtuVH7HZ/i2wS9KRb39bTPHUWdU9AVwq6Yzyz+NiJsGHqaXax7yXAz8cx1omjG6+ra8jo3yEvQouAz4FbJb0WNn2l7bXj2NN0dxnge+Vg4YdwKfHuZ6u2d4kaR3wKMVdS79gAj3O3SlJ9wKXAzMkjQA3AV8D7pe0guIfpqvHr8KJI4+mR0RURJ50jIioiAR2RERFJLAjIioigR0RUREJ7IiIikhgR0RURAI7IqIi/j/yTWkUcJCrVwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here we define the initial condition, and the diffusion matrix for the 2D problem\n",
    "\n",
    "\n",
    "def make_D_matrix(the_problem, Diff, decreasing_factor):\n",
    "    Diff_low = Diff * decreasing_factor\n",
    "    n_x, n_y = the_problem.nx, the_problem.ny\n",
    "    width_x, width_y = the_problem.wx, the_problem.wy\n",
    "    D_matrix = Diff * np.ones((n_y, n_x))\n",
    "\n",
    "    c_init = np.zeros((n_y, n_x))\n",
    "\n",
    "    x = np.linspace(0, width_x, n_x)\n",
    "    y = np.linspace(0, width_y, n_y)\n",
    "\n",
    "    for i in range(n_y):\n",
    "        for j in range(n_x):\n",
    "            if j == 0:\n",
    "                c_init[i, j] = c0  # Initial condition\n",
    "    #\n",
    "    # overwrite the center of the image with a low diffusivity\n",
    "    #\n",
    "\n",
    "    for i in range(n_y):\n",
    "        for j in range(n_x):\n",
    "            if (\n",
    "                abs(x[j] - width_x / 2) <= 0.2 * width_x\n",
    "                and abs(y[i] - width_y / 2) <= 0.2 * width_y\n",
    "            ):\n",
    "                D_matrix[i, j] = Diff_low\n",
    "                # here we define a square of low diffusivity in the middle\n",
    "    return x, y, D_matrix, c_init\n",
    "\n",
    "\n",
    "width_x = 10  # dm\n",
    "width_y = 10  # dm\n",
    "n_x = 21\n",
    "n_y = n_x\n",
    "poro = 0.4\n",
    "\n",
    "the_prob = Problem_Def(n_x, n_y, poro, width_x, width_y)\n",
    "x, y, D_matrix, c_init = make_D_matrix(the_prob, Diff, decreasing_factor)\n",
    "fig, ax = plt.subplots()\n",
    "# This generates a colormap of diffusion.\n",
    "cm = cmap.get_cmap(\"magma\")\n",
    "plt.contourf(x, y, D_matrix, cmap=cm)\n",
    "plt.colorbar()\n",
    "\n",
    "# \"magma\" refers to a colormap example. You can chose other ones\n",
    "# https://matplotlib.org/examples/color/colormaps_reference.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "lines_to_next_cell": 2,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-3f0b9e511535468c",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x115d89278>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAAD8CAYAAABErA6HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFJlJREFUeJzt3X+s3fV93/HnK3YM4UcaUydZazvBqE4WQjZIPZIGDbGQH05T4WwJqZnoSJTNqgRJmrarYNpgIorGqippq7khV4mbdE1xEaBgVVYpMmGsy6A2hF82dTEmhRvTEsrvjkBs3vvjfE0Ol3vvOb7nxvfje58P6eh8P5/v5/P9fI4ML3/9Od/v96SqkCS141VzPQFJ0ssZzJLUGINZkhpjMEtSYwxmSWqMwSxJjRkYzEk2JXk0yb19dSckuTHJ/d370p/sNCVp7kyWgxP2J8nvJ9mT5O4k7+zbdyDJnd1ryzDjDXPG/HVg7YS6i4FtVbUa2NaVJWm++jqvzMF+HwJWd68NwJf79j1XVad2r3OGGWxgMFfVLcDjE6rXAd/otr8BfGSYwSTpSDRFDvZbB/xR9dwKvC7Jz8x0vMUz7PfGqnoEoKoeSfKGqRom2UDvbxCOIj+/+vU/BcCSE07o7T+utwry3HPPAfD0008D8OSTT85wapLmg2efffaxqnr9TPuf8upj69k6MFTbvz3w/E7gh31VY1U1dgjDLQce7iuPd3WPAEcn2QHsB66oqm8NOthMg3lo3YcbAzhx8dF13UfPBuDN538cgCXvOReAu+6+G4Bt27YBcP311/+kpyapYbfccsvfjtL/2TrAf3ntiUO1/fdP7P5hVa0ZYbhMUnfweRdvqqp9SU4CbkpyT1U9MN3BZnpVxt8fPE3v3h+d4XEkaT4YB1b2lVcA+wCq6uD7XuBm4LRBB5tpMG8BLui2LwA8vZW0kG0B/l13dca7gae6Zd6lSY4CSLIMOAPYNehgA5cyklwFnAUsSzIOXAZcAVyd5FPAQ8C5M/00ktS6KXLw1QBVdSWwFfhFYA/w/4BPdl3fBnwlyYv0ToSvqKrRg7mqzpti19mD+krSfDBNDh7cX8CFk9R/B3jHoY7nnX+S1BiDWZIaYzBLUmMMZklqjMEsSY0xmCWpMQazJDXGYJakxhjMktQYg1mSGmMwS1JjDGZJaozBLEmNMZglqTEGsyQ1xmCWpMYYzJLUGINZkhpjMEtSYwxmSWqMwSxJjTGYJakxBrMkDSHJ2iS7k+xJcvEk+9+cZFuSu5PcnGRF374LktzfvS4YNJbBLEkDJFkEbAQ+BJwMnJfk5AnNfgf4o6r6Z8DlwH/r+p4AXAa8CzgduCzJ0unGM5glabDTgT1VtbeqXgA2A+smtDkZ2NZtf7tv/weBG6vq8ap6ArgRWDvdYAazJMGyJDv6Xhsm7F8OPNxXHu/q+t0FfLTb/tfA8Ul+esi+L7P4UGcvSUeC4094DWd+9JThGl+5+7GqWjNNi0xSVxPKvwn8jySfAG4Bvg/sH7LvyxjMkjTYOLCyr7wC2NffoKr2Af8GIMlxwEer6qkk48BZE/rePN1gLmVI0mDbgdVJViVZAqwHtvQ3SLIsycFMvQTY1G3fAHwgydLuS78PdHVTMpglaYCq2g9cRC9Q7wOurqqdSS5Pck7X7Cxgd5K/Ad4IfKHr+zjweXrhvh24vKubkksZkjSEqtoKbJ1Qd2nf9jXANVP03cSPz6AH8oxZkhpjMEtSYwxmSWqMwSxJjRkpmJN8LsnOJPcmuSrJ0bM1MUlaqGYczEmWA58B1lTVKcAietf2SZJGMOpSxmLgNUkWA8cw4U4YSdKhm3EwV9X36T3m7iHgEeCpqvqLie2SbDj4YJBn6sDMZypJC8QoSxlL6T3WbhXws8CxSc6f2K6qxqpqTVWtOT6LZj5TSVogRlnKeB/wYFX9oKp+BFwHvGd2piVJC9cowfwQ8O4kxyQJcDa9e8glSSMYZY35Nnr3hd8B3NMda2yW5iVJC9ZIDzGqqsvo/ZaVJGmWeOefJDXGYJakxhjMktQYg1mSGmMwS1JjDGZJaozBLEmNMZglqTEGsyQ1xmCWpCEkWZtkd5I9SS6eZP+XktzZvf4myZN9+w707dsyaKyRbsmWpIUgySJgI/B+YBzYnmRLVe062KaqPtfX/tPAaX2HeK6qTh12PM+YJWmw04E9VbW3ql4ANtN7Hv1UzgOumulgBrMkwbKDv7TUvTZM2L8ceLivPN7VvUKSN9P7AZGb+qqP7o57a5KPDJqMSxmS5qUlJ5zAm8//+HCNr7z2sapaM02LTFJXU7RdD1xT9bLf0ntTVe1LchJwU5J7quqBqQbzjFmSBhsHVvaVVzD1j0+vZ8IyRlXt6973Ajfz8vXnVzCYJWmw7cDqJKuSLKEXvq+4uiLJW4GlwP/tq1ua5KhuexlwBrBrYt9+LmVI0gBVtT/JRcANwCJgU1XtTHI5sKOqDob0ecDmqupf5ngb8JUkL9I7Gb6i/2qOyRjMkjSEqtoKbJ1Qd+mE8n+dpN93gHccylguZUhSYwxmSWqMwSxJjTGYJakxBrMkNcZglqTGGMyS1BiDWZIaYzBLUmMMZklqjMEsSY0xmCWpMQazJDXGYJakxhjMktQYg1mSGjNSMCd5XZJrkvx1kvuS/MJsTUySFqpRf8Hk94A/r6qPdb+DdcwszEmSFrQZB3OS1wJnAp8AqKoXgBdmZ1qStHCNspRxEvAD4A+TfDfJV5McO7FRkg1JdiTZ8UwdGGE4SVoYRgnmxcA7gS9X1WnAPwIXT2xUVWNVtaaq1hyfRSMMJ0kLwyjBPA6MV9VtXfkaekEtSRrBjIO5qv4OeDjJW7uqs4FdszIrSWpMkrVJdifZk+QVqwNdm48n2ZVkZ5I/6au/IMn93euCQWONelXGp4Fvdldk7AU+OeLxJKk5SRYBG4H301st2J5kS1Xt6muzGrgEOKOqnkjyhq7+BOAyYA1QwO1d3yemGm+kYK6qO7vBJGk+Ox3YU1V7AZJsBtbx8lWC/wBsPBi4VfVoV/9B4MaqerzreyOwFrhqqsG880+SBlsOPNxXHu/q+r0FeEuS/5Pk1iRrD6Hvy4y6lCFJTcpxS1nynnOHbP3Ly5Ls6KsYq6qx/sNN0qkmlBcDq4GzgBXA/05yypB9X3EgSVroHquq6ZZlx4GVfeUVwL5J2txaVT8CHkyym15Qj9ML6/6+N083GZcyJGmw7cDqJKu6ix3WA1smtPkW8K8Akiyjt7SxF7gB+ECSpUmWAh/o6qbkGbMkDVBV+5NcRC9QFwGbqmpnksuBHVW1hR8H8C7gAPAfq+ofAJJ8nl64A1x+8IvAqRjMkjSEqtoKbJ1Qd2nfdgG/3r0m9t0EbBp2LJcyJKkxBrMkNcZglqTGGMyS1BiDWZIaYzBLUmMMZklqjMEsSY0xmCWpMQazJDXGYJakxhjMktQYg1mSGmMwS1JjDGZJaozBLEmNMZglqTEGsyQ1xmCWpMYYzJLUGINZkhpjMEtSYwxmSRpCkrVJdifZk+Tiadp9LEklWdOVT0zyXJI7u9eVg8ZaPJsTl6T5KMkiYCPwfmAc2J5kS1XtmtDueOAzwG0TDvFAVZ067HieMUvSYKcDe6pqb1W9AGwG1k3S7vPAbwM/HGUwg1mSBlsOPNxXHu/qXpLkNGBlVf3ZJP1XJflukv+V5F8OGsylDEnz0nPPPcddd989bPNlSXb0lceqaqyvnEn61Es7k1cBXwI+MUm7R4A3VdU/JPl54FtJ3l5VT081GYNZkuCxqlozzf5xYGVfeQWwr698PHAKcHMSgH8CbElyTlXtAJ4HqKrbkzwAvAXo/4vgZVzKkKTBtgOrk6xKsgRYD2w5uLOqnqqqZVV1YlWdCNwKnFNVO5K8vvvykCQnAauBvdMNNnIwJ1nUrZ1Mtq4iSUe8qtoPXATcANwHXF1VO5NcnuScAd3PBO5OchdwDfCrVfX4dB1mYynjs91EXzsLx5KkJlXVVmDrhLpLp2h7Vt/2tcC1hzLWSGfMSVYAHwa+OspxJEk/NupSxu8CvwW8OFWDJBuS7Eiy45k6MOJwkjT/zTiYk/wS8GhV3T5du6oaq6o1VbXm+N76tyRpGqOcMZ8BnJPke/Tugnlvkj+elVlJ0gI242CuqkuqakV3ach64KaqOn/WZiZJC5TXMUtSY2blzr+quhm4eTaOJUkLnWfMktQYg1mSGmMwS1JjDGZJaozBLEmNMZglqTEGsyQ1xmCWpMYYzJLUGINZkhpjMEtSYwxmSWqMwSxJjTGYJakxBrMkNcZglqTGGMySNIQka5PsTrInycWT7P/VJPckuTPJXyY5uW/fJV2/3Uk+OGgsg1mSBkiyCNgIfAg4GTivP3g7f1JV76iqU4HfBr7Y9T2Z3u+ivh1YC/xBd7wpGcySNNjpwJ6q2ltVLwCbgXX9Darq6b7isUB12+uAzVX1fFU9COzpjjelWfnNP0lqzdNPP822bduGbb4syY6+8lhVjfWVlwMP95XHgXdNPEiSC4FfB5YA7+3re+uEvsunm4zBLEnwWFWtmWZ/JqmrV1RUbQQ2Jvm3wH8GLhi2bz+XMiRpsHFgZV95BbBvmvabgY/MsK/BLElD2A6sTrIqyRJ6X+Zt6W+QZHVf8cPA/d32FmB9kqOSrAJWA3813WAuZUjSAFW1P8lFwA3AImBTVe1Mcjmwo6q2ABcleR/wI+AJessYdO2uBnYB+4ELq+rAdOMZzJI0hKraCmydUHdp3/Znp+n7BeALw47lUoYkNcZglqTGGMyS1BiDWZIaYzBLUmMMZklqjMEsSY0xmCWpMQazJDVmxsGcZGWSbye5L8nOJFPe9SJJGt4ot2TvB36jqu5Icjxwe5Ibq2rXLM1NkhakGZ8xV9UjVXVHt/0McB8DHv4sSRpsVh5ilORE4DTgtkn2bQA2AJzwKp+ZJEmDjPzlX5LjgGuBX5vwm1cAVNVYVa2pqjXHT//7g5IkRgzmJK+mF8rfrKrrZmdKkrSwjXJVRoCvAfdV1Rdnb0qStLCNcsZ8BvArwHuT3Nm9fnGW5iVJC9aMv42rqr9k8l9/lSSNwDv/JKkxBrMkNcZglqTGGMyS1BiDWZKGkGRtkt1J9iS5eJL9Zya5I8n+JB+bsO9A39VrWwaN5T3SkjRAkkXARuD9wDiwPcmWCQ9tewj4BPCbkxziuao6ddjxDGZJGux0YE9V7QVIshlYB7wUzFX1vW7fi6MOZjBLmpeefPJJrr/++mGbL0uyo688VlVjfeXlwMN95XHgXYcwnaO74+8Hrqiqb03X2GCWJHisqtZMs3+ym+nqEI7/pqral+Qk4KYk91TVA1M19ss/SRpsHFjZV14B7Bu2c1Xt6973AjfTe0zylAxmSRpsO7A6yaokS4D1wMCrKwCSLE1yVLe9jN5zhqb9pSeDWZIGqKr9wEXADfR+renqqtqZ5PIk5wAk+RdJxoFzga8k2dl1fxuwI8ldwLfprTFPG8yuMUvSEKpqK7B1Qt2lfdvb6S1xTOz3HeAdhzKWZ8yS1BiDWZIaYzBLUmMMZklqjMEsSY0xmCWpMQazJDXGYJakxhjMktQYg1mSGmMwS1JjDGZJaozBLEmNMZglqTEGsyQ1xmCWpMYYzJLUGINZkhpjMEtSYwxmSWqMwSxJjTGYJakxIwVzkrVJdifZk+Ti2ZqUJLVmUN4lOSrJn3b7b0tyYt++S7r63Uk+OGisGQdzkkXARuBDwMnAeUlOnunxJKlVQ+bdp4AnqurngC8B/73rezKwHng7sBb4g+54UxrljPl0YE9V7a2qF4DNwLoRjidJrRom79YB3+i2rwHOTpKufnNVPV9VDwJ7uuNNafEIE10OPNxXHgfeNbFRkg3Ahq74/FuuvPZeAK68tqv65RGm0IRlwGNzPYmfAD/XkWM+fiaAt47S+dlnn73hlltuWTZk86OT7Ogrj1XVWF95mLx7qU1V7U/yFPDTXf2tE/oun24yowRzJqmrV1T0PtwYQJIdVbVmhDGbMx8/E/i5jiTz8TNB73ON0r+q1s7WXBgu76ZqM1RW9htlKWMcWNlXXgHsG+F4ktSqYfLupTZJFgM/BTw+ZN+XGSWYtwOrk6xKsoTe4vaWEY4nSa0aJu+2ABd02x8Dbqqq6urXd1dtrAJWA3813WAzXsro1lAuAm4AFgGbqmrngG5jA/YfiebjZwI/15FkPn4maOhzTZV3SS4HdlTVFuBrwP9MsofemfL6ru/OJFcDu4D9wIVVdWC68dILdElSK7zzT5IaYzBLUmMOSzDPx1u3k6xM8u0k9yXZmeSzcz2n2ZJkUZLvJvmzuZ7LbEnyuiTXJPnr7s/sF+Z6TrMhyee6//7uTXJVkqPnek6HKsmmJI8mubev7oQkNya5v3tfOpdzPNx+4sE8j2/d3g/8RlW9DXg3cOE8+VwAnwXum+tJzLLfA/68qv4p8M+ZB58vyXLgM8CaqjqF3pdS6+d2VjPydXq3Kve7GNhWVauBbV15wTgcZ8zz8tbtqnqkqu7otp+h9z/6tHfzHAmSrAA+DHx1rucyW5K8FjiT3rfmVNULVfXk3M5q1iwGXtNdN3sMR+C9BFV1C72rGPr13978DeAjh3VSc+xwBPNktzIe8QHWr3uK1GnAbXM7k1nxu8BvAS/O9URm0UnAD4A/7JZovprk2Lme1Kiq6vvA7wAPAY8AT1XVX8ztrGbNG6vqEeidBAFvmOP5HFaHI5gP+XbEI0mS44BrgV+rqqfnej6jSPJLwKNVdftcz2WWLQbeCXy5qk4D/pF58E/jbt11HbAK+Fng2CTnz+2sNBsORzDP21u3k7yaXih/s6qum+v5zIIzgHOSfI/ektN7k/zx3E5pVowD41V18F8019AL6iPd+4AHq+oHVfUj4DrgPXM8p9ny90l+BqB7f3SO53NYHY5gnpe3bneP8/sacF9VfXGu5zMbquqSqlpRVSfS+3O6qaqO+DOwqvo74OEkB59Wdja9u7COdA8B705yTPff49nMgy81O/23N18AXD+HcznsRnm63FBmeOv2keAM4FeAe5Lc2dX9p6raOodz0tQ+DXyzOznYC3xyjuczsqq6Lck1wB30rhL6Lg3dxjysJFcBZwHLkowDlwFXAFcn+RS9v4DOnbsZHn7eki1JjfHOP0lqjMEsSY0xmCWpMQazJDXGYJakxhjMktQYg1mSGvP/AYcJJqWkOvfgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here we plot the initial condition using the colormap again\n",
    "fig, ax = plt.subplots()\n",
    "# This generates a colormap of diffusion.\n",
    "cm = cmap.get_cmap(\"RdGy_r\")\n",
    "plt.contourf(x, y, c_init, cmap=cm)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Steady state the old-fashioned way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-bbb5e477a2fcd825",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x11639cb38>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAAD8CAYAAABErA6HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHS1JREFUeJzt3X+s3fV93/Hnyz+pIYl/BQuundoEQ0PIFlqP0KJFLJTEWSLI1nQyVSISRbMqQZKm7SqYNpgcRcmqKD+kudushLZbUzxEosSKrBDEj3VdG2qHkASbGC7Gg2uIARuTdjBfX+e9P87X5Phw7z3fc77f7/l+vt/zekhHPuf783Owed/3fX9+fBURmJlZOhbU3QAzMzuTA7OZWWIcmM3MEuPAbGaWGAdmM7PEODCbmSWmb2CWdLuk5yQ90rVtpaR7JD2e/bmi2maamY2PPBnznwGbe7bdDNwbERuBe7PPZmZWAuWZYCJpPfDtiLg0+3wAuCoinpV0HvBARFxcZUPNzMbFoiHPWxMRzwJkwfncuQ6UtBXYCrAY/doqLR7ylnDWguFL4osXD3/uoqULhzpv4dLBv+sw5yxYumSg47V46cD30KIB7jHIsUBosL+bmZmZgY4/depUpdcf9HiA6enpgc85efLkwOcMe68yzh22vaedOHHihYh447DnX7BgWbxCvr/7n8b03RHRWxmozbCBObeI2AHsADhvwdL46OK1Q1/ronMG+x++28S5y4Y+d+WFK4c6b8WFg/+besObJwY+Z9n69QMdv/i8wY4HWLAm/zlatW6ga5/QYH+vx48fr/T4o0ePDnT8kSNHBjoe4Jlnnhn4nKmpqYHPATh8+PBQ5xW5Jwz3Hbs99thj/6fI+a9wirzx5rPTB1cXuVfZhk0jj2QlDLI/nyuvSWZm423YwLwLuCF7fwPwrXKaY2ZmeYbL3QH8LXCxpClJHwM+B1wj6XHgmuyzmZmVoG+NOSKun2PX1SW3xczM8Mw/M7PkODCbmSXGgdnMLDEOzGZmiXFgNjNLjAOzmVliKp+SbYN56YnDQ03LHsTJZw8NNS3brEnOWrAg/zIOx6pty6CcMZuZJWZsAvPh516uuwlmZrmMTWA2M2sKB2Yzs8Q4MOdwbHK4noEXJ58vuSVmNg4cmM3MEuPA3AIvHzpUdxPOEEefrrsJZqWTtFnSAUmTkl7zAGpJvyzpXkk/kvSApLVd+26Q9Hj2uqH33F4OzGZmfUhaCGwH3gtcAlwv6ZKewz4P/LeI+EfANuCz2bkrgduAdwCXA7dJWjHf/RyYzUas6LPwrBaXA5MRcTAipoGdwHU9x1wC3Ju9v79r/3uAeyLiWES8CNwDzPvgVwdmMzNYLWlv12trz/4JoLtGN5Vt6/ZD4Ley9/8CeJ2kVTnPPYOnZJvZrIo8ITsFixcvYOLcZfkOPsYLEbFpniM0y7bo+fyHwH+S9BHgr4DDwEzOc8/gwGzWAMMGycOHD5fckrE1Bazr+rwWOKMmFRHPAP8SQNI5wG9FxEuSpoCres59YL6buZRh1uXo0aN1N8HStAfYKGmDpCXAFmBX9wGSVks6HVNvAW7P3t8NvFvSiqzT793Ztjk1KjA/9g/TdTdhJF56wlmOWUoiYga4iU5AfRS4MyL2Sdom6drssKuAA5IeA9YAn8nOPQZ8mk5w3wNsy7bNyaWMnI5NHmPlhSsHPu/FyedZceEbK2iRpeDIkSN1N8FGJCJ2A7t7tt3a9f4u4K45zr2dX2TQfTUqYy7KK8yZWROMVWA2M2sCB2ZrtePHj9fdBLOBOTCbjdAws/6aPp7YBufA3BKpLWRk9fMY5ubyqAwza6VFSxfmH0n1k2rbMihnzAMY5YL5VY9lPvnsoUqvb83m8km9HJjNrHReQa+YsQvMHstsTeLMdTyNXWA2K8ugs/5GmUW646/ZHJgH5AezmlnVHJjN7Awun9SvcYF5XFaYA68yN+4cIMdXocAs6VOS9kl6RNIdks4qq2E2OE8yKaYtazG7vtx8QwdmSRPAJ4BNEXEpsJDO4tHJq2tkRpPrzD8/cqjuJpiNjaKljEXAL0laBCyj51ErbTVsB6CNL4/rtUEMPSU7Ig5L+jzwFPAK8N2I+G7vcdnTZrcCvN4zwM1ycX25uIVLFzf2IRVFShkrgOuADcD5wNmSPtR7XETsiIhNEbFpmRrX1zinUQ2bcwdgmlJ9cknR+nIZPxD820FxRSLlbwJPRsTzEXES+AbwG+U0q3qeAWhmqSoSmJ8CrpC0TJKAq+k8pND6SKUT0AsZpWnYrNWjMdpj6MAcEQ/SefDgQ8CPs2vtKKldjeBOwNFr4hNJ/Ku9DapQ0TcibouIX4mISyPiwxFxoqyGzacNk0yqypo9lrnZ6uz0c4djOtrTGzeEMurMo8ia3QFo/biM0S5jHZjNUuKM1U5zYK5RKp2AVp1R1JfLyJb9Q6E/SZslHZA0KenmWfZ/UdLD2esxSce79p3q2rer373GfsbH4edeZuLcZYWucWzyWP5ni9nIVNVRmOoY5hS0taNT0kJgO3ANMAXskbQrIvafPiYiPtV1/MeBy7ou8UpEvD3v/ZwxN4TrzNWqewEjZ6zJuxyYjIiDETEN7KQzwW4u1wN3DHuzsc+YyzJs1vzi5POlTxt9+dAhlq1fX+o1LU3u9JvbwqWLecObJ/IevlrS3q7POyKie/jvBPB01+cp4B2zXUjSL9OZEX1f1+azsuvPAJ+LiG/O15jGZsxlDpkraxZg1dO0q8iaPcmkOnl/ra97QomzdQBeOL10RPbqnZOhWc6JOa61BbgrIk51bXtTRGwCfgf4kqQ3z9eYxgbmVDUxOOflpT/z15errrU6Ux65KWBd1+e1zL2a5hZ6yhgR8Uz250HgAc6sP7+GA3MFUpgR6IkmzTBMtlpmUHa2nNseYKOkDZKW0Am+rxldIeliYAXwt13bVkhamr1fDVwJ7O89t5sDc6bsRY2GCc5NyJrbqOyOvypLGClnym0dkQEQETPATcDddNYEujMi9knaJunarkOvB3ZGRHeZ4y3AXkk/BO6nU2OeNzC7869Cw3QIltkZmLcT8OSzh1h8Xv/jBhFHn0ar1vU/sEHKHCbXtqA8DiJiN7C7Z9utPZ//wyzn/Q3wtkHu5Yy5SxVLgVZV1nDWPL+6FjuqKmusIii7jJGuRgfmpixmNGhw9ozA9OTJlqsqYThTHj+NDsxVqGoB/SqCc56s2Z2A8xv1xJJUgrKz5bQ5MI9QXcE5j7zjmT1kbnZVlDCcKY8vB+ZZVPnYqRSG0tlgyur0GyRLrTIoV5Ett3lERh08KqMGg4zWyDNK46UnDs879bQtU7TL7tArq4yRJyilEpTHyYKlSxr7777xGXNVHYBVP6zVmXN1mvj4qdOqDsquLTdD4wPzOChjlEaeTkCvm/Fa/coYZWfLZuDAPK+UsuZ+wXmU45rdAZhfSiUM/4BoDgdmG0t56steED8fd/yVz4G5jzZlzaMuZ8TRp/sf1GD9ApKzZRtWKwJz1TMAqw7OVp4mdvw5KFuvVgTmpisza7ZyFO30czC0IhyYc2pK1lxGOSOPJncA1v18v27Olm02DsyJSClr9rC5YsYpGLrjrxoOzANoStZsxZQxdjkPZ8s2l9ZMyX7sH6a56JwldTejkGGftN2rjVO0m9ipZ/XS4qWlPwBiVJwxD8hZc7ryBO+q68upZKmptMOG48A8hBSCs0dnNFvTVo+z0XJgTkwqixuNawfgqOrLbeD/FtVxYB5SClmzWS9ny+3QqsDclGcAjsKoxjPbL+QNil5v2fppVWAetVSeD1iXJk8yaSNny+1RKDBLWi7pLkk/kfSopF8vq2HWnzsAB5PSjD9rHkmbJR2QNCnp5jmO+VeS9kvaJ+kvu7bfIOnx7HVDv3sVHcf8ZeA7EfFBSUuAZQWv1ziHn3uZiXPH7mu30ig6/try1Otx6/iTtBDYDlwDTAF7JO2KiP1dx2wEbgGujIgXJZ2bbV8J3AZsAgL4fnbui3Pdb+iMWdLrgXcCXwWIiOmIGMtZAFWUNMooZ4xy8fy5tH3pzxS4hDESlwOTEXEwIqaBncB1Pcf8a2D76YAbEc9l298D3BMRx7J99wCb57tZkVLGBcDzwJ9K+oGkr0g6u/cgSVsl7ZW09+X4eYHb5eMOwPz6dQCO65C5YTg4Nt7q03Eqe23t2T8BdGcZU9m2bhcBF0n635K+J2nzAOeeoUgpYxHwq8DHI+JBSV8Gbgb+ffdBEbED2AFw3oKlUeB+SXNJo15NmLJdRRnDPxDmpkVLWLBmfd7DX4iITfNdbpZtvfFsEbARuApYC/wvSZfmPPcMRTLmKWAqIh7MPt9FJ1CPLY9ttlFyUB6pKWBd1+e1QG+hfQr4VkScjIgngQN0AnWec88wdGCOiJ8CT0u6ONt0NbB/nlNsQHnqzB6ZYTYSe4CNkjZkAx22ALt6jvkm8M8AJK2mU9o4CNwNvFvSCkkrgHdn2+ZUdFTGx4GvZQ09CHy04PUazyWNdkptFEKd2XJq/y1GISJmJN1EJ6AuBG6PiH2StgF7I2IXvwjA+4FTwL+JiKMAkj5NJ7gDbIuIebOuQoE5Ih6mMwTEzKzVImI3sLtn261d7wP4/ezVe+7twO157+WZfxVwrTk9/SaX9BvDXFSZHX+uLbdfKwOzh8zZKDlQWtlaGZhTUFbWXHSiSQqTTKw8/iEwHhyYrRAvZGRWPgfmCjWh1uzlP5vD2fL4cGA2q1ib1l8ex6FydWjNU7LH2YuTz7PiwjfW3QyrkLPlISxaglat639cgpwxV6wJ5QwzS4sDs83LK8yZjZ4Ds5lZYhyYzcwS48Bs1kfdIxHc8Td+HJjNzBLT2sCc0noZRUdmlPH8PzNrjtYGZjOzpnJgNjNLjAOzmVliPCXbzFoptIATWlJ3M4bijNnMLDGtDcwXndPMn5SzWXnhyrqbYGYj1NrAnBI/NdvMBuHAbNbH+eefX3cTbMw4MJslbu3atXU3wUbMgdnMLDEOzGZmOUjaLOmApElJN89z3AclhaRN2ef1kl6R9HD2+i/97uVxzBVresff4vPW190Es9pJWghsB64BpoA9knZFxP6e414HfAJ4sOcST0TE2/PezxlzC/h5f2aVuxyYjIiDETEN7ASum+W4TwN/DPy/IjdzYDar2MTEROFrpNIB2KQRKjMzMxw/fjzXC1gtaW/Xa2vP5SaAp7s+T2XbXiXpMmBdRHx7luZskPQDSf9T0j/t13aXMirUhDLGsvXr626C5bR27Vovml+dFyJi0zz7Ncu2eHWntAD4IvCRWY57FnhTRByV9GvANyW9NSJ+NtfNnDFbIQvWrK+7CWajMAWs6/q8Fuh+tM3rgEuBByQdAq4AdknaFBEnIuIoQER8H3gCuGi+mzkwV6SsbLnodOw3vLn4r9GWjlRKGmNoD7BR0gZJS4AtwK7TOyPipYhYHRHrI2I98D3g2ojYK+mNWechki4ANgIH57tZKwNzm9bJsPQ5WLZfRMwANwF3A48Cd0bEPknbJF3b5/R3Aj+S9EPgLuB3I2LexxK5xlyBJtSWx82qVas4evRobfefmJjg8OHDpVzLteZ6RMRuYHfPtlvnOPaqrvdfB74+yL1amTGbDWrNmjV1N8HsVQ7MJXO23E6pDROrs3yS2n+LNiocmCUtzMbnzTZ2zwrwOsxm46mMjPmTdIrhY6+ObNmz/saXOx3bq1BglrQWeB/wlXKaY9ZeZcwAtPFQdFTGl4A/ojO4elbZ1MatAK8fwSCQuobKNbG23G/WX5MWMFq+fPnpqbVjxSM05nbq1KnG/psYOmOW9H7guWwmy5wiYkdEbIqITcvUzr7GKoJyGfXlFCaXaNW6/ge1QJ1lBZc02qdIpLwSuDabfrgTeJekvyilVQ3SxEzZZjeKIXNtKWd4ZEa1hg7MEXFLRKzNph9uAe6LiA+V1jLryx1/g1m1alWh81MORs6a26WdtYURqSpb9jA5G4aDc3uUEpgj4oGIeH8Z1yrCa2T8Qr/6chnLfXpluTPlDYxtKWdYdZwxD8m1ZUuRs+Z2cGBOTCpljCYNlStTvw7AlOvMo+b/FtVxYB5CCtmyO/6arcpyhrPm5nNgHlAKQdlmt3z58r7HFB2ZYTYKrQnMbej4K6uMMYqOv1HLE3RTkEq2mko7bDheKH8AzpbHw5o1azhy5Mic+88//3yeeeaZOffnVebi+bMZ9+naMzMztT4coYjWZMxNN0i2XHV9eVw7/soyTtmqOwCr4cCcU1Oy5VGVMZo8hjmlOnPVY5rH6YdEmzgwJyClbNk6ig6bc0C0IloRmKvu+GtKtmzN6STs5qzZerUiMFep6qBcZrZcRhmjzPpy25f8dNZsVXFgtrGUp848yidnNzlrdgdg+RyY59GmbLlMTe74G7WUsuaU2tJEkjZLOiBpUtLNs+z/XUk/lvSwpL+WdEnXvluy8w5Iek+/ezkwj4kmTipJwSjXzvCqc+mStBDYDrwXuAS4vjvwZv4yIt4WEW8H/hj4QnbuJXTWrH8rsBn4k+x6c2p8YK6q469J2XJZ2jJ+ObUOwEEy1SaXNFrucmAyIg5GxDSdpzZd131ARPys6+PZQGTvrwN2RsSJiHgSmMyuN6fGB+a2yxOU2zgFezZlB9yyxjPnyZpTCs42lAng6a7PU9m2M0i6UdITdDLmTwxybjdPyZ5FldlyKst6Wn79pmhXocrp2lVM1S5rmnqZZmZmBvl7Wy1pb9fnHRGxo+uzZjknXrMhYjuwXdLvAP8OuCHvud2cMY/QoEG5jGw5r7aUMepSdtYMzpxH7IWI2NT12tGzfwroHv+5FpjvJ9FO4ANDnuvA3CuV5/iVVVcuu4zhERnFpBKcXWse2B5go6QNkpbQ6czb1X2ApI1dH98HPJ693wVskbRU0gZgI/B3893MgXkEqipfjHKIXFvlrTPnGdOcd4RGKsHZ8ouIGeAm4G7gUeDOiNgnaZuka7PDbpK0T9LDwO/TKWMQEfuAO4H9wHeAGyPi1Hz3a3SNuewRGVVky8ME5VGWMNpq+fLlHD9+fOT3rarWWkXNedyXBR1UROwGdvdsu7Xr/SfnOfczwGfy3ssZc4WqCsp55S1juL6cT5kzAYcpJThzHh8OzJmys+Uqg3ITsuW2r5Mxn6pKGpB2cPbU7PI4MFcghSFxbRm7PAoprc+cR5nB2Z2AaXJgLtmwQbkJ2bJHZOQvZ1SZNUPambMV19jAXGbHX1lljKqDchVcX65OU4Kzs+b0NDYwp2YUQTlvtuwyxvhw5txOjR4uN06a0OHXZKtWrar1icoeula+6enp5KaJ5zX2GXMZZYwUOvvstapaZW6UC+g3jUdmlGPsA3Od/GDV9htFoCqjnOE6c1ocmM0S4eBop411YHYZw9rCnYDt0sjAXNVTS0apitEY4BEZTees2aChgTkVzpYtj6Z0iPmHQjqGDsyS1km6X9Kj2VJ3c66sZGdKpdOviZNLUnueX0pczmiPIhnzDPAHEfEW4ArgxlmeGpusqh+2ajYsZ642dGCOiGcj4qHs/d/TWTx6bH5kj2r6tSeWpCnVscxFs+Yyfig0pXSTslJqzJLWA5cBD86yb6ukvZL2vhw/L+N2Zq3nrHm8FZ6SLekc4OvA70XEz3r3Zw813AFw3oKl8z4Z1qytUnyKdNudPHmysdPcC2XMkhbTCcpfi4hvlNOk6hWtLzdxFbmivOSn2egUGZUh4KvAoxHxhfKaZMPyGOZimrZg/lw8OqP5imTMVwIfBt4l6eHs9c9Latec2jC5xCwP15nH19A15oj4a0AltsV6eESG1cFLkNbPM/8G5Nl+ZlY1B+YRaXLHn81u0LHMoxzf6zpzs41dYPaMP2sS15nTIWmzpAOSJiXdPMv+d0p6SNKMpA/27DvV1Re3q9+9/GgpM7M+JC0EtgPXAFPAHkm7ImJ/12FPAR8B/nCWS7wSEW/Pe7+xy5iLGGV9ueqOvyYuYGTN0cJp2ZcDkxFxMCKmgZ3Add0HRMShiPgRUHiKswPzCLi+bE0zhiWU1aeXjsheW3v2TwBPd32eYrC1gc7Krvs9SR/od7BLGS3hySXWa2JigsOHD9fdjNpMT08P8v1fiIhN8+yfbWjwIEtMvCkinpF0AXCfpB9HxBNzHeyM2WyEhvkVfwyz1xRNAeu6Pq8Fci9+EhHPZH8eBB6gs+jbnByYzcz62wNslLRB0hJgC9B3dAWApBWSlmbvV9OZNb1/vnMcmK3V/MQTK0NEzAA3AXfTWXv+zojYJ2mbpGsBJP0TSVPAbwP/VdK+7PS3AHsl/RC4H/hcz2iO1xirGrPHMJvZsCJiN7C7Z9utXe/30Clx9J73N8DbBrmXM+acxnGpT+sv1SeZWLM1KjCPy8pyXrzIbLw1KjCbVa0tazJbszkwm7WYFzNqJgdmswaoYyyzx0/Xx4HZzCwxYzVcziwFfmL2aExPTzf2SSzOmM3MEuPA3AKpLWCkVev6H2Rmc3JgNjNLjANzDp71Z2aj5MBsZpYYB2Yzs8SMTWD2ynJm1hRjE5jNzJrCgTkxo1hZzk/INkubZ/6ZWSudPHmysTMsnTGbmSXGgdnMLDEOzGZmiXFgNjNLjAOzmVliHJjNzBLjwGxmlphCgVnSZkkHJE1KurmsRpmZpaZfvJO0VNL/yPY/KGl9175bsu0HJL2n372GDsySFgLbgfcClwDXS7pk2OuZmaUqZ7z7GPBiRFwIfBH4j9m5lwBbgLcCm4E/ya43pyIZ8+XAZEQcjIhpYCdwXYHrmZmlKk+8uw748+z9XcDVkpRt3xkRJyLiSWAyu96cikzJngCe7vo8Bbyj9yBJW4Gt2ccTn50++MjQdzw29JnFzv3JvHtXAy8UuHqq/L2ao43fCeDiIiefOHHi7scee2x1zsPPkrS36/OOiNjR9TlPvHv1mIiYkfQSsCrb/r2ec+ddFKdIYNYs2+I1GzpfbgeApL0RsanAPZPTxu8E/l5N0sbvBJ3vVeT8iNhcVlvIF+/mOiZXrOxWpJQxBXQ/dXMt0MwVQ8zM5pcn3r16jKRFwBvo/K4+cKwsEpj3ABslbZC0hE5xe1eB65mZpSpPvNsF3JC9/yBwX0REtn1LNmpjA7AR+Lv5bjZ0KSOrodwE3A0sBG6PiH19TtvRZ38TtfE7gb9Xk7TxO0FC32uueCdpG7A3InYBXwX+u6RJOpnyluzcfZLuBPYDM8CNEXFqvvupE9DNzCwVnvlnZpYYB2Yzs8SMJDC3ceq2pHWS7pf0qKR9kj5Zd5vKImmhpB9I+nbdbSmLpOWS7pL0k+zv7NfrblMZJH0q+/f3iKQ7JJ1Vd5sGJel2Sc9JeqRr20pJ90h6PPtzRZ1tHLXKA3OLp27PAH8QEW8BrgBubMn3Avgk8GjdjSjZl4HvRMSvAP+YFnw/SRPAJ4BNEXEpnU6pLfW2aih/RmeqcrebgXsjYiNwb/Z5bIwiY27l1O2IeDYiHsre/z2d/9Grf8R1xSStBd4HfKXutpRF0uuBd9LpNScipiPieL2tKs0i4JeycbPLaOBcgoj4K147N7d7evOfAx8YaaNqNorAPNtUxsYHsG7ZKlKXAQ/W25JSfAn4I+DndTekRBcAzwN/mpVoviLp7LobVVREHAY+DzwFPAu8FBHfrbdVpVkTEc9CJwkCzq25PSM1isA88HTEJpF0DvB14Pci4md1t6cISe8HnouI79fdlpItAn4V+M8RcRnwf2nBr8ZZ3fU6YANwPnC2pA/V2yorwygCc2unbktaTCcofy0ivlF3e0pwJXCtpEN0Sk7vkvQX9TapFFPAVESc/o3mLjqBuul+E3gyIp6PiJPAN4DfqLlNZTki6TyA7M/nam7PSI0iMLdy6na2nN9XgUcj4gt1t6cMEXFLRKyNiPV0/p7ui4jGZ2AR8VPgaUmnVyu7ms4srKZ7CrhC0rLs3+PVtKBTM9M9vfkG4Fs1tmXkiqwul8uQU7eb4Ergw8CPJT2cbfu3EbG7xjbZ3D4OfC1LDg4CH625PYVFxIOS7gIeojNK6AckNI05L0l3AFcBqyVNAbcBnwPulPQxOj+Afru+Fo6ep2SbmSXGM//MzBLjwGxmlhgHZjOzxDgwm5klxoHZzCwxDsxmZolxYDYzS8z/B6CiTyxm+kk5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here we give you the asymptotic solution to the problem\n",
    "# we are using everything we have done before\n",
    "\n",
    "### Asymptotic behavior\n",
    "\n",
    "Qsource = np.zeros((n_y, n_x))\n",
    "A, b = build_2D_matrix(bc_dict, the_prob, D_matrix, Qsource)\n",
    "v = np.linalg.solve(A, b)\n",
    "n = n_x * n_y\n",
    "# array v contains the solution\n",
    "# we convert it in a matrix:\n",
    "\n",
    "c = vec2mat(v, n_y, n_x)\n",
    "\n",
    "# and we plot the matrix\n",
    "plt.contourf(x, y, c, 20, cmap=cm)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repeat the steady state calculation using a sparse matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x1164bccf8>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWYAAAD8CAYAAABErA6HAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHS1JREFUeJzt3X+s3fV93/Hnyz+pIYl/BQuundoEQ0PIFlqP0KJFLJTEWSLI1nQyVSISRbMqQZKm7SqYNpgcRcmqKD+kudushLZbUzxEosSKrBDEj3VdG2qHkASbGC7Gg2uIARuTdjBfX+e9P87X5Phw7z3fc77f7/l+vt/zekhHPuf783Owed/3fX9+fBURmJlZOhbU3QAzMzuTA7OZWWIcmM3MEuPAbGaWGAdmM7PEODCbmSWmb2CWdLuk5yQ90rVtpaR7JD2e/bmi2maamY2PPBnznwGbe7bdDNwbERuBe7PPZmZWAuWZYCJpPfDtiLg0+3wAuCoinpV0HvBARFxcZUPNzMbFoiHPWxMRzwJkwfncuQ6UtBXYCrAY/doqLR7ylnDWguFL4osXD3/uoqULhzpv4dLBv+sw5yxYumSg47V46cD30KIB7jHIsUBosL+bmZmZgY4/depUpdcf9HiA6enpgc85efLkwOcMe68yzh22vaedOHHihYh447DnX7BgWbxCvr/7n8b03RHRWxmozbCBObeI2AHsADhvwdL46OK1Q1/ronMG+x++28S5y4Y+d+WFK4c6b8WFg/+besObJwY+Z9n69QMdv/i8wY4HWLAm/zlatW6ga5/QYH+vx48fr/T4o0ePDnT8kSNHBjoe4Jlnnhn4nKmpqYHPATh8+PBQ5xW5Jwz3Hbs99thj/6fI+a9wirzx5rPTB1cXuVfZhk0jj2QlDLI/nyuvSWZm423YwLwLuCF7fwPwrXKaY2ZmeYbL3QH8LXCxpClJHwM+B1wj6XHgmuyzmZmVoG+NOSKun2PX1SW3xczM8Mw/M7PkODCbmSXGgdnMLDEOzGZmiXFgNjNLjAOzmVliKp+SbYN56YnDQ03LHsTJZw8NNS3brEnOWrAg/zIOx6pty6CcMZuZJWZsAvPh516uuwlmZrmMTWA2M2sKB2Yzs8Q4MOdwbHK4noEXJ58vuSVmNg4cmM3MEuPA3AIvHzpUdxPOEEefrrsJZqWTtFnSAUmTkl7zAGpJvyzpXkk/kvSApLVd+26Q9Hj2uqH33F4OzGZmfUhaCGwH3gtcAlwv6ZKewz4P/LeI+EfANuCz2bkrgduAdwCXA7dJWjHf/RyYzUas6LPwrBaXA5MRcTAipoGdwHU9x1wC3Ju9v79r/3uAeyLiWES8CNwDzPvgVwdmMzNYLWlv12trz/4JoLtGN5Vt6/ZD4Ley9/8CeJ2kVTnPPYOnZJvZrIo8ITsFixcvYOLcZfkOPsYLEbFpniM0y7bo+fyHwH+S9BHgr4DDwEzOc8/gwGzWAMMGycOHD5fckrE1Bazr+rwWOKMmFRHPAP8SQNI5wG9FxEuSpoCres59YL6buZRh1uXo0aN1N8HStAfYKGmDpCXAFmBX9wGSVks6HVNvAW7P3t8NvFvSiqzT793Ztjk1KjA/9g/TdTdhJF56wlmOWUoiYga4iU5AfRS4MyL2Sdom6drssKuAA5IeA9YAn8nOPQZ8mk5w3wNsy7bNyaWMnI5NHmPlhSsHPu/FyedZceEbK2iRpeDIkSN1N8FGJCJ2A7t7tt3a9f4u4K45zr2dX2TQfTUqYy7KK8yZWROMVWA2M2sCB2ZrtePHj9fdBLOBOTCbjdAws/6aPp7YBufA3BKpLWRk9fMY5ubyqAwza6VFSxfmH0n1k2rbMihnzAMY5YL5VY9lPvnsoUqvb83m8km9HJjNrHReQa+YsQvMHstsTeLMdTyNXWA2K8ugs/5GmUW646/ZHJgH5AezmlnVHJjN7Awun9SvcYF5XFaYA68yN+4cIMdXocAs6VOS9kl6RNIdks4qq2E2OE8yKaYtazG7vtx8QwdmSRPAJ4BNEXEpsJDO4tHJq2tkRpPrzD8/cqjuJpiNjaKljEXAL0laBCyj51ErbTVsB6CNL4/rtUEMPSU7Ig5L+jzwFPAK8N2I+G7vcdnTZrcCvN4zwM1ycX25uIVLFzf2IRVFShkrgOuADcD5wNmSPtR7XETsiIhNEbFpmRrX1zinUQ2bcwdgmlJ9cknR+nIZPxD820FxRSLlbwJPRsTzEXES+AbwG+U0q3qeAWhmqSoSmJ8CrpC0TJKAq+k8pND6SKUT0AsZpWnYrNWjMdpj6MAcEQ/SefDgQ8CPs2vtKKldjeBOwNFr4hNJ/Ku9DapQ0TcibouIX4mISyPiwxFxoqyGzacNk0yqypo9lrnZ6uz0c4djOtrTGzeEMurMo8ia3QFo/biM0S5jHZjNUuKM1U5zYK5RKp2AVp1R1JfLyJb9Q6E/SZslHZA0KenmWfZ/UdLD2esxSce79p3q2rer373GfsbH4edeZuLcZYWucWzyWP5ni9nIVNVRmOoY5hS0taNT0kJgO3ANMAXskbQrIvafPiYiPtV1/MeBy7ou8UpEvD3v/ZwxN4TrzNWqewEjZ6zJuxyYjIiDETEN7KQzwW4u1wN3DHuzsc+YyzJs1vzi5POlTxt9+dAhlq1fX+o1LU3u9JvbwqWLecObJ/IevlrS3q7POyKie/jvBPB01+cp4B2zXUjSL9OZEX1f1+azsuvPAJ+LiG/O15jGZsxlDpkraxZg1dO0q8iaPcmkOnl/ra97QomzdQBeOL10RPbqnZOhWc6JOa61BbgrIk51bXtTRGwCfgf4kqQ3z9eYxgbmVDUxOOflpT/z15errrU6Ux65KWBd1+e1zL2a5hZ6yhgR8Uz250HgAc6sP7+GA3MFUpgR6IkmzTBMtlpmUHa2nNseYKOkDZKW0Am+rxldIeliYAXwt13bVkhamr1fDVwJ7O89t5sDc6bsRY2GCc5NyJrbqOyOvypLGClnym0dkQEQETPATcDddNYEujMi9knaJunarkOvB3ZGRHeZ4y3AXkk/BO6nU2OeNzC7869Cw3QIltkZmLcT8OSzh1h8Xv/jBhFHn0ar1vU/sEHKHCbXtqA8DiJiN7C7Z9utPZ//wyzn/Q3wtkHu5Yy5SxVLgVZV1nDWPL+6FjuqKmusIii7jJGuRgfmpixmNGhw9ozA9OTJlqsqYThTHj+NDsxVqGoB/SqCc56s2Z2A8xv1xJJUgrKz5bQ5MI9QXcE5j7zjmT1kbnZVlDCcKY8vB+ZZVPnYqRSG0tlgyur0GyRLrTIoV5Ett3lERh08KqMGg4zWyDNK46UnDs879bQtU7TL7tArq4yRJyilEpTHyYKlSxr7777xGXNVHYBVP6zVmXN1mvj4qdOqDsquLTdD4wPzOChjlEaeTkCvm/Fa/coYZWfLZuDAPK+UsuZ+wXmU45rdAZhfSiUM/4BoDgdmG0t56steED8fd/yVz4G5jzZlzaMuZ8TRp/sf1GD9ApKzZRtWKwJz1TMAqw7OVp4mdvw5KFuvVgTmpisza7ZyFO30czC0IhyYc2pK1lxGOSOPJncA1v18v27Olm02DsyJSClr9rC5YsYpGLrjrxoOzANoStZsxZQxdjkPZ8s2l9ZMyX7sH6a56JwldTejkGGftN2rjVO0m9ipZ/XS4qWlPwBiVJwxD8hZc7ryBO+q68upZKmptMOG48A8hBSCs0dnNFvTVo+z0XJgTkwqixuNawfgqOrLbeD/FtVxYB5SClmzWS9ny+3QqsDclGcAjsKoxjPbL+QNil5v2fppVWAetVSeD1iXJk8yaSNny+1RKDBLWi7pLkk/kfSopF8vq2HWnzsAB5PSjD9rHkmbJR2QNCnp5jmO+VeS9kvaJ+kvu7bfIOnx7HVDv3sVHcf8ZeA7EfFBSUuAZQWv1ziHn3uZiXPH7mu30ig6/try1Otx6/iTtBDYDlwDTAF7JO2KiP1dx2wEbgGujIgXJZ2bbV8J3AZsAgL4fnbui3Pdb+iMWdLrgXcCXwWIiOmIGMtZAFWUNMooZ4xy8fy5tH3pzxS4hDESlwOTEXEwIqaBncB1Pcf8a2D76YAbEc9l298D3BMRx7J99wCb57tZkVLGBcDzwJ9K+oGkr0g6u/cgSVsl7ZW09+X4eYHb5eMOwPz6dQCO65C5YTg4Nt7q03Eqe23t2T8BdGcZU9m2bhcBF0n635K+J2nzAOeeoUgpYxHwq8DHI+JBSV8Gbgb+ffdBEbED2AFw3oKlUeB+SXNJo15NmLJdRRnDPxDmpkVLWLBmfd7DX4iITfNdbpZtvfFsEbARuApYC/wvSZfmPPcMRTLmKWAqIh7MPt9FJ1CPLY9ttlFyUB6pKWBd1+e1QG+hfQr4VkScjIgngQN0AnWec88wdGCOiJ8CT0u6ONt0NbB/nlNsQHnqzB6ZYTYSe4CNkjZkAx22ALt6jvkm8M8AJK2mU9o4CNwNvFvSCkkrgHdn2+ZUdFTGx4GvZQ09CHy04PUazyWNdkptFEKd2XJq/y1GISJmJN1EJ6AuBG6PiH2StgF7I2IXvwjA+4FTwL+JiKMAkj5NJ7gDbIuIebOuQoE5Ih6mMwTEzKzVImI3sLtn261d7wP4/ezVe+7twO157+WZfxVwrTk9/SaX9BvDXFSZHX+uLbdfKwOzh8zZKDlQWtlaGZhTUFbWXHSiSQqTTKw8/iEwHhyYrRAvZGRWPgfmCjWh1uzlP5vD2fL4cGA2q1ib1l8ex6FydWjNU7LH2YuTz7PiwjfW3QyrkLPlISxaglat639cgpwxV6wJ5QwzS4sDs83LK8yZjZ4Ds5lZYhyYzcwS48Bs1kfdIxHc8Td+HJjNzBLT2sCc0noZRUdmlPH8PzNrjtYGZjOzpnJgNjNLjAOzmVliPCXbzFoptIATWlJ3M4bijNnMLDGtDcwXndPMn5SzWXnhyrqbYGYj1NrAnBI/NdvMBuHAbNbH+eefX3cTbMw4MJslbu3atXU3wUbMgdnMLDEOzGZmOUjaLOmApElJN89z3AclhaRN2ef1kl6R9HD2+i/97uVxzBVresff4vPW190Es9pJWghsB64BpoA9knZFxP6e414HfAJ4sOcST0TE2/PezxlzC/h5f2aVuxyYjIiDETEN7ASum+W4TwN/DPy/IjdzYDar2MTEROFrpNIB2KQRKjMzMxw/fjzXC1gtaW/Xa2vP5SaAp7s+T2XbXiXpMmBdRHx7luZskPQDSf9T0j/t13aXMirUhDLGsvXr626C5bR27Vovml+dFyJi0zz7Ncu2eHWntAD4IvCRWY57FnhTRByV9GvANyW9NSJ+NtfNnDFbIQvWrK+7CWajMAWs6/q8Fuh+tM3rgEuBByQdAq4AdknaFBEnIuIoQER8H3gCuGi+mzkwV6SsbLnodOw3vLn4r9GWjlRKGmNoD7BR0gZJS4AtwK7TOyPipYhYHRHrI2I98D3g2ojYK+mNWechki4ANgIH57tZKwNzm9bJsPQ5WLZfRMwANwF3A48Cd0bEPknbJF3b5/R3Aj+S9EPgLuB3I2LexxK5xlyBJtSWx82qVas4evRobfefmJjg8OHDpVzLteZ6RMRuYHfPtlvnOPaqrvdfB74+yL1amTGbDWrNmjV1N8HsVQ7MJXO23E6pDROrs3yS2n+LNiocmCUtzMbnzTZ2zwrwOsxm46mMjPmTdIrhY6+ObNmz/saXOx3bq1BglrQWeB/wlXKaY9ZeZcwAtPFQdFTGl4A/ojO4elbZ1MatAK8fwSCQuobKNbG23G/WX5MWMFq+fPnpqbVjxSM05nbq1KnG/psYOmOW9H7guWwmy5wiYkdEbIqITcvUzr7GKoJyGfXlFCaXaNW6/ge1QJ1lBZc02qdIpLwSuDabfrgTeJekvyilVQ3SxEzZZjeKIXNtKWd4ZEa1hg7MEXFLRKzNph9uAe6LiA+V1jLryx1/g1m1alWh81MORs6a26WdtYURqSpb9jA5G4aDc3uUEpgj4oGIeH8Z1yrCa2T8Qr/6chnLfXpluTPlDYxtKWdYdZwxD8m1ZUuRs+Z2cGBOTCpljCYNlStTvw7AlOvMo+b/FtVxYB5CCtmyO/6arcpyhrPm5nNgHlAKQdlmt3z58r7HFB2ZYTYKrQnMbej4K6uMMYqOv1HLE3RTkEq2mko7bDheKH8AzpbHw5o1azhy5Mic+88//3yeeeaZOffnVebi+bMZ9+naMzMztT4coYjWZMxNN0i2XHV9eVw7/soyTtmqOwCr4cCcU1Oy5VGVMZo8hjmlOnPVY5rH6YdEmzgwJyClbNk6ig6bc0C0IloRmKvu+GtKtmzN6STs5qzZerUiMFep6qBcZrZcRhmjzPpy25f8dNZsVXFgtrGUp848yidnNzlrdgdg+RyY59GmbLlMTe74G7WUsuaU2tJEkjZLOiBpUtLNs+z/XUk/lvSwpL+WdEnXvluy8w5Iek+/ezkwj4kmTipJwSjXzvCqc+mStBDYDrwXuAS4vjvwZv4yIt4WEW8H/hj4QnbuJXTWrH8rsBn4k+x6c2p8YK6q469J2XJZ2jJ+ObUOwEEy1SaXNFrucmAyIg5GxDSdpzZd131ARPys6+PZQGTvrwN2RsSJiHgSmMyuN6fGB+a2yxOU2zgFezZlB9yyxjPnyZpTCs42lAng6a7PU9m2M0i6UdITdDLmTwxybjdPyZ5FldlyKst6Wn79pmhXocrp2lVM1S5rmnqZZmZmBvl7Wy1pb9fnHRGxo+uzZjknXrMhYjuwXdLvAP8OuCHvud2cMY/QoEG5jGw5r7aUMepSdtYMzpxH7IWI2NT12tGzfwroHv+5FpjvJ9FO4ANDnuvA3CuV5/iVVVcuu4zhERnFpBKcXWse2B5go6QNkpbQ6czb1X2ApI1dH98HPJ693wVskbRU0gZgI/B3893MgXkEqipfjHKIXFvlrTPnGdOcd4RGKsHZ8ouIGeAm4G7gUeDOiNgnaZuka7PDbpK0T9LDwO/TKWMQEfuAO4H9wHeAGyPi1Hz3a3SNuewRGVVky8ME5VGWMNpq+fLlHD9+fOT3rarWWkXNedyXBR1UROwGdvdsu7Xr/SfnOfczwGfy3ssZc4WqCsp55S1juL6cT5kzAYcpJThzHh8OzJmys+Uqg3ITsuW2r5Mxn6pKGpB2cPbU7PI4MFcghSFxbRm7PAoprc+cR5nB2Z2AaXJgLtmwQbkJ2bJHZOQvZ1SZNUPambMV19jAXGbHX1lljKqDchVcX65OU4Kzs+b0NDYwp2YUQTlvtuwyxvhw5txOjR4uN06a0OHXZKtWrar1icoeula+6enp5KaJ5zX2GXMZZYwUOvvstapaZW6UC+g3jUdmlGPsA3Od/GDV9htFoCqjnOE6c1ocmM0S4eBop411YHYZw9rCnYDt0sjAXNVTS0apitEY4BEZTees2aChgTkVzpYtj6Z0iPmHQjqGDsyS1km6X9Kj2VJ3c66sZGdKpdOviZNLUnueX0pczmiPIhnzDPAHEfEW4ArgxlmeGpusqh+2ajYsZ642dGCOiGcj4qHs/d/TWTx6bH5kj2r6tSeWpCnVscxFs+Yyfig0pXSTslJqzJLWA5cBD86yb6ukvZL2vhw/L+N2Zq3nrHm8FZ6SLekc4OvA70XEz3r3Zw813AFw3oKl8z4Z1qytUnyKdNudPHmysdPcC2XMkhbTCcpfi4hvlNOk6hWtLzdxFbmivOSn2egUGZUh4KvAoxHxhfKaZMPyGOZimrZg/lw8OqP5imTMVwIfBt4l6eHs9c9Latec2jC5xCwP15nH19A15oj4a0AltsV6eESG1cFLkNbPM/8G5Nl+ZlY1B+YRaXLHn81u0LHMoxzf6zpzs41dYPaMP2sS15nTIWmzpAOSJiXdPMv+d0p6SNKMpA/27DvV1Re3q9+9/GgpM7M+JC0EtgPXAFPAHkm7ImJ/12FPAR8B/nCWS7wSEW/Pe7+xy5iLGGV9ueqOvyYuYGTN0cJp2ZcDkxFxMCKmgZ3Add0HRMShiPgRUHiKswPzCLi+bE0zhiWU1aeXjsheW3v2TwBPd32eYrC1gc7Krvs9SR/od7BLGS3hySXWa2JigsOHD9fdjNpMT08P8v1fiIhN8+yfbWjwIEtMvCkinpF0AXCfpB9HxBNzHeyM2WyEhvkVfwyz1xRNAeu6Pq8Fci9+EhHPZH8eBB6gs+jbnByYzcz62wNslLRB0hJgC9B3dAWApBWSlmbvV9OZNb1/vnMcmK3V/MQTK0NEzAA3AXfTWXv+zojYJ2mbpGsBJP0TSVPAbwP/VdK+7PS3AHsl/RC4H/hcz2iO1xirGrPHMJvZsCJiN7C7Z9utXe/30Clx9J73N8DbBrmXM+acxnGpT+sv1SeZWLM1KjCPy8pyXrzIbLw1KjCbVa0tazJbszkwm7WYFzNqJgdmswaoYyyzx0/Xx4HZzCwxYzVcziwFfmL2aExPTzf2SSzOmM3MEuPA3AKpLWCkVev6H2Rmc3JgNjNLjANzDp71Z2aj5MBsZpYYB2Yzs8SMTWD2ynJm1hRjE5jNzJrCgTkxo1hZzk/INkubZ/6ZWSudPHmysTMsnTGbmSXGgdnMLDEOzGZmiXFgNjNLjAOzmVliHJjNzBLjwGxmlphCgVnSZkkHJE1KurmsRpmZpaZfvJO0VNL/yPY/KGl9175bsu0HJL2n372GDsySFgLbgfcClwDXS7pk2OuZmaUqZ7z7GPBiRFwIfBH4j9m5lwBbgLcCm4E/ya43pyIZ8+XAZEQcjIhpYCdwXYHrmZmlKk+8uw748+z9XcDVkpRt3xkRJyLiSWAyu96cikzJngCe7vo8Bbyj9yBJW4Gt2ccTn50++MjQdzw29JnFzv3JvHtXAy8UuHqq/L2ao43fCeDiIiefOHHi7scee2x1zsPPkrS36/OOiNjR9TlPvHv1mIiYkfQSsCrb/r2ec+ddFKdIYNYs2+I1GzpfbgeApL0RsanAPZPTxu8E/l5N0sbvBJ3vVeT8iNhcVlvIF+/mOiZXrOxWpJQxBXQ/dXMt0MwVQ8zM5pcn3r16jKRFwBvo/K4+cKwsEpj3ABslbZC0hE5xe1eB65mZpSpPvNsF3JC9/yBwX0REtn1LNmpjA7AR+Lv5bjZ0KSOrodwE3A0sBG6PiH19TtvRZ38TtfE7gb9Xk7TxO0FC32uueCdpG7A3InYBXwX+u6RJOpnyluzcfZLuBPYDM8CNEXFqvvupE9DNzCwVnvlnZpYYB2Yzs8SMJDC3ceq2pHWS7pf0qKR9kj5Zd5vKImmhpB9I+nbdbSmLpOWS7pL0k+zv7NfrblMZJH0q+/f3iKQ7JJ1Vd5sGJel2Sc9JeqRr20pJ90h6PPtzRZ1tHLXKA3OLp27PAH8QEW8BrgBubMn3Avgk8GjdjSjZl4HvRMSvAP+YFnw/SRPAJ4BNEXEpnU6pLfW2aih/RmeqcrebgXsjYiNwb/Z5bIwiY27l1O2IeDYiHsre/z2d/9Grf8R1xSStBd4HfKXutpRF0uuBd9LpNScipiPieL2tKs0i4JeycbPLaOBcgoj4K147N7d7evOfAx8YaaNqNorAPNtUxsYHsG7ZKlKXAQ/W25JSfAn4I+DndTekRBcAzwN/mpVoviLp7LobVVREHAY+DzwFPAu8FBHfrbdVpVkTEc9CJwkCzq25PSM1isA88HTEJpF0DvB14Pci4md1t6cISe8HnouI79fdlpItAn4V+M8RcRnwf2nBr8ZZ3fU6YANwPnC2pA/V2yorwygCc2unbktaTCcofy0ivlF3e0pwJXCtpEN0Sk7vkvQX9TapFFPAVESc/o3mLjqBuul+E3gyIp6PiJPAN4DfqLlNZTki6TyA7M/nam7PSI0iMLdy6na2nN9XgUcj4gt1t6cMEXFLRKyNiPV0/p7ui4jGZ2AR8VPgaUmnVyu7ms4srKZ7CrhC0rLs3+PVtKBTM9M9vfkG4Fs1tmXkiqwul8uQU7eb4Ergw8CPJT2cbfu3EbG7xjbZ3D4OfC1LDg4CH625PYVFxIOS7gIeojNK6AckNI05L0l3AFcBqyVNAbcBnwPulPQxOj+Afru+Fo6ep2SbmSXGM//MzBLjwGxmlhgHZjOzxDgwm5klxoHZzCwxDsxmZolxYDYzS8z/B6CiTyxm+kk5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Here we give you the asymptotic solution to the problem\n",
    "# we are using everything we have done before\n",
    "\n",
    "### Asymptotic behavior\n",
    "\n",
    "poro = 0.4\n",
    "prob = Problem_Def(n_x, n_y, poro, width_x, width_y)\n",
    "Qsource = np.zeros((n_y, n_x))\n",
    "A, b = build_2D_matrix(bc_dict, prob, D_matrix, Qsource)\n",
    "A = csr_matrix(A, copy=True)\n",
    "\n",
    "v = spsolve(A, b)\n",
    "n = n_x * n_y\n",
    "# # array v contains the solution\n",
    "# # we convert it in a matrix:\n",
    "\n",
    "c = vec2mat(v, n_y, n_x)\n",
    "\n",
    "# # and we plot the matrix\n",
    "plt.contourf(x, y, c, 20, cmap=cm)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now do the 2d problem\n",
    "\n",
    "As before, we will save a series of snapshots for plotting at different timesteps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-e8a4fef275e77cec",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "dt = 1.0  #days\n",
    "Tf = 2000.0  # days\n",
    "number_of_fig = 9\n",
    "plot_times = np.linspace(0, Tf, number_of_fig, dtype=np.int, endpoint=True)\n",
    "plot_timesteps = (plot_times / dt).astype(np.int)\n",
    "nTstp = int(Tf / dt)  # number of timesteps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up the problem domain and build A and b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "width_x = 10  # dm\n",
    "width_y = 10  # dm\n",
    "n_x = 25\n",
    "n_y = n_x\n",
    "n = n_x * n_y\n",
    "poro = 0.4\n",
    "\n",
    "the_prob = Problem_Def(n_x, n_y, poro, width_x, width_y)\n",
    "x, y, D_matrix, c_init = make_D_matrix(the_prob, Diff, decreasing_factor)\n",
    "Qsource = np.zeros((n_y, n_x))\n",
    "A, b = build_2D_matrix(bc_dict, the_prob, D_matrix, Qsource)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare dense and sparse\n",
    "\n",
    "You should see about a factor of 6 speedup when you set \n",
    "`try_sparse = True`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "lines_to_next_cell": 2,
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-80a44bbf4fb48e87",
     "locked": false,
     "points": 10,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "with try_sparse=False, elapsed time is 18.92159104347229\n"
     ]
    }
   ],
   "source": [
    "try_sparse = False\n",
    "v = mat2vec(c_init, n_y, n_x)\n",
    "Adelta = np.zeros((n, n))\n",
    "for i in range(n):\n",
    "    Adelta[i, i] = poro / dt\n",
    "Aa = A + Adelta\n",
    "Bdelta = np.zeros(n)\n",
    "if try_sparse:\n",
    "    Aa = csr_matrix(Aa, copy=True)\n",
    "\n",
    "save_figs = list()\n",
    "\n",
    "clock_start = time.time()\n",
    "fig_timesteps = []\n",
    "fig_count = 0\n",
    "for t in range(nTstp):\n",
    "    capture_time = plot_timesteps[fig_count]\n",
    "    for i in range(n):\n",
    "        Bdelta[i] = v[i] * poro / dt\n",
    "    bb = b + Bdelta\n",
    "    if try_sparse:\n",
    "        v = spsolve(Aa, bb)\n",
    "    else:\n",
    "        v = np.linalg.solve(Aa, bb)\n",
    "\n",
    "    if t > capture_time:\n",
    "        save_figs.append(vec2mat(v, n_y, n_x))\n",
    "        fig_count += 1\n",
    "\n",
    "clock_stop = time.time()\n",
    "print(f\"with try_sparse={try_sparse}, elapsed time is {clock_stop - clock_start}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot your saved snapshots on a grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-be8acbb200b99a1f",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn0AAAJdCAYAAABH38H7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xe8HFX9//HXJ8lNMATSCJEkwJUmTUGJgAWNgg0LfL+igoBg+eJXv1ix4NcC9vK1f20/bCACioiCfC2oEAuoEGx0DHAxjUtIgxBIIZ/fH+dM7ty9O7uzu7P9/Xw89nHv7pyZOTt7ZuYzZ845Y+6OiIiIiPS2ce3OgIiIiIg0n4I+ERERkT6goE9ERESkDyjoExEREekDCvpERERE+oCCPhEREZE+oKAvg5ktNDM3s7PbtP7BuH43s8F25KFdzGyymX3EzG41s4dT2+HgAtexIFluUcsUaQUzOzeW3XPbnZdmMLOh+P1ObXdeOlWnlAH9Vt1nQrNXEIOms/Kmd3drXm46QyqQPNfdh9qYlU71A+DF8f+HgeH4/+Y8M2v7Sl5m9jZgGvATd/9bm/OyAFgADLn7ue3MSzOY2bHAwcDf3P0n7c6PZItB3CCw0N0XtjUzUqimB30lhqsn6QtJELwQGMpIsxm4PfV/XzCzfRkJ+I539x/UsZg821cE4G3A7oRy0tagjxDwnQX8Fji3QroVhGPDiuZnqVDHAqcA5wEK+jrbqcCz4v8LK6S7E3gEWNfk/EhBWhr0uftjW7m+bubuy4B9252PNnhC/LuqzoBPpKe5+3uB97Y7HyLufmS78yC1UZs+6TST49/1bc2FiIhIr3H3pr6AswEPq6ppvicn8wFPrJL2/Jju12WmbUe4hXMtsIZQFX0P8F3g4ArLXBiXeXaZaUm+FtQyP+GWjVd4DaXSDqY+H8xYx1Tgg8BfgAcI7d/+CXwN2KNC3rblH9gB+ChwW5x/FXAFcFiDv3tN2z1dTjJe5+ZYZy3bd0G6XAJ7Ad8GlgAbgaXAN4C5VdY5nnAr5JeE5gubgJXx/fGANbgd9wO+AtwCPEgIhm8Hvg+8DBhXYJkfitvkVGAi8C7g78BDhNs3VwEvyJHnw4DvAIvjvA/E/H8beF6F+Y4l3PZbHrfjGuB3wH8CA9X2M8CA/wD+HNf5IPBH4KQ6yptn7YvAnsA5wN2xrKTL1dT4u18A3AisTm3/C4HDy+RlsFI+kt+kTDnP3CcI5fuHwLKYx/uB3wCvAcZnzJNsk4Xx/ZHA/xHK8yPArYTbz9vVWIYX5Ph+C1Lph2igHKaXCewMfA64A9iQ/l1T6Z9E2Dfuid9zDWHfeRswKc+2qva9K6R5JvDT+Ps8TNi3PwZMid9/1HErqwwAxxH2hdXxe/4NeCtljg9VfqtknZVeg+V+qyq/w8z4O9wZv+c9wJeBWan0uxPOXXfH3+FfwGeBHarkeSrwPsJ+v4ZQ3pcAF1Fmf+v3V/NXUGfQF+e9Mc77PxXSbE84ETpwSsm0uallOOFEsjb1/lHgzRnLXUjxQd8XgXtT86+O75PX9am0g+V2stT0A2LBTtI8TDjRJe8fAV6WkbckzQmEIDGZ/6GSbfX8On/zmrc78M64Ddal0qS3zRdzrLeW7bsgle7ZhADB4zbcnJq2jIzAD5gN/CmV1ku+pwOXARPr3I7vidsh6zd2YFqBZX4opjk99b02pbaNA1uB12bMPz7+Bun8rS8pV2vLzDeFcOJLz7curit5fy0wvcJ+9hFCwOjx91tXsrwPZZS3R1PrS5eVezP2xVeltsdD8fsNlTvexdeDhH0xvf3eUpKXXeM6k+PYptK8AK9MpT+XCkEf4eSaXt8aYEvqs99Q5kSayvtCQqC1NTV/+re4iozAMSM/T4vf4WFGynHp93tageUwSfN6Ro4H2/adkrRvK/lua+O6kvd/B3aptK0qfO8FyXIypr+5zLo3xv9viXlzqgR9hOAp2bfXlJS/82o85rwybrNkG6wv81vtWua3OrXC7/BqRs5V61PfMfme04CnEC4ukn0xfQz+Q1Z5I1xgpo/5Wxh9jNwKvLee42+vvpq/gsaCvnczcuIte8UCnJQqTFNSn49n5ICxFjiRePIF9mD0SeaFZZa7kIKDvhrnH0ylGyyZtgNwV5y2FDg62T7AQYTaDSecbA6qsP7VwM2EoGccoabkKYRaP487dK1Xio1u91OTdTdQ5vJs3wUl2+EyYN84bSLwitTB47tl5p8IXBen3xB/g8lx2vaEA91wnP75Or7DG1P5u4xUDR0wA3guobZvxwK3/VBqeywFjiHWsAGPT5WrB4GpZeb/VGr53wL2SU3bOS7v+2Xm+3Gc55+EC5Ed4ufbAS8l1A448OMK+9nq+J1PAR4Tp80DLmfkhLh3he98as598cG4jeenpqe/538Sgq7DiAE5Yb96HPAFwkloC/CkMus5myqBREx3LhlBHyFQSvL6/4DHpsrk2xg5mZb7HZL1r4nb6+PATnHajsCHUssuG3DVm++Cy2H6t7oNeA4jx8f0b/XiVNqfAI9L7dsnM7L/X0NJ0JHnt6JC0EcIhJMLjiuTfBHa2R9HuOOymupB32pCEPV24rGAUKv2jdR3e04dv9VCMs5fefef1PrXAH8l3j0CBgi14cnF4P/G5fwGOCC175/OyMXK6zP2yyTI/SHh7uCE1PHmw4yU92Nr3Qa9+mr+CkZf+ZZeMZS+vlgy75zUjlH2thDhNpoD55d8/srUesfUWMWdKzlB3lhLoU8td0GF793o/IOpdIMl094TP98EHFhm3h0IVeQOXFFh/fcBO5eZ/oRUmqfX+Hs3ut1PpfVB31WUv0365jh9Q3IwSU37rzjtJjJuPwCHEE7yG8tt5wp5m87ICecict4iLmDbDzFysbBvmemzGKmtObFk2j6M7KufquG7vijOs4LsGtV5jNSCHVwyLdnPHHh2mXknES4aHXhfhe98aoU8pvfFIVIXl3WUzaRW5ptlpp1NA0Ef8BhCsODAhRnzvjn1XeZnrD/zZA/8KE7/VR3fvWy+iyyHcXryHdYB8yqs5+aY7veUqUkCXpJa1nG1/lZUDvp+HafdTJlbyIQL8W1lrsK2zCy7wKI4/Rt1/FbJflW2HJT5rcbkIZW/e4GZZaZ/OJXmpozt8N04vVzTrR/GaWMuylNp3h7T/K3WbdCrr1Z35Jhd5TU1ndjdlxNOyBCuvEYxs10I7U4gtOtLe2X8+0d3/2XpvO6+hXDlCnCgmT2hNE0HS77bJe5+U+lEd38Q+HR8+0Izm1qaJjrH3e8rM/+NhKAR4Il15q2btvvH3X1rmc8vi38fA+xdMu318e9X4/Yew91vIBzUJxIO4nkdRwjcNwPv8Hj0yqGobX+Ju99WZv6VhFoWGFsuTiHUFq+ihnE5GdmO53vosT6Guy8Fro5vn5+xnGvc/erSD919I+HCsFye6/Fld2+kk9H/xb/PKCAvpZ5LqAWGEJSU81VGhno5ISPNRuAzGdOSfaKIbVlNPeUw7fxYdsYwsycC+8e3H3H3R8us56eE2nzI3lY1M7MZhNpHCE2XNpZZ99WEYLSaJYTAqJzL499W/FaVfMPdV5X5PH2M+ly57UDGvhu34b/Ht5+ssO5k2xxkZrPzZLbXtXrIlnoGXv4ucBTwb2a2vbs/lJr2KsItreWEK6e0+fFv6edpVxNqJ8bH9DfWkb+WMrOJjOwAlb7br+LfcYRq7zEnRELD1yzLCbekZlRIU043bves7bA89f+27WBmOzDyG3zEzD5YYdnJfLvXkJ+nxb83uHstY7EVte2rlQsYWy6SPP/K3R+pks+0JPg5zcxeXSFdcuGStR3ryXM9rqmWwMz2AN5ECPT3JATwpRfY8wrIS6nk91/i7neUS+Duj5rZVYRb//PLpQFurhDYFrktq2n0N630WyXffQthXMQsvwIOJXtb1eNJhFv+VFn3QuCIKsu6PuOCFVr7W1VyXcbn6XF7r6+SZnrJ509lZJ+6yixXaLE7Giu45YMz1+NSwtXpFEJkn67RS2r/LihT8HeOf8vWHgC4+yNmdj+hlnHnrHQdZgbhhA0VvhuhLUwi67uVraGKtsS/AznzVbqurtnuFWrqtqQOJunt8FhGDjh5D6iTqycZtXwIPdxqUdS2r6dc1JxnMxsAdopvp1JS058hazs2oyyXM6ZmPM3M/o1wS35S6uMHGOnMMZFwAtu+gLyUqvr7R8mxoZHfvxXnjkZ/00q/VfLd78+oYUpU21b1mJX6f3lmquq/I7Su3DciK49bakhTWt7mpP7PW4NXyzG4Z3X8OH2xZu/H8e22moB4a+qg+DarehvCgTbXqmrPXdtVyrNn/N8qvbzdx6f+P9zdLcfr7DrWU++2aee2r2WZ6e14fM7teGrB+a3VmNuACTObSWhrNYnQLGUBoXPPVHef7WFw+pe3II+9vO/VIvO3SmnHtkpXS1Vabs8/krQBybHj4ZzHDXM9Tg7ogqAvSoK655jZ3Ph/Usv3t3Lt2hi5yts1a6Fmth2hpxOE7uJ5JQeT7SqkyVNrUY/VqfVnfreSabV8t0Y1c7t3ivQtgma0SUxu6Q7WOF87t33NeY63gZPHN3VC285GHU3o5boGeIm7/9bdHy5J08ynElX9/aPk1nI37ntFSbbVLDObVCFd1rZKaqDqOQekayDnZKSpNq3f3Rv/PsbM9mprTrpMtwR9VxGq2ccBrzKzcYT2fJBdy7co/q30mJgFjFQbZ7UpKGdN/Fv24BrbfO1XYf7k6q7mKzl33wT8I76t9N2Oin+3EgZvbpVmbve86t6+uRbuvoYwvhSEoQeKdm38Oz92Vsqrnds+yfNzY2CZV9Lu6uVxv261pFlIEWUlOR7c7u4bMtIclfF5EXlJfv95ZrZPuQRmNp6RTkXN2PcqKXJbNyrZVhMYecZsOcnvVbqtKp4DosMyPv8rI8eoBRXmrzSt2TrptyrnWka2YTOOwT2rK4K+2F7vgvj2ZELPp7mEGq8LM2b7fvz7VDN7XulEM5tAeJoFwE0ZtYVZ/h7/vixj+jsZ3aan1APx77Qa1pmWfLfjzOzA0olmNoUwxiHAz9y9lQ/DbuZ2z6vR7ZvHOfHvkWZW8aATe5rV4oeE7zAB+LzlbKVMe7f9uYT9cSYjPYTzSLbjPoQBgTOZ2faxI1ORiiwryX62T7nA18wOZuRitRl5+RWh9zRk9959AyM1SBfVuZ56tWK/zMXd/8HIhdv7YzA8ipkdzUjgVrqtknPAHDM7vMy8OxOeDlNu3asZ6Vh3RrkybWbPpHonjmbqmN+qnDjqRNKT/F1ZFzmJOo7BPasrgr4oqdF7AvCJ+P+V7p7VG+dHjPT+utjMXhUbjmNmj4vTnxqnv7vM/JUkB4Dnm9mHzGzHuNydzOzjwPsJA8VmSU62J5pZPY1Lk0fVDAA/N7MXJrUksa3jLwk9bzfFvLRSM7d7Xo1u3zy+zsj3PN/MPmpm2676zWyymS0wsy8TBhfOLQbpybZ5JfDjGDAky55uZi8ys8uSshe1bdu7+2Lgf5Jlm9k3zWzbMDdmNsvMXmlmPy6Z7zJG2ux+0sy+lj6Am9lEMzvMzD5F6CRSdMefpKwcZ2alPQRrdSWhhmQGcEHSFCV+h1fE6ZUa3id5OcDMnlYhXVnxVvLZ8e0JZvb1ZJiKWB7fTBggGuAHcUihVkq+3xFmtm+L113Oe+LfI4BL4j6CmQ2Y2YmMHOevJQzenHYtI52WzjWz+RaMM7MFhJ63lc6vZxFqqg4ELk/2FTObYGb/TthX11SYv9mS3+roVJOqTnMG4SJnR+APZvba9PBk8Xz872Z2Ka2/wOlczR4IkNoGZx71OJ4yy7ohtSwnNP6utO65hMKbpN/I6MfUPErJI5FS8y4kY3BKQiPSq1LL2Upoa5c8tuidVeY/KTXvJsKt6yHgD6k0g6k0g2WWcWCcL0nzMKMfPfUIJQOKpuZN0iyosO0y85/jN29ku59KxoCkNaw/z/ZdkKSpsqzMbUXoefqbkjK5jrGPrdpc5/d4L6Mfw7aBfI9hq3fbD1FhsNeY5tyY5tyM/eLLJfl7kOqPYZtMOCin51vPSPvV9OdzS+atWk6pMJAu4dmnyW+1hdCbcghqew52Ku0nS/KbfqzXXYSavrLljlCze1tq3tVJXkjty5V+gzi99DFsqxn9WKurqPIYtgrfb0FW/nOU5+mE9mxJPlamvt/hqXSNlsOqx7dU2rczel9Nnt2avP8HMCdj3ucz+pFtDzEyaPQdhNuOmduKkcespded9PK+MTX9tlq+fyrNqdR5LCWMS5p8l+SRmMlvNS+VLvO3qvY7kO/58hXLG2H4m7vLlPf0I/ucOgYT79VXpw3OPJswpEGWdPu9Bxip3i3Lw2Cv84F3EJ5E8DDhBLOEMPTLIe7+pVq/hIeBPF9EuFq7jZEd/0rgue6eNbBpMv/3CLep/0A4ke9CGEMo99hdHm7NHUA4UP+NcMKaRKhV+jrhcTaX1PK9itKs7V7D+hvevjnXcz+hzc8xwCWE7zeJMJjzMuDnhEcJDda5/E8Qeqh/A1gcPzbCQ9kvIgxh9EDJPG3b9u7+qLufThh77wLCA9MHCPvHzYRHs41pEuHuG9z9BEJbs/MJwdE4wjBN9xGClHcTHqOWZxiLWvL8O8K+/GtCwD6bUFZ2r3N5ZxJGGbiOsO0HCL/dxwknqMwhOjwMnn0k8E3CyXT7VF6m1JCHdxCawPyI0OloCuEkeDXwWsIxqlKNY1N4aAv7TEIzhGWEjg7J96ulHWiRefo8YX/5HmEfmUz43f5E2IcO9fCQgHLz/pJQS3gFIWAbH5fxScLTeO4tN19q/i8Qgpqfxfm3I/zuHwUOZ6Q9XaW7Rk3h7v8k7I+XE4LzmYz8Vh0z1Ju7/5UwyPbphH34fkbGxfwnofnX8YwM5Nz3LEbLIiIi0iHM7AJCzfC33f117c6P9IZuatMnIiLS82K71qR26hftzIv0FgV9IiIiLWZmHzaz081st1RHvO3N7JWEW/HbEZoPlXYiEambbu+KiIi0mJn9hNAeGEJHmwcJQ6QklTHLgBd4c4a1kj7VMQ0yRURE+sjnCR17nkbobDaDEPjdQegc8mUPY/qJFEY1fSIiIiJ9QG36RERERPqAgj4RERGRPqCgT0RERKQPKOgTERER6QMK+kRERET6QF8FfXEQzPVmNr7deekUZuZmtle78yHlqcwWy8yGzOyodudDylN5L5aO71Kqp4O+0gO8u//L3ae4+6MtzMMCM1vaxOWbmX3KzFbF16fNzKrPKZ2oT8rsu8zsJjN70MzuNrN3lUwfMrOH48l/vZldWTL97WZ2r5mtM7Nvm9mkZuVVmqtPyvvZZrY5VZ7Xm9keqekHm9kNZrYh/j04NU3HdylUTwd9feI04FjgIOCJwIuBN7Q1RyKVGfBqYDrwAuB0Mzu+JM1L4sl/irs/b9uMZs8HzgSOBAaBPYAPtSTXIvX7Qao8T3H3uwDMbCJwGfA9wv5wHnBZ/Bx0fJeiuXtPvoDzga3Aw8B64N2Ek4QDE2KahcBHgWtjmp8CM4ELgAeA64HB1DL3BX4FrAZuB16RmnY0cAthRPVlwDuB7eP6t8blrwfmEILtM4E7gVXAxcCMuJwkj6cRRmtfAZxR4XteC5yWev864E8V0r8rLnM58Nq4rr3itBcBf43ffQlwdmq+/wPeXLKsfxAOSEYYXf4+YF38/MB2l4Fue/VLmS3zvb8E/G/q/RBwVEbaC4GPp94fCdxbYdknA/fEPL8vvWzgUOCPwNqY5y8DE+O0rwCfLVnWT4G3xf/fE7fZg3G7Htnu8tNtr34p78DZwPcypj0v5sVSn/2L8Pg10PFdr4Jfbc9AU79cyckj44CyGNgTmBoPCHcARxEeUfdd4Dsx7fZxR3lNnPZk4H7ggDh9BXBE/H868OT4/wJgaUm+3gb8CZgHTAL+H3BRSR4viut8ArCS7JPgOuCw1Pv5wIMZaV8ADAMHxmVfWHJQWBDXN45wVTkMHBunvQL4c2pZBxEOhhOB5wM3EJ4bacB+wC7t/v278dUPZbZkuUY4Ef1nyTYYjsu4EjgoNe3vwCtT73eK655ZZtn7E07iz4x5/hywhZGg7xDg8LhtBoFbGQnqDiWcOMel1rMBmA08Pm7XOanvv2e7y043vvqhvBOCvnWEQPRm4I2paW8Hfl6S/gpiEImO73oV/NLt3XDAuNPd1wE/B+5091+7+xbgh8CTYroXA0Pu/h133+LufwF+BBwXp28G9jezHd19TZye5Q3A+9x9qbtvJBwUjjOz9LOQP+TuD7n7jcB3gBMyljWFcGBIrAOmZLT7eEX8vje5+0Nxvdu4+0J3v9Hdt7r7PwgHtWfFyZcBe5vZ3vH9yYRbFpvid9+BcJVt7n6ru6+o8P2lMd1eZtPOJpyEvpP67ETCiXV34Grgl2Y2LU4rV94hlL9SxwFXuPvvYp4/QKjRAcDdb3D3P8VtM0Q4sT8rTrsuLvvImPx4YKG7DwOPEgKB/c1swN2H3P3OHN9V6tPt5f1iQqA0C/gP4INmlqQtLc/E9ztkTNfxXRqioC9c7SQeLvN+Svx/d+AwM1ubvAgnp8fG6S8j3D64x8x+a2ZPrbDO3YEfp5ZzK+FEMjuVZknq/3sItxzKWQ/smHq/I7De3cs9VHlOmeVuY2aHmdnVZrbSzNYB/0mo4SAe+C4GTjKzcYQD3Plx2lWEW2NfAYbN7BwzS+dJitXtZRYAMzud0LbvRbF8AeDu17j7w+6+wd0/Qbj9ekScXK68Q7hlV2pUeY8nwlWp9e9jZlfETiEPAB8nlvfoPOCk+P9JjJT3xYSaoLOB+8zs+2ZW8btKQ7q6vLv7Le6+3N0fdfdrgS8yEoiWlmfi+wczpuv4Lg3p9aCv3I5RryXAb919Wuo1xd3fCODu17v7McDOwE8IO1BWHpYALyxZ1nbuviyVZtfU/7sRbjWVczOhKj5xUPysnBVllpt2IXA5sKu7TwW+TqjOT5xHOIgeCWxw9z8mE9z9S+5+CHAAsA+hbYnUrh/KLGb2WmKHDHev1nPSGSmH5cr7sLuvGjNXSXk3s8mE9mCJrwG3AXu7+47AfzO6vH8POMbMDiLU1PxkW4bcL3T3ZxCCAwc+VeU7SHl9Ud5LlJbnJ5bU3D2RkWO4ju9SqF4P+oYJvfuKcAWwj5mdbGYD8fUUM9vPzCaa2YlmNtXdNxMayiZDDgwDM81sampZXwc+Zma7A5jZLDM7pmR9HzCzyWZ2AKGNyg8y8vVd4B1mNjfWNpwBnJuR9mLgVDPbP54AzyqZvgOw2t0fMbNDgVelJ8aDwFbgs8SrwJj/p8SryAHgIeCR1PeX2vR8mTWzEwm1as/12IsxNW03M3t6zN92FoZz2Qm4Jib5LvC6WIanA+8nu7xfArzYzJ4Re0N+mNHHvB3i915vZvsCb0zPHIPR6wll/Ufu/nDM4+PN7DkWhop5hFDbpPJen34o78eY2XQLDgXeQridCqHN4qPAW8xsUqz9Brgq/tXxXYrlHdCwsFkv4BhCT6i1hJ5ag4xtJPz6VPqPAuem3h8FLE69fzyhl9NKwm2iq4CDCY1dfwGsYaRH2TNS8307pl/LSM+wdxB6lz1I6CH28Zg2yWPSM+xe4N0VvqMBnyY0El4d/7cK6c+MyyzXu+s4wi2BBwkH0C9T0uuMcJJ1YI/UZ0cSenStJzScvgCY0u7fvxtffVJm7ya0E1qfen09TjsglqXkVuxvgPkl87+DcKJ+gNCWalKFdZ0St2e53rvPJNT0rQd+TwgK/1Ay/0nxuz079dkTgevidlgd95U57S473fjqk/J+UVz2+lje3lIy/UmEjhIPA38BnpSapuO7XoW+LP6o0iHMbJBwUhzw0FC5o5jZqwlDCDyj3XmRztDpZbYRZvZMwm3eQXffWi299L4eL+86vve4Xr+9KwWKtwzeBJzT7ryINFu8nfVW4JsK+KTX6fjeH6oGfRYec3Sfmd2U+myGmf3KzP4Z/05vbjal3Sw8CWEl4bbahW3OjkhTmdl+hFt9uwBfaHN2RJpKx/f6lYuRSqabmX3JzBab2T/M7MmtzuOo/FS7vRtvb6wHvuvuB8bPPk1oEPpJMzsTmO7u72l6bkVEREQ6RLkYqWT60cCbCcMFHQZ80d0Pa20uR1St6XP33xEakKYdQ+jeTfx7bMH5EhEREeloGTFS2jGEgNDd/U/ANDPbpTW5G2tC9SRlzfY4Ire7rzCznQvMk4iIiEguZnZIs5bt7jc0uIi5jB40e2n8rC1PNak36MvNzE4jdG1nADtkr6lT2G7WNMY9Zgc22wQeeugh1qxZA8CmTZu2zbd58+ZmZ016zMaNG+9391nNXEdpeZ5pA9umbTdupOJ8YCD8P2HS+G2fjZ80UPb/cZMmjix/YNLI/xNGPg8LG3nvNrqSfsuW0Z0IH3300YrTS98n0vtgotK+WC59PWnK6fdjQCvKM4wu09tvv/0h++67b9l0y/5y45jP0mUeRsp9Il3+E+myX+59en/YlsfUfgFl9g0YtX8kSvcTKF/2S/eXSmkrfQ6Vy3ueMl3r/lLv/lWqE/a3Aw8cc3d0mxtuuCFzf5iELZrSpHDGzNYThu5JnOPutXR2KffIvLYNm1LvVho2s11iLd8uwH1ZCePGOQdgl3GT/PwjDmPw+Qcz/dhTuW/CTK699lquu+46li1bxtKlIwPzL1+ed3BzkeCOO+64p3qqxpSW59cMzNs2bZ8pIyeduTtPBmDGXjO2fTZ9r5Hj1dQ95277f/Lg4Lb/B3YZ+X/c7JH/AWzmyGD7G230CW7t2rUV369aNfqBFcPDw5RTbr9L75elli1bljktz/yV9PsxoBXlGUaX6fnz5/uiRYvKpvvvSXuO+Sxd5mGk3CfS5T+R3g9g9L4Ao/eHRHq/gLH7BozePxKl+wmM3TeyPoOx+00ia/+ByuU2z76QZ5+qdZl5dML+llX2AMwsc3+YbgOkj8VF+sSmu2539/kNLGIpo5+UMo/8T28pXL1DtlxOGPSU+PeyCmlzmzevOT+aiIiISBtcDrw69uI9HFiXNI9rh6o1fWZ2EbAA2MnMlhIe7fJJ4GIzex1hNPWXNzOTItL5bDaKAAAgAElEQVS95s2bV1hthIhIJ8mIkQYA3P3rwM8IPXcXAxsIj+xrm6pBn7ufkDHpyILzIiItMnv27Iq3qESkPebOnVvzLV5pnwoxUjLdgf9qUXaqavsTOWbPnt3uLIiIiIj0vLYHfSIiIiLSfAr6RDrU1uGhdmdBRFpMt3almRT0iYiIiPQBBX0iIiIifUBBn4iIiEgfUNAnIiIi0gc6NuibM2dOu7MgIiIi0jM6NugTEZHutnrx6nZnQURSFPSJiIj0seXLl7c7C9IiCvpERERE+oCCPhERkQKoxkw6nYI+Ealo3rx5Dc2/dOnSgnIi/WDD0FDVNHpajUh9FPSJiIiI9AEFfSIiIiJ9QEGfSBea5JvanQUREekyCvpEOsjmFUPtzoKIdAm1l5VaKegTEZGu56uWtDsLIh1PQZ+IiEif0jAz/UVBn4iIiEgfUNAn0kSd8uzRmTNntjsL0geW3beh3VkQkQoU9ImISNusu3NZ1TT1dnDK28t92rRpdS2/aMuWVd8WIo2Y0MjMZvZ24PWAAzcCr3H3R4rImEi/2DA0xOTBwYaWMW3aNNauXVtMhkR6mPaV3rPduHHsM2VicxbeGTdrClN3TZ+ZzQXeAsx39wOB8cDxRWVMREREmkedOPpPo7d3JwCPMbMJwGRAJUgkhzy3tBrRKberRCQoekw9jdEn9ag76HP3ZcBngH8BK4B17n5laTozO83MFpnZog2+tf6cinSARsrzmsUrm5Qrkfqly/TKld1TRrcOD7U7CyJdp5Hbu9OBY4DHAXOA7c3spNJ07n6Ou8939/mTTf1GpLupPLeGbju1TrpMz5o1q93ZEZEmauSsdRRwt7uvdPfNwKXA04rJloiU0hMHpF9sGBqqaz7tI/npwqo/NRL0/Qs43Mwmm5kBRwK3FpMtEWm22bNntzsL0ofUzEGkfRpp0/dn4BLgL4ThWsYB5xSULxERkYblHatPpB80NE6fu58FnFVQXkRERMravGKIgV0G250NhoeHm7LcWgZmbrTnrm7t9i+1RBcRERHpAwr6RPpEnufvzpkzpwU5ESmGhm0RqY2CPpEuVdpWqXRA5k4YoFkDyEqrtasHb7fcMu2WfEpzKOgTkY6ik1L/afYTakq1+oJIFz/SKRT0iXQw3b6SflXvWH3l9FIP3kYCSF1QiYI+ERHpCptXDI35TBdGIvkp6BPpI3k6c4gUafXi1S1fZ7c9maOW4VrqpVo+AQV9Ij2tWtulPE/lyLqd1IoTlXS2O9b3zm1TkX6goE+ki1SrweiEHrsirVbuFm+31fY1k2r5JNHQEzlEpL0m+SY22sSKaaZNm8batWsLXW+1Wr56GpvrxCR51Ptkjjz7StqqVavGfJb1NI6ssptnP8hbY15vBw7tV5Kmmj6RDlOusXqRKrXrKz1B1HOiqXWe5cuX68TUZ9YsXpkrXS09eJtd21d0wLds2TIFfNJyCvpE+kCRt30rnajqCfhEoPVj9RWhkYAvLwV8UiTd3hXpMc24nZvIOlkp2JNOsHV4iHGzB0d95quWYDN3BUbf4q20n5Te2i1Xy9fsgE/j8UkzqKZPpE/UWttXetIpIuDTrVwpQrObQFSjgE+6lWr6RPrQzJkzyzZUz1JUwCfSbOVq+2qRp5avnEr7goI96RQK+kS6XK29EouS9wSlk5F0ulqaQ5Qrzwr4pFvo9q5ID2q040b6BFTupKWAT1olqwdvpVu8pT156+3FW1rLV0vAl7d37tKlSxvqrKF9TGqhmj6RDld6uyrdML1WzezkkaYTkbRbo7d5S9Ua8FWjmr3iDAyMY+7Ok5uz8NY/RbCpVNMnIjWpdrLSCUk62STPfnRcLe1cFfBJN1LQJ9IBahmENq8ixuar9fm6OiFJI2odq6+ZvXjTt3bzlus8T6rRrVxpJ93eFZHcKp2wdEISgGX3bWjerbYaFX2LN1FuPyhy0PKE9ikpWkM1fWY2zcwuMbPbzOxWM3tqURkTkfwq3bJqBZ2cpJmaURNer1YEfKrVk2Zp9PbuF4FfuPu+wEHArY1nSUS6iU5OUo+8z9+tphm3eGu5tduMgE+kWeoO+sxsR+CZwLcA3H2Tuze/W6CI5FbkM3cbaXwu/W314t7oAlnLPlDPowkV8EmzNVLTtwewEviOmf3VzL5pZtuXJjKz08xskZkt2uBbG1idSPupPI+mk1T3S5fplSuLqX3rRPWO1VeJnkUt3aaRoG8C8GTga+7+JOAh4MzSRO5+jrvPd/f5k02dhaW7qTxLr0mX6VmzZrU7O31FtXvSao2ctZYCS939z/H9JYQgUEREpGWaOXRLIm/tnZ5WI52s7qDP3e8FlpjZ4+NHRwK3FJIrEWmZ2bNn1zWfTlrSSu3owVupjOvxhNKNGr0/9WbgAjP7B3Aw8PHGsyQiRSqyM4dIs9U6QHO3UcAn7dRQ0Ofuf4ttQZ7o7se6+5qiMibS6yqd3KrdrirXKL3WsfpmzpxZU3oREeluaokuIiJdrxXt+irJc2tXtXzSbgr6RPqMbveKNKbWZ1KLdAoFfSIi0rO2Dg/VlD79NI4iqZZPOoGCPhGpSk/jkH5WrfxXm66ATzqFgj6RLlBrbYWIiEgpBX0iIiIifUBBn4iIiEgfUNAnIiIi0gcU9ImIiIj0AQV9IiLSFdrx/N1GqeeudBIFfSIi0hUmDw62OwsiXU1Bn4iISJPMmTOn3VkQ2UZBn4iIiEgfmNDuDIiIiIjUa8Kk8czYa0ZzFn5bcxbbLqrpExEREekDCvpEusC42YPtzoKIiHQ5BX0iUtW8efPanQWRtqlW/qtNV2cO6RQK+kREpGfVWks+e/bs5mREpAMo6BMREanB3Llza55HtX3SCRT0ifSZtWvXtjsLIoUb2GWwretXEwjpBgr6RDpQPSewjTaxpvSrVq2qeR0iUj/V9km7NRz0mdl4M/urmV1RRIZE+sXUPWu/RZSwmbvmTquaPekmjewX3UCBn7RTETV9bwVuLWA5ItIGw8PD7c6CSEeqNUDTLV7pdA0FfWY2D3gR8M1isiMi3UI1FtJKkwcHM6e1oj1faUCX1ZkjT+A3Z84c7T/SFo3W9H0BeDewtYC8iIiI9A0FftJqdQd9ZvZi4D53v6FKutPMbJGZLdrgig2lu6k8j6aTVvdLl+mVK1e2OztNU0s72Lwaqe1LqNZPWqmRmr6nAy81syHg+8BzzOx7pYnc/Rx3n+/u8yebOgtLd+vn8qz2Sr0pXaZnzZrVlHXM2GtGU5bbarXsA7XuLwr8pBXqPmu5+3vdfZ67DwLHA1e5+0mF5UxEGtaKnrs6WUk7NaM9X7ueyqFaP2m2/qqqEBGRjjB9r+bUKhatNAjL26EjSVtPDbmCP2mWQoI+d1/o7i8uYlkiUrtaB2auV9YJTCcoSczdeXLhy6zUc7eSWp+7W69qj2Wrt2mEgj8pmmr6RHpUqwdl1glKGlXrwMytevRatdo+CIFfM2r9kvVr35IiKOgT6QD11mQ0W7mTWLUTl05O0q1mzpy57f9q7fqy9oM8tX4K/qRdFPSJdLjSW1SNDD3Rqto/nZik3bJu7dbbFKJcma438EvmVfAnraagT0TGSJ+M6qntA52YpHPUe6FUWttXdOBXaf48tH9JrRT0iXS5cjUXrajRy3uy0olJGpHV9KFSe75WdeBIVAr8VOsnnURBn4hU1eiTB3RiklapJ+CbNm3atv/T7frKySrHlfYF1fpJp1DQJ9In0rV/q1atqnl+PXJK+lG5Dh3NDvwaqfWT1jKzF5jZ7Wa22MzOLDN9NzO72sz+amb/MLOj25HPhII+kR5T1K3drGEp8qatRCcnaUTWrd1ytXzp9nz1duIoMvBrdq2fLqxax8zGA18BXgjsD5xgZvuXJHs/cLG7P4nw9LKvtjaXoynoE+kDndTGL6GTk6TVOkZfoyrtE9Vu8SbqCfxAtX495FBgsbvf5e6bgO8Dx5SkcWDH+P9UYHkL8zeGgj6RPlPp1m7eQWiz1HOSUvDXf/I+gq2W8Sur1fI1qtbn8RYV+OVZVhbtV4XYycwWpV6npabNBZak3i+Nn6WdDZxkZkuBnwFvbmpuq5jQzpWLyFi1PGUgz+2qZtTyzZ07l2XLlmVOnzdvHkuXLq1pmXPmzGH58rZeBEsXqPcpHLXe2p05c+aYC6TZs2czPDw86rNK5bbafpAEfpX2pbzLytIP+9X4SQPNfJbz/e4+P2OalfnMS96fAJzr7p81s6cC55vZge6+tdBc5qSaPpEuUq3motWPXquk3gfNi9Sq2bV81VQqt430yC1qWdqvmmYpkC5o8xh7+/Z1wMUA7v5HYDtgp5bkrgwFfSI9rFoQWFprUU6jg89K79pnSn0dI0R6xPXA3mb2ODObSOiocXlJmn8BRwKY2X6EoG9lS3OZoqBPpI/UM1SLiHQ/1fYVz923AKcDvwRuJfTSvdnMPmxmL43JzgD+w8z+DlwEnOrupbeAW0Zt+kREpGlm7DWjsGWVa8/X7lu7RajWRlY6l7v/jNBBI/3ZB1P/3wI8vdX5yqKaPhER6Ti19Nytpt7x+TqRntohjVDQJ9LBWv0M0U6gE5M0Wyd1eBJpJQV9Ij2i9ETWCSe2InsuSu8qcmDmdt3abVUP3kbpoqq/KegT6VLNuGXV6+N5SW/px5pwkUYo6BPpE+q5KyLS3xT0iYhIx6v3SRxFq/VxbHm14rFsCd3i7V8K+kREpGf1Us9dkUbVHfSZ2a5mdrWZ3WpmN5vZW4vMmIiI9J4mPiNVRKpopKZvC3CGu+8HHA78l5ntX0y2RKTZ8jyCTaQd6h2jr9sGZW4n3eLtT3UHfe6+wt3/Ev9/kPAIEj2MU6RJ+umEphOSiEjxCmnTZ2aDwJOAP5eZdpqZLTKzRRt8axGrE2mbRsqzbmtJJ0qX6ZUr2/Yc+JppuBaR2jUc9JnZFOBHwNvc/YHS6e5+jrvPd/f5k039RqS7FVWeixyMVqQR6TI9a5YuTLpFJw34LN2joSjMzAYIAd8F7n5pMVkSkUZ1wtM4RGREJwZpakbRfxrpvWvAt4Bb3f1zxWVJpL8U+WB5EalMF0TSzxqp6Xs6cDLwHDP7W3wdXVC+RKQGOpFJL2v2wMydsv/UMkCzSD0m1Duju/8BsALzIiJNokewSSvM3XlyzfM0s32rBmYWGU09K0SaaMZeM9qdBREREUBBn4iISN9SZ47+oqBPRES6Xj8NXi5SLwV9IiIiXagTh4GRzqagT6SD5O2lqAbqIiJSq7p774qIiIi02/hJA3rKUU6q6RMRERHpAwr6RKSipUuXNjS/2h1JLfSEGpHmUdAnIiJdZdzswXZnoSwNfyKdTkGfiIiISB9Q0CciItLHVEPZPxT0iYiIiPSBjg36li9f3u4siIhIA/TsaZHO0rFBn4iIiIgUR0GfiIiISB9Q0CciIiLSBxT0iYiIiPQBBX0iIiIifUBBn0iH6tSnDohI88ydO7fdWZAe1vagb3h4uN1ZEBEREel5bQ/6REREJFi2bFm7syA9rKGgz8xeYGa3m9liMzuzqEyJSHO1soZ96dKlLVuXiIhkqzvoM7PxwFeAFwL7AyeY2f6NZEYnBxEREZHmaKSm71Bgsbvf5e6bgO8Dx+SZcfLg4Kj3qs4WERERaa5Ggr65wJLU+6XxMxERERHpMBMamNfKfOZjEpmdBpwW326cf8Wvb+KKXzew2rrsBNzf6pW2ed39tl6Axzd7BaXl+ROb7rpp28TVqYTJ/7c1LSv99vv24z7c9PIMY8r0ejO7vULy0dtidcnU0vfNK/9p7SwblShftdnJzCrla/eW5aSHNRL0LQV2Tb2fBywvTeTu5wDnAJjZInef38A669Ku9bZz3f223mTdzV5HJ5Tndq6739bbznW3ojzD6DJdTTt/hyydmCdQvmrVqfnqNY3c3r0e2NvMHmdmE4HjgcuLyZaIiIiIFKnumj5332JmpwO/BMYD33b3mwvLmYiIiIgUppHbu7j7z4Cf1TBLrlsITdCu9bZz3f223nasu5++a7+ut53rbud3zqI85ad81aZT89VTzH1M3wsRERGRrnDwbo/137z75KYse6c3f+aGXmprqMewiYiIiPSBhm7vioiIiLTTuEkTxzz0oReZ2e7Aa4BnAHMIw+StAK4h9Ku4p9oyFPSJiIiIdDAzey9wNjDA6DGR9wWeA5xpZh9y949XWo6CPhEREZEOZWbvBz4MbAQuBH5DGCvZCGMkPwc4DviomU1w9w9nLUtBn4iIiEgHMrPHA2cBQ8DR7l7uOTffNrOPEUZT+aCZ/cDdyz5ZRx05RERERDrTmwhjIb88I+ADwN1vBV5OqP07PSudgj4RERGRzvRs4Fp3v6FaQndfBPwJWJCVpueCPjPbzczWm9n4duelF5jZkJkd1e589AuV39Yxs1PN7A/tzkc/U3lvHTNbYGZL250PqdluQC3P4b4e2DVrYtcHfaVBibv/y92nuPujLcxDU3cmM/t5PDAmr01mdmNq+pCZPZyafmXJ/G83s3vNbJ2ZfdvMJjUrr1KbPim/zzazq2P5GyozfTBO32Bmt5VeZFQqv9Xmlc7S7+XdzHY2s4vMbHmcfo2ZHVaSt60lx/tTUtNnmNmPzewhM7vHzF7VrO8hHWM7YH0N6dcDj8ma2PVBXz9w9xfGA+MUd58CXAv8sCTZS1Jpnpd8aGbPB84EjgQGgT2AD7Uo6yIADwHfBt6VMf0i4K/ATOB9wCVmNgtyld/MeUXapFJ5n0KoiTkEmAGcB/yfmU1JpVmePt67+3mpaV8BNgGzgROBr5nZAc34EtIxVhF+77xmA/dnTezqoM/MzidUff40XhG9O175u5lNiGkWmtlHzezamOanZjbTzC4wswfM7HozG0wtc18z+5WZrTaz283sFalpR5vZLWb2oJktM7N3mtn2wM+BOakrszlmNs7MzjSzO81slZldbGYz4nKSPJ4Wr/hWmNkZOb/zIHAEcH7OzXQK8C13v9nd1wAfAU6tsPyT4xXkKjN7X8m0Q83sj2a2Nub5y2Y2MU77ipl9tiT9T83sbfH/98Rt9mDcrkfmzH/P6pfy6+7Xufv5wF1ltsE+wJOBs9z9YXf/EXAj8LKYJLP85pi3dF0zzezyuN2uA/Ysmf5FM1sSp99gZkfEzx9roSZxZirtIWa20swGzGwvM/uthZqb+83sB1nbop+pvIO73+Xun3P3Fe7+qLufA0wEHp9j+21PKNsfcPf17v4H4HKg7PPHzOwxZnauma0xs1uAp5RMT77vg3E7/Vv8fFLcnk9Ipd3Zwt2kWWa2k5ldYeE8sNrMfm9mXR1LdLi/Ac+sIf0CwoVwWV39Q7n7ycC/GKnl+nRG0uMJO8ZcwoH+j8B3CFdatxK6Qyc71a8I4+DsDJwAfNVGrqS+BbzB3XcADgSucveHgBcy+upsOfAW4FjgWYSRs9cQrtLSng3sDTyPMLBinltTrwZ+7+53l3x+QTwJXWlmB6U+PwD4e+r934HZ6RNYwsz2B75G2FZzCLUn81JJHgXeDuwEPJVQ+/KmOO084IRk5zezneL0iyx0OT8deErcds8ndD/va31afksdANzl7g+mPvt7/DyZnlV+q81b6ivAI8AuwGvjK+164GDCdr0Q+KGZbefu9wILgVek0p4EfN/dNxMC0SuB6YT95X+rfOe+pPI+lpkdTAj6Fqc+3tnMhs3sbjP7fPyeAPsAj7r7Ham0lcr7WYTttyfhmHtKyfQ7CRUIUwm1598zs13cfSPwfUIZT5wA/NrdVwJnEMaIm0WoVfpvRg8WLMX6IbC3mb2gWkIzO5pwN+TSrDRdHfTV4Dvufqe7ryNc5d3p7r929y2EDfqkmO7FwJC7f8fdt7j7X4AfEQY9BNgM7G9mO7r7mjg9yxuA97n70rgTnQ0cl1zRRh9y94fc/UbCQe2EHN/l1cC5JZ+dSLj1tTtwNfBLM5sWp00B1qXSJv/vUGbZxwFXuPvvYp4/AGxNJrr7De7+p7hthoD/RzhI4u7XxWUnNXjHAwvdfZgQLE4ibLsBdx9y9ztzfFcJeqn8liotn8T3O2RMT5ffavNuY6GjwMuAD8Y830S4UNnG3b/n7qvitvssocwmNTDnEU+CcVknMFLbvpmw781x90diDYzUr5fL+zZmtiOhDH0ofleA2wgXHrsQBtw9BPhcnJa7vEevAD7m7qvdfQnwpfREd/+huy93963u/gPgn8ChcfJ5wKtSNXgnM7q87wLs7u6b3f337q6gr3kuAK4D3mXxzlo5ZjZACMhviPOU1S9B33Dq/4fLvE/aU+wOHBarrdea2VpCQPXYOP1lwNHAPfF2zlMrrHN34Mep5dxKCH7S9+aXpP6/h3CFmcnMnhHzckn6c3e/Jt7e2uDunwDWEq7gIDTq3DGVPPk/XTuSmJPOU7wqXpVa/z6xWv9eM3sA+Dih1i+x7cQY/54fl7MYeBvhQHqfmX3fzCp+VxmlJ8pvhtLySXz/YMb0dPmtNm/aLMJg9KV53sbMzjCzWy3cpl1LqAFJyvdlhABiD+C5wLp4oQPwbsLYWNeZ2c1mVlqDKLXp5fIOhFuvwE+BP8VjNgDufq+73xIDsbsJZSsJYmsp71ByPGdseX+1mf0t9Z0PJJZ3d/8zoW3is8xsX2Avwq1kgP8h1ExeaWZ3mdmZtXx3qU0MrA939yPdfVOVdEe6+6HxwqWsXgj6irzCWAL81t2npV5T3P2NAO5+vbsfQ7iV8BPg4gp5WAK8sGRZ27n7slSadLfq3YDlVfJ3CnCpu1fryeOEkxDAzUD6du9BwLC7rxozV3hw87Y8mdlkwi3exNcIV6J7u/uOhGp9S03/HnBMvL28H2EbhQy5X+juzyAcXB34VJXv0C/6qfyWczOwh5mlaysOip8n07PKb7V501YCW8rkGQAL7ffeQ6gdme7u0wi1KAbg7o8QtteJjK71SE7U/+Hucwg1Rl81s73yff2+0+/lHQu9z38CLCOUl0rSx/I7gAlmtndqelZ5h5LjOaPL++7ANwjNbmbG8n4To4/nyUX8ycAlcR/A3R909zPcfQ/gJcA7TG20u0YvBH3DhHvYRbgC2MdCZ4aB+HqKme1nZhPN7EQzm+qhHc8DhCvBJA8zzWxqallfBz4Wdy5iA9hjStb3ATObHNugvAbIbAAerwxfTsmtXQvjXD095m87M3sX4Wrtmpjku8DrzGx/M5sOvL90GSmXAC82s2fEauQPM7qM7BC/9/p49ffG9MzuvpTQLup84Efu/nDM4+PN7DnxYPcI4Wq9ZUM0dLieL78WGslvR3hQuMVyOhHAQ/ukvwFnxc//DXgi4TYdVCi/OebdxsOQIJcCZ8c878/oNk47EILClYQT6wcZW6vyXUInkpcSLnCS7/dyM0vavq4hnKhVvsvr6/Ju4RbcJYRj4KvdfWvJvAviMd3MbFfgk4Ra5uTOy6XAh81sezN7OnAM2Z36Lgbea2bTY/l8c2ra9oRyujKu9zWEmr6084F/IwR+303l8cUWOi8ZI9tV5b3JzOypZvYNM1tkZv+Mf79hZofXspxeCPo+Abw/VlG/s5EFeWgQ/jxCe7TlwL2EGqlkXLCTgSELtzb/k3gr08OjUS4C7or5mAN8kVAdfqWZPUgYJfswRvstoZr8N8Bn3P1Ksh1LqHm4uuTzHQg1cGsIV44vIFyxrop5+wXw6TjfPfF1Vsb3vxn4L0LD6BVxmenxrN4JvIpwO+EblD/onQc8gdEHokmEg9f9hG26M6GWUPqj/D6TcJL7GaG24WFCx4fE8cB8Qnn7JHCchwbjecpv5rxlnE64NXgvIXD8TmraLwntx+6I63iE0bfGcPdrCG1c/+KhTWviKcCfzWw9YZu91cd2tJKg38v70whtEZ8HrLWRHsRJc5wnEzquPEQYmusmQieTxJsIY7DdF7/DG+Nxu5wPEcry3XH96drpW4DPxnUNE47Z16RnjhfxfyEEh79PTdob+DXhdvMfga+6+8KMPEgBzOyThN/ntYQ2rXvEv68FrjWzj+deltpftp6FIQfuBgY8NE7uGWb2TEItyGDpVaz0hl4uv9WY2VXAhe7+zXbnRVqjz8v7twk9nd/f7rxU8uS9d/NrPv+epix78ktOv8Hd5zdl4TlYGIboIsKIFx8ljCRwL6EzzbMJdz92B45394vLL2XEhGoJRPKKty7eCnxTAZ/0GjN7CqEmpvS2n0jPicHuvzPSW1ra43TCbfjD3D096PJdhNrqywk1wqcz0m41U9XbuxYee3Sfmd2U+myGhQEx/xn/Tq/5a0hPMbP9CL2GdwG+0ObsiBTKzM4j3NJ6m48eF1Ck55jZRwiBxP+oqULbHUToSFP2KRuxOcslhPbMVeVp03cuoZ1Y2pnAb9x9b0L7BnXZroGHceqsl24VuPut7r69uz/N3R9od36keXqx/Fbj7qe4+1R3P7fdeZHW6tPy/oHYE/pj7c6LMEBo41nJesIg31VVDfrc/XfA6pKPj2FkYNPzCJ0MRERERKQ4dwIvtIxH3cXPX8jop7pkqrdN32x3XwHg7ivMbOeshGZ2GnAawAB2yEwbGDV9u3Hl486BgbGfT5g0fsxn4ycN5PoMYNyksYGwDUwa+9mEjIC5zOee8cjBLVvKXxQ++mh2z/aseapN27Qpc7zGbTZv3lw1Td5lNZK+VN585bFx48b73X1WYQsso1p5TtRSrqF82YbsslxLGYcayzmULetQe3mH7DJfaZ4804sq960o8/WU81aUZxhdprfffvtD9t1337Lplv3lxorLqbXMQ+3lvtK0Wsr+tmk1HOsTRe4D1ebrxH2glrS1lvsDDywdNWbEDTfc0JL9oQNdAHwMuNTMzvDU06wsjNn4WcJwO7lGxGh6Rw4PD5Q+B2CXcZP8NQPzRk3fZ0r5nWvuzpPHfDZjrxljPpu+19gyMHXPuWWXOXlwcMxnA7uM/Wzc7LGfAeCyPtQAACAASURBVNjMXcd8tjHjqShr166t6XOAVavKjZccDA8PZ05bvrz6GKFLly6tmgZg2bJl1RPVsdwsefKe1x133HFP9VSNqVaeE7WUayhftqF8+YbayjiUL+dQW1mH2st7pWmVyjtULvNQXLlvRZmvp5y3ojzD6DI9f/58X7RoUdl0/z1pz4rLqbXMQ+3lHoor+1B7+Ydi9wGo/7gP7dkHain/tZb7rLIHYGYt2R860OcITexeQhhHdylhSLVdCM/6NkKP3s9lLSCt3nH6hs1sF4D49746lyMiIiIiZXh49NpRwPsIw7bsSnhG8q7x/fuA53qFR7Sl1Rv0Xc7IaPanEEcMFxEREZHiuPsWd/+ku+8FTAcGCY+L3Ct+nvuJKFVv75rZRcACYKdYrXgWYeT7i83sdcC/CI8HExEREZEmiaNj1D1CRtWgz91PyJikByyLiIhIW9nApIptN2WEnsghIiIi0oHMLO/g2Obug9USKegTERER6UxbCT10S00Dpsb/lwO5xsdR0CciPanR4YSKUuSwRCLSX9w9c6wkM3sc8CVCAPi8PMurt/euiIhIx9u8YqjdWRBpivhc5OMIY/Z9OM88CvpERHLqlNpDEREAd98I/Ap4RZ70CvpEREREutcWYHaehG0P+u5Y39izW9ttUr5BsEVEREQKZWazgGOBJXnSqyOHiIiISAcys7MyJk0gPIrtGEIv3jPzLE9Bn4h0HfWI7W1rFq9k+l6z2p0NKaE2rW3xwSrT1wFnufun8yxMQV8ZW4eHGDd7sN3Z6Hk6cYuIiFT07IzPHwXWALcX+uzdbrTuzmVM3XNuu7MhIiLSNqqZ637u/rsil9f2jhwi0rvWrl3b7ixIn9gwNNTuLIh0vJ6s6ZP6LFu2rN1ZEAFgeHi43VkQEekYZrYf8EbgcMITONYCfwa+6u635l2Ogj6pm24diIiINJeZvR74KjC+ZNIhwBvM7L/c/Rt5lqXbu31AwZlIeardbr5l921odxZEupaZHQZ8DXgAeAewH7Bj/Pv2+PlXzewpeZanmr4uoFtdIiLSDM268NHoDIV5J6Gn7hElt3HvAO4wsyuBv8Z0r6y2sK6v6VuzeGVD8+th3CK9R7XbvWvdnaqdlb5yBHBZVrs9d78NuBx4ep6FdX3Q16mmTZvW7iyISIEUSPYPX5XriVYNW7VqVd3zqiatb0wD/lklzT+BnfIsTEGfSBfKqu3QsBUiIj3lfmBmlTQ7A6vzLExBXw1adfUnkmi0+UJCzRikE61enOs8JXUoumZaNd1t8w9gfpU084Eb8yysr4I+1YIUp9EDgG5NiEirZF30bB0u/7lIB/km8ICZlX3MmJnNIYzZ9808C2so6DOzt5vZzWZ2k5ldZGbbNbI8aR8NXdFe7arx6MaTXlEXDK0o87q4kXqfSqNRGwTA3S919yPdvewBy92Xu/tz3P2HeZZXd9AXo863APPd/UDCoIHH17u8TtPqk2G9DXp1UhGRXlRP0wbdzamNLvYbZ2YvMLPbzWyxmZ2ZkeYVZnZLrCS7sNV5TGt0nL4JwGPMbDMwGWhqBLJ68Wpm7DWjmasQERlD7ZlEpJSZjQe+AjwXWApcb2aXu/stqTR7A+8Fnu7ua8xs5zrX9XjgdMJj2KYCq4A/ER7DVq137zZ11/TFqsbPAP8CVgDr3P3KepdXNI3lFDTjZKUToHQylU8RaZFDgcXufpe7bwK+DxxTkuY/gK+4+xoAd7+v1pWY2YmEDh1vAp4M7BnX/RbgJjN7Td5l1V3TZ2bTCV/ucYRGhD80s5Pc/Xsl6U4DTgPYsUMfALJ5xRADuwy2OxvSBRotz8vu28DcnScXnS2RuqXL9G677dbm3OSz7s5lTN2zbLv2vlCtWU+/XfjYhImMmz3YrMXvZGaLUu/Pcfdz4v9zgfSwHkuBw0rm3wfAzK4hNIM7291/kXflZnYA8C3gbsJj135HePTaZ4CFwBeAc8zsFnf/c7XlNdKR4yjgbndf6e6bgUuBp5Umcvdz3H2+u8+fbH3VWbhrtLpdRze3Q+yk8tzMsfraOTxRKxqwqy3TiHSZnjVrVruzIx2o34LIEvcn+0d8nZOaZmXSe8n7CcDewALgBOCbZlbL0xvOIMRqL3H3n7v7Q/HzDe7+c+BZwCPAu/IsrJGz1r+Aw81sspkZcCRQ9jEhItJ+GqtPRBK68CnEUmDX1Pt5jO3bsJTwGLXN7n43cDshCMzr2cAvstrtufu9tOIxbLEa8RLgL4RBAccB51ScqQOod1dj+vyKry2KGqC5kzTy+KlWU5nvLlnHeF30tE43382p0fXA3mb2ODObSBjB5PKSND8hBG6Y2U6E27131bCOxwJ3VEmzBMjVy7Wh+1Pufpa77+vuB7r7ye6+sZHl1atZJ8Vyw7aUu+01yTc1Zf0iUptOCtD66MTXNL14wSO9w923EHrU/pJwp/Nid7/ZzD5sZi+NyX4JrDKzW4CrgXe5ey1XvQ8BA1XS7A3kahfTmT0rpGfpRNh5tg4PNbMRdKGKKD+6rdVZ6h2Kq5c7czTSrjXPhU8t+0AnXUh1Inf/GfCzks8+mPrfgXfEVz2WA5k9rMxsAfBSQmePqtSzIurUqv9KO38RPbh0AuxNasZQDJ3wirHsvg3tzkKmbnwqjfSVa4AjzKy0tu+5ZnYp8GvgPuBDeRbWsUFfEQcJjdVXLJ0Am6ueR7HVWsY79eJGpN/oKUyS0w8IQV16KBgHngq8GLgMOMzdV+RZWF/e3t0wNMTkwcF2Z0NEWkw12/1Bx/hs2ge6i7svBPYv+fjZwHrgjtQQLrl0bE1fp1DVf1BELZ+uUKWZiq6JVs12Z6jUmUM13aNpH+gP7v57d/9rrQEfKOgbpdcPCKV0xdc8d6wvtke3Hj7fuFaWd13gtJ/Kv8hYXXd7t96eXlI/Xe11tlp7MXbyYwdb8TQOkX7VigsfXfAUy8yuzpvU3RdUS9R1QV+WNYtXMn2vsY8QyjohNtLmw1ctwWbuWj1hG3VaoKYDQf9Zu3ZtoctrZRnqtP2nG9yxfhP7TJlY17yNXMz38tAtIsARlH/cW8IY++i3TLq9W6LcLd5mt+urpxdXq06AOvl1t7y3uHq97aqaMvSnVt3izbrAaeWFT7Vjda37gI79ncHdJ7j7+NIXMB04CriB8HS0agM4Awr6yuqUtn3NvNWV5wBQ1E6vWr7RmjFmWaUG7eVOfEWW8VpPbM14BFuRJ6hGl6XyXrxqbVpr6dDRCcf3SvuAmjhIHu7+gLtfTQj8ngGcmWe+vg76mnkVWO5EWO6zrJ2/mYMytzLgk9pUGquvmx5JlRUI1nuya1UNh8p97ynyOL/R6rt9nVcn7AN5liedw93XAT8HTs2Tvi+CvnoGaS69Gmz09lcttSG9FPCp1qMzNLu2r1GdcrKT7tVIbV8jx/dab+0WXdOtfUCAB6jwqLa0vgj6Kql0FVjppOirljS87lp3/monv04L+KS1qp30qtV41HPiy1ujDc2p1a5EtRudpVqzhmpPpMlT0521D3TD8C317AfNCPga2Q90kd96ZvYY4GjCUzuq6pneu63QigfTZ+34jdR2gK72JKil13qRvdSb0VmpUrlvV8Cnk55AsbV89ZYpHfN7g5mdkjFpArArcCKwJ/DpPMvrqaAva9iWaiqdCOsZ06x0h8/blq/dAV/RNR06ARavWhlv5ph9edoz1dqpo95G650Y8ElnyDtMV63H9nLlv6jeufXsB0XuA3mWKW3zbSoPy7IV+A7wgTwL66mgr5JOGsupHwI+6X15T3hF39ZVwNe/8l7Yt/N434rmDc0I+KRjvSbj80eBNcAid891axe6NOhrxlM52vmA7mYEfO0O9lTLV79Gy/f/b+/e4+0o63uPf74QQMMtECCQBA2CqEhVbPCGFxRs1VrQVj1wBNF6Dq2tl1pbS2sVpbW2trXSc6ieHFHxUq0CnlLrtQpFraUEpCpEEDFCEggYCBcRksjv/DGzZWVnXWbNfdZ836/XfpG9ZvbMszfPWs9vfs9tUoM3rq7XMYQBHPBZ9dr+oA/tCvj8fminiPhImdfr/USOaU072H3+U1/WMR0O+GyUspduKWsWb9FhDNDtgM/1vrhJkzmgWP0vc0LHNJOYptXFutTFMs8KSY8bM/ZvOw76BjQxw2tYAzjuCW9Sw9d0wGftkGeZomlkadymHcZQR8CXZZa7zY6q3wdZTJvlq+vBJ+t1rRNeRDKub6JWB31l71yQ9wNgVCZk2mVb5r/5pw34xpkm2KvyDe6nvfaos8HLktHO28iVGfDZbCgr2131zjQO+KwoSftL+gdJGyVtk/TA/C/g7UDMe22oTo7pGyfvDN45Rcf2FUnvVz2Go+o3twO+7mhyDOugWdlpw3W/+8aNZ50/c7fI53ybAj7rhL8F/hvwI+B7JBM45luRfl0y6WIzF/SVoawGcfCDYVKWr8qAr44nOTd601l/670sO2Bh7p+vexbjuOVa8mb5ZiXgs+yK1vsyVPHAU8Z41lGqDPj83uiE44B/iYgTRp0g6W3A2yPiOZMuVijok7QI+ABwJMkaMr8REd8scs2qtWlG1yQO+Kwu06xZVvVafA74+i3r7PUuLN+SVROfoTO13eaCXUtbSL6F9gOuKOtiRcf0nQ18ISIeDTweWFO8SLMvS5bPAV+/ZZnFmFUbBrPPV0e9qSvg83vA5qtjLVZ36/bGTsC2Mi+Wi6S9gGcC5wJExJaIyDXQ4bp7tuQtRmXGzeQtsjl3UQ74rE3yLkk0StuXIzLLsuNSGaqsw35/dMo7gEsnnHMJyWSOiYpk+h4B3AZ8SNK3JH1A0u7zT5J0uqTVklbfGyMnlLRekVle4xrCabJ8WRq8upak6GvA16b6PM0MxjzZvjK6S5qaoW7ZDdbp224rdw3IqpW9ZmVV8mT5xvFs3f6IiLMi4msTzrk0Is7Kcr0iQd8C4InA+yLiKOAnwBlDCrMqIlZGxMqFKm+FmHHdX5M+CLI2gFWs25d3r9Es3J1Vvarqc1/MWsA3C++FwTq9//75Vz6wRNFsN1RTj71GZTdJenjWryzXKzKRYx2wLiIuS78/nyFBX1/lTfnnzfI54LNJ5g9oL3sWY5Z1KPPw2CUrWxuWLKr6AWiaa1qr3QAo47kTMxG5g76IuEXSTZIeFRHXkkwrvibv9fooawDVhoDPwZ4V1fV1KK0+0+w/nWUW77QzeLMObcjycF9l704WXqOy8z7MjkGfgANJJtAuAb4BXJ/lYkXX6Xsd8HFJu5JEo68qeL1ea2uj5Te1DVPVAPY82vre6Zvr7tnC4XuMXtOxi8atUTlfGXurj9OGPaatXhHx6lHHJO0CvA34HWDkeYMKBX0RcRWwssg1bLwms3wO9tqv6A40Vcma3ejykhR+f1jZHKTZNCJiK/BWSccDf0ayc8dYvd2Ro4kFO8tO85f9AeFGzLrIDaV1/eEnC2f5bIx/B07JcmJvg74uqCPT4UCvvaYZ11SXvF2609aztmf5rJ/KGNLgz1yrwENJdu6YyGtOVGzUGI/5b/xpn8rKeIrzh0//VLU7R55lKsqow85mdM/6W+9tugitVGZd9vuid1YBz85yojN9PeWAzyw/v39svjLW5xvFme9+k7QTcDSwArg+IrbbizedX5GJM30dVPQpzg1WPdq4vWAXtHUhZjOzuklaBPwb8E3gE8Dlki6UtCA9frCk90l6YpbrzWzQ15XtecxmTdPrklk/+TM/UefDkBMItfgz4GnABSTL5H0JeBHwW+nxden3p2W52MwGfV1XVTrfb1Izs+r54cdK8gLgPyLipRFxDvArwA9Il2eJiCDJBD4ry8Uc9JmZWeuM21+9q/JuvWa9diDw9blvIuJnwKXAEQPn3ARk2kbGQV+POMs3m7rYrVW0kXMjaWY9sZ5kq7VBtwB7D3y/M7Awy8Uc9JlZpfywYWaW2yeB50s6YOC1+9k+fnsiSbZvIgd9ZlappUuX7vDa8uXLC12z6M+bmXXEO4HvAZ+T9BRJmjugxKuAZwIXZbmY1+nrkaVLlzrrMoPauP2Ume1o1Gfw8uXLPWTBRlkD7AIsJdlu7afA3HpgdwB7AdeQzPKdyJm+llq2rN59gc3M2qRtWxBOa8mS+cOwzHJ5gCTIW5t+bSQJ9q4FVgN/CDw5IjLtETizmT5nP4Zzts+qtmTJEi9XYbXzZ36izqyh25PqRcShZV7Pmb4OKjqeadgYKyvf4Xvs2nQROslZbjOzajjoM7NOanIyhx+cbL7FixdXdm0/CNk4kl4h6eIs585s925bLF68eOhG3EXT4kVT+HONllPz/bL3oe1pPDx43cxssnT/3acCi0nW5JvvROBZkl4JBEBEnDfsWg76WmJYA7hs2bLKtmOb4zEZ7dXGgeyLFi1i8+YdxwuPeriZM209y1r3HTh2z7IDMq0hO7PqmMHbt/dFaCfu12wOp5H0S8CFPLj4cow4NYBz534McNA3qE0Zj7zKemMP66pyIGhd0bcGznZU1ySOUQ89o5Q5qamOJIC10pnAVuAdwI3Az4ac86L065WTLtbboK8JeT4Amnqju/u3G9o6YzFrXR8XsLW9kXOW3NrMD0Mz4xeAD0fEO0adIGkF8KKI+Miki3kiR4vkGZhe9WD2pUuXetC6DbVo0aKmi/Bz3qGjHTxjvZiy67HfFzNhd5K1+UpROOiTtLOkb0n6bBkFsunV8cZ24GdF5alDnrVoZZh2OE9senAb091iy8jz5j/4ZJ3Bm/fzNM/7wYFfLwSjx/ptp4xM3xtItgmxKWV942d5o9cV+Dn4mx0LV6yo9PpZdySYVHezNnRu3GbHNJOY8gxxqLruD5p2Z44q6nGV7w23CZU7BPi7CeecnZ43UaGgT9Jy4FeADxS5Th5FZjbWMYljVNfXpA+AUW/OtgR+4Dd5V1Vd77NkOUbVnbLqbt2Bn98Llte4ujOuHufNfvuhqJsi4saIuHvCOXdFxI1Zrld0Isd7gTcDe446QdLpwOkAe9U0b6Sswe2DT4O7HLRi5HmTjFvOYthg8FEDcLMMbK9r8G5fB7E3UZ9Hmaaed322+jSTOjyAfTqDdfphD3tYw6XpvmGf96MmNtX9Oer3RvdIOjPrqRHxdkn7Aq+b+37+SbkzfZJeCNwaEVeMOy8iVkXEyohYuVCzMW9kpyUrGru3M37NmsX6PMrguKZRJk3mGJbZzpvtmybD4axGdoN1ev/92zkbfJQmZq/XNYGpimzfpOtaK71tii+A/eZ9v50irdYxwAmS1gKfBJ4j6WN5LtTGGV9VjvnIMsaj6Btz+fLlnuDRYWUuzFxnlq/oQPYuBn5+D/Rb1iBw1Od+nZM65jjw65RnZ/w6Nj3/xnnfbyd30BcRfxQRyyNiBXAS8NWIOCXv9erSxm6uad/0bWv43Oj1y7SZjrIHsret/lu5sj7wtHWNSihvH94q62+Z13YbUJ2IuDTrV3r+fYPfzzfb/VM51Tmza86wN82k9H6bZjX6TV+vrA1e1oecSWNWxy1bkVfegezgwM+ya3Li3nx5sn1VdfNaP5US9EXEJRHxwjKuNUsGPwjmP/llHes0qZt2msCvjoWcLZum9x8t88Emy1plZTd20K7Az3U/m6brPZRT97MEeaOyfW0K/PxA1B2SHi7pLZLOl/Tl9L9vkTTV7KuZy/QVTfkX/UCYpuurrEHubcv6WTs0OZShj4Gf1aesrt3BDPe4CXpZMt3DPvun7eYt8l7Iy++L9pP0GuA64E+BFwPPSf97FnCdpN/Meq1WB31lPxVOagRHBXyjur60+OCp7t+WQe5eqNOKmt/AZW3wHPhZndo4hhvGj3P1zjU2SNLxwP8GfkIS9B0PHAk8F/hz4D7g7yUdl+V6rQ762mja5VomdX1N2whmaQCzLutSVQPowK9aWbIc0zR2RdagzKMNgZ+Dv3bKMomjSJavymENo16btpsX8vf45FXG+8Gf+5X5feAe4OiIODMiLo6INRHx1Yh4K3B0evwPslysk0FfmctZzGli8sacPIFfF7J+lk/R+p03ow3F1qAso8GbFPiV2eiVXffd6LVLm7J8fQn8rBJHA5+KiB8MOxgR3wc+nZ43USeDvjza/gFQVfajqcWc3QDanDIDPyi30XNDN5ua/Lyvekkj8IzenlkI3DbhnFvT8ybqTdA3zrjMR56uryzjnYapIvCDdu3iYbMja/dWHkUeGhz4za4yu3an/WwfNpljmvqeZ+2+uh+e/V5opZuAZ00451iSRZknmqmgr+rFOsvcfq3uWV1NBH7O9pVvUh2ftmt3XMM3f6JSkbX6yuzegnLH+GW5Xlau881rQ6/OqGCwC9281joXAE+V9L8k7T14QNIiSecATwEuzHKxmQr68ig7yzeNabp5oZuBn7VHFeNWy8j2tSXws2aVMVZ7VMDX5JjtrNrSzVukDfBDTyXeBVwD/A7wI0mXSPqEpEuAHwGvAa5Oz5uoF0FfFU9+0y7XUnT9vjldC/z8IVCfaev5/IeaMjPZw4zLbldZ57PyA0+1ii7BVWZPTpkP9NM+5JTdzeuM32yLiLuApwOrgF2AZwAvS/+7IH39mPS8iXoR9I0yTZav6FIto9Td7QXO+PVR3ZmOabu3ivDEDoNiWb6sn+/TDHHIM7617PX7IF/g5/dBu0TE5oj4LWARcBTJGL4nAIsi4reyBnzQ86BvlKq6deuY1QXtC/yc7dveuIzHuC6ucZmOcVm+YY1ek5mOUYo0eG0K/Fzf69eGcXzTKPvhx0Md+iEitkbEtyPiaxHxnYjYOu01HPTNM6wxrLrbK4+6GhY/8XVb1gxfnXW8imyfzba8XbtVP/BMUna2rwgHfgZJf/BMGPWhUMXA3mnH802yePFiNm3aVOo1ly9fzrp160q9ptk0Fi1axObNm5suBpA0eOvXr2+6GDPt8D12beS+bcjy1VXXly5dyoYNG4YeK/sz321IO0i6OOupEXHspJM6F/RVsRuHjVfGm3/ch5UVU3QCR1fUWYfc4NXLn+tmIz0D0JjjAiLrxToX9Fl5nP2oTlNZj65bsmQJGzdubLoYNgO6sExL3er4zG/iAX/btm2t6VUoW0QMjdMk7QX8IvBu4Abg5CzX85i+AV0Zz9cEj+1rVp4xTX1r9Nq0YLMnc5Qn7wSmYbqa5c6q7M9pf+63V0TcFREXA8eTLOlyRpaf62XQ17fGsC3cEJqZNcsTOmZLRNwJfB54ZZbzWxv0FV3IE9oxwNcsqzzjmpzpmI4bvH7wg731zF3Aw7Kc2Nqgr25dbAyLrF2Wh1P93eFGrxyu8+Uo4yG+Kk0N4cm7dJF7TGyQpIcCLwBuzXK+J3K0XNUD2+uezOFZvO3jcavWRbPck1Pkcz/LzPNpPvc9k71Zkk4bcWgBcDDwcuBQkgkdE81E0FfmnoyDhjWGZa/RZ2bZTXpoKLvBs+rl3YXGrCc+yPhlWR4APgS8NcvFcgd9kg4GPgIcmN50VUScnfd6dami2+t+9Wd5Dj/11c8NXzc5q92MUZ/xXRzC01Wu+6V61YjXfwbcAayOiExdu1As07cNeFNEXClpT+AKSV+OiGsKXNPMeqCKXWiq4gcdM2tKRHykzOvlnsgRETdHxJXpv+8G1gCzO8hihtU9o9EDkcsx7RaDXcl0VLX36CDP4rW+83ugn0qZvStpBXAUcFkZ16tbVxpD6742z2Ic5LGr1mazPIkjC6/cYHkVDvok7QFcAPxuRNw15PjpklZLWn1vPFD0dpn1/UOhSn3+AGiqPptVZbBO33bbbU0Xx8wqVCjok7QLScD38Yi4cNg5EbEqIlZGxMqFKhZjelPu6fU5QKtCmfXZqlFFnZ/l99Fgnd5/f08aMptluVstSQLOBdZExHvKK1I71L12mRfqNDN7kPebNitfkVTFMcCpwHMkXZV+vaCkcpn1SlNZ7C4uzFzWg44HslsdFi1a1HQRRvJ7oH+KzN79ekQoIh4XEU9Ivz5XZuHK5qfA0aZ98xft7nKG0vrA9bwdujhZr45Z7NY/HpQ0Bc9otLqVtTBzFxs9m31Nj9PuYqY7K8/gtWFmYhs2s76Zdo0+MzPrLkmPBY4FHgccADwU2ATcCHwduDRdM3ksB30V2bx5c9NFMDOzHOrq1enSzjRWv3S3s1cBrwEelb48bA/ePwS2SLoAOCci/n3UNTvfvet9Sc1sPi/bMru8BmvC40Vnl6SdJL0G+AHwt8BewMeA1wLPBZ4IHEEyofYlwLtINsd4GfB1SZ+VdPiwa/c+09eFsU5Llixh48aNTRfDzMwsk6VLl7Jhw4ami9FVbwbeCVwEvB/4UkQMy/DN+QyApGXAKcDrge8Au80/sfOZPjOzvLxkRfW6svVgX/k90Er/BRwdES+OiC9OCPh+LiLWR8RfAocCfzzsHAd9PVBVt5S7u8zMrM8kPU/StZKul3TGmPNeIikkrZx0zYj4fERcmbdMEXFfRPzNsGMO+uzn/MRnbeE1ysys7STtDJwDPJ9kjN3Jko4Yct6eJF2ul9Vbwh056DMzs87zckXWgCcB10fEDRGxBfgkcOKQ8/4UeDdwX52FG2YmJ3J4dpdZOyxatMjLF5k1ZPny5axbt67pYsyyZcBNA9+vA548eIKko4CDI+Kzkn5/2htIujjrqRFx7KSTZjLoK2qWV2lvE8/uMjOzon72s59V+XC5n6TVA9+viohV6b815PyfT7qQtBPJkiuvLHD/Z4y4z89vw/C1+4Zy0GdmneMHhtnm9VetRX4cEaMmX6wDBlfyXg4MfjDtCRwJXCIJ4EDgIkknRMRgIDlSRAyN09Jxgr8I/CVJtvGkLNfzmD4zMzObyCs27OBy4JGSDpG0K0ngddHcwYi4MyL2i4gVEbEC+A8gc8A3TkTcHRGXAMeTjC0cukTLfI0HfYfvsWvTRSjkfnW7/GZmZja9iNhGskvGF4E1wKci4mpJZ0k6oaYyXIaA8QAAIABJREFU3A18DnhFlvPdvWtmZmaWQ0R8jiToGnztbSPOPbaiYmwlmVQyUeOZPjOzrnD3lpm1iaRDSPbfXZvlfGf6zMzMzFpI0odGHFoALAWenv77rVmu56DPzMxm1i4HrWi6CGZFnMr4JVvuBN4YEedluZi7d81sJrWlK3bp0qVNF8HMuusRwCFDvp4A/AZwB/D4rBdzps/MzMyshSLixjGHvyPpi+l/r4uI90+6njN9ZmZmZh0UEbcAnwV+O8v5DvrMzMzMuusO4LAsJxYK+iQ9T9K1kq6XdEaRa5mZmZlZdpL2Bk4ANmU5P/eYPkk7A+cAzyXZf+5ySRdFxDV5r2lmZmZmCUlnjji0gGRB5l8FFgNnZblekYkcTwKuj4gb0oJ9EjgRcNBnZmZmVtzQ3T0G/Bj4w4j4qywXKxL0LQNuGvh+HfDk+SdJOh04Pf32/ndtueG7251w+4irD3v9e9MXMrUfyR+mCU3du2/3BXhU1TeYWJ/nTFOvIU/d7tv/3z6+hyuvz7BDnb5H0rVjTh/9t5i2zkORz/RBTdaNcVyu6ewnaVy5Hl5bSdrl2SNe30oSd62PiAeyXqxI0DdsscDY4YWIVcAqAEmrI2JlgXvm0tR9m7x33+47d++q79GG+tzkvft23ybvXUd9hu3r9CRN/n8YpY1lApdrWm0tV9Mi4tIyr1dkIsc64OCB75cDG4oVx8zMzMyqUCTTdznwyHSz3/XAScB/L6VUZmZmZgaApGcCrwCOAvYi2X7tKuDDEfH1rNfJHfRFxDZJrwW+COwMfDAirp7wY5m6ECrQ1H2bvHff7tvEvfv0u/b1vk3eu8nfeRSXKTuXazptLVetJD2HZFWU90bExvS195OMu50/hO4o4DckvT8iMi3OrIgdhuGZmZmZdcKRRx4Z559/fiXXfsxjHnNFnWMNJf0L8FjgkIgISW8A3gNcB7wT+DqwEVgCHAO8BXg08MaIOHvS9b33rpmZmVk7PAH4cjyYkftNkvkST46IuwbOWwuslXQRyVJ5pwMTgz5vw2ZmZmbWDvsCNw98fwhwwbyA7+ci4m7gAuARWS7uoM/MzMysHbYBBw58v4lkTb5xtpBxGzYHfWZmZmbtcB1wnKSF6fcXAr8qaddhJ0taQLL37meyXNxBn5mZmVk7fJxkx7NPS3okyUSNzcCXJR0jaScASTtJehrwFZLlW96S5eKeyGFmZmadtW3bNjZtytS72QVnk2y99gLgWuAe4C5gKfA1YJukTcBiHozhNgDfBlZMuriDPjMzM7MWiIifkXTn/jrwa8BjSBZjXjvv1HV5rt/5oE/Sw0imK++d/rGsIpJeCfyPiHh602WZVa7P9ZF0LPCxiFjedFnMzAZFxAUks3JL1bkxfZLWSjp+7vuIuDEi9qizgZR0rKRcUXbG6z9b0sWS7pS0dsjxtZJ+Kume9OtL846/UdIt6c9/UNJuA8dWpNe+V9L3Bv+WVr++12dJDxuox3NfIelNA2V7YN7x0wZ+fl9Jn5H0E0k/kuStIM2s8yQdOPms6XUu6OuJnwAfBP5gzDm/mgYHe0TEL829KOmXgTOA40j69x8BvGPg5z4BfItkPMBbgPMl7V9u8c22M7I+DwS5e0TEHsAvAA+w/RPuhsFzIuK8gWPnkCxXsAR4OfA+SY+t7DcxM6vHOkn/LukPJB1a1kU7FfRJ+ijwMOCf0yf+N6eZq0inLSPpEkl/lv6x7pH0z5IWS/q4pLskXS5pxcA1Hy3py5Jul3StpJcNHHuBpGsk3S1pvaTfl7Q78Hlg6UDmYWk6k+YMST+QtEnSpyTtm15nroynS9og6ea5TMYwEfGfEfFR4IYcf6bTgHMj4uqIuAP4U+CVaTkOB54InBkRP03Tx98Bfn3E33uxpIvSv9t/AofOO362pJvS41dIekb6+oFpJnHxwLm/KOk2SbtIOkzSv6WZnx9L+sccv2fnuT4P9Qrg0ohYm+HvtztJ3X1rRNyTbjp+EXDqiPMfKunDku6QdA1w9Lzjc7/v3enf6cXp67ulf89fGDj3ACXZ9v0l7Sfps5I2p+d9TekMOzOznM4GDgL+Avi+pG9Leoekxxe5aKc+mCLiVOBGHsxyvXvEqSeRfPAvIwlUvgl8iGSl6zXAmfDzRuPLwD8ABwAnA3+vBzMF5wK/GRF7AkcCX42InwDPZ/vswwbg9cCLgGeRzLK5gyQLMejZwCOBXwLOULGu1Y+nQdSX5lWCxwL/NfD9fwFL0gDsscAN6Qreg8dHZUbOAe4jqXi/kX4Nupxky5h9Sf6Gn5b0kIi4BbgEeNnAuacAn4yIrSSB6JeAfYDlwP/K9ivPFtfnoV4BnDfvtQMkbZT0Q0l/m/6eAIcDP4uI6wbOHVefzyT5+x0K/DLJA9KgHwDPAPYmyY5/TNJBEXE/8EmSOjznZOBfI+I24E0kg6r3J8k4/jE7boxuZpZZRLwpIg4BfhE4CxDwJ8C30ofT90h6uiRNc91OBX1T+FBE/CAi7iTJYvwgIv41IrYBnwaOSs97IbA2Ij4UEdsi4kqSbqWXpMe3AkdI2isi7kiPj/KbwFsiYl3aSLwdeMlcxib1joj4SUR8h6TRPjnn7/dykq7bhwMXA1+UtCg9tgfJmj1z5v6955Bjc8f3nH8DSTuTZFHelpb5u8xrjCPiYxGxKf3b/Q2wG/Co9PB5pI1keq2TgY+mx7amZV8aEfelGRobbdbrMwBppngJMLhz+vdIHiwOAp5D8gH4nvRY5vqcehnwzoi4PSJuAv5u8GBEfDoiNkTEAxHxj8D3gSelh88D/vtABu9Utq/PBwEPj4itEfG1gX0zzcxyi4irIuLtEfELJO3rm4GNJA/mlwIbJK2S9HxJu0y63qwGfRsH/v3TId/vkf774cCT026ZzZI2kwRUcwMof51krZwfpd2RTx1zz4cDnxm4zhrgZySN2JybBv79I5IMytQi4htp9+y9EfEukoUbn5EevodkevecuX/fPeTY3PG72dH+JLO755f55yS9SdKatJt2M0mGZL/08D+RBBiPAJ4L3BkR/5keezPJU8t/Srpa0vwMom1vpuvzgNNI9pi8Z+6FiLglIq5JA7EfktSduSB2mvpMWr5x9fkVkq4a+J2PJK3PEXEZydjEZ0l6NHAYSVcywF8B1wNfknSDpDOm+q3NzDKIiOsj4q8j4mnAwcBvkwzROg34F+A2Sf8g6aWjrtHFoK/MJ+ibgH+LiEUDX3tExGsAIuLyiDiRpKvs/wGfGlOGm4Dnz7vWQyJi/cA5Bw/8+2EkCyqWIUiCKICrgcHu3scDGyNiU3rsEZL2nHf86iHXvI1kD8D5ZQZ+npX5Q5LsyT4RsYgkyyKAiLiP5O/1crbPisw15P8zIpaSZJT+XtJh0/7SM8L1mWS8HfBSduzanW+wrl8HLFCyav2cUfUZkk3MR9XnhwP/F3gtsDitz98duBc8mL0+FTg/reNExN1pV8wjgF8Ffk/ScRN+DzOz3CLi5oh4fyQTOZeQDI35CsmWbCPHyXcx6NtIMiO1DJ8FDpd0qpIJBrtIOlrSYyTtKunlkvZOx6HdRZLpmCvDYkl7D1zr/cA708aDdID3ifPu91ZJC9MxVq9ixP+YdBD9Q4Bdkm/1EKX77ilZ4uKYtHwPkfQHJNmIb6Q//hHg1ZKOkLQPyRiADwOkY5+uAs5Mf/bFwOMYshZQJEuGXAi8PS3zEWw/BmpPkqDwNpKG923smHX5CMkkkhOAjw38fi+VNLc22h0kDXlf16TrdX0e8GKSjPXF83722LTOS9LBJIOa/wkgHY94IXCWpN0lHQOcyMADxjyfAv5I0j5p/XvdwLHdSerhbel9X0WS6Rv00bScp5DU7bkyvlDJ5CTx4N+1r/XZzGoWEZvT4Va/TtJL92ujzu1i0Pcu4E/SLpjfL3KhdELDL5EMlN8A3AL8JcnYNEie6NdKugv4LdIxahHxPZKlT25Iy7GUZKbNRSRdPHcD/wE8ed4t/42kG+grwF9HxJcY7pkk3XafI8lG/JRk4gMkwdb7SIKl9cDzSDIym9KyfQF4N0nj+aP068yBa58ErEx//i+Al6SD0Yd5LUnX4S0kgeOHBo59kWR82XXpPe5j+64zIuIbJMtvXBnbz8Y8GrhM0j0kf7M3pF13fdT3+jznNOAjQ8bCPZFk4spPgH8nyb69fuD4bwMPBW5Nf4fXRMSoTN87SOrqD9P7D2afrwH+Jr3XRpKlY74x+MMRsQ64kiQ4/NrAoUcC/0rS3fxN4O8j4pIRZTAzq0w69Ov/jToujzeunpIlNX4I7JIOvu8NSV8F/iEiPtB0WawcPa/PHySZ6fwnTZfFzBKPfvSj49xzz63k2k9/+tOviIiVlVw8A0lnTj5rxx+LiLcPO9D5bdisvSQdTZKpmd8taNY5abD7azw4W9rMrGpvIxlbPDieeZy5894+7KCDPquEpPNI1nl7Q2y/LqBZ50j6U+CNwLt6PBTBzJrxX8BnyriQg74apOPZplpAsesiYv7CtzYjelqf3wq8telymFkvXRURZ5VxoS5O5DAzMzOzKTnoMzMzM+uBWrt399tvv1ixYkWdt7QeueKKK34cEfvXdT/XZ6tS3fXZzGZfrUHfihUrWL16dZ23tB6R9KPJZ5XH9dmqVHd9NrPZ5+5dMzMzsx7w7F0zMzOzdjobuLysi03M9En6oKRbJX134LV9JX1Z0vfT/+5TVoHMzMzMDCLi9yLiE2VdL0v37odJ9ncddAbwlYh4JMm+m2eUVSAzMzOzvpJ0kqQXFfj5JZL+btixiUFfRFwK3D7v5ROB89J/z+28YGZmZmbFHABcIGm1pNMl7ZXlhyQ9SdI5wA3AKcPOyTumb0lE3AwQETdLOmBMIU4HTgfYiwX88W6HDj3v8D12HXmzZQcsHHls38P2HXlsn8NGr3aw96HLRh5bOGIZjl0OGv46wE5LRh/T4oOHvn6/Rv/OmzdvznVs06ZNI49t3Lhx5DGADRs2jD0OsG7duonnrF+/vpTrZC1TnbLWZxhfp6E99Rry1e1R9Rq6V7f7Xq/NrD0i4u8kfR/4G+B9wHslXQ58E/gesAn4KbAvcCBwNHAMsALYCvwfmtp7NyJWAasADtppt6j6fmZVcn02M7OqRcTnJX0BeCHwGuB44OnsuAXmXDt0M/BnwKqIGPn0mTfo2yjpoDTLdxBwa87rmJmZmeW2bdu2idn+LoqIAP4Z+GdJewLPBJ5A0v37UODHwI3ApRFxTZZr5g36LgJOA/4i/e8/5byOmZmZmY0REXcD/5J+5ZZlyZZPkPQjP0rSOkmvJgn2npv2OT83/d7MzMzMWmpipi8iTh5x6LiSy2JmZmZmFfE2bGZmZmY90Pmg7/br5y8hWJ2tN6+t7V5mZmZmZepE0Lf+1nubLoJllGUtMzMzM6tfJ4I+a4esC8+amZlZ+zjoMzMzM+uBmQ767rj+tqaLYGZmZtYKMx30mfXNnT/wmEozMxvOQZ9ZxziDbWZmeTjoM+uJe9eubboIZmbWoNYEfdfds6XW+7kbzKyY2HRT00UwM7MptCboMzNrK68/aWazwEGfmZmZWQ846DOzmbZhw4ami2Bm1goO+swqUmScap17Svdd3TvNOAg1s6Y46Buij7Mc3RD129ab1zZdhLE2bdrUdBHMzDrPQd8MccNoXbZ58+ami1CI96Y2s7Zz0FeSBzaureU+XW8YwY2jmZlZExz0WSZZAjUva2FN2LhxY9NFMDPrBAd9PeGG0czMrN8c9JnNGO82Y2Zmw8x80OfN6W0WuV7Xx8MWzGxWFAr6JL1R0tWSvivpE5IeUlbB6pAnI9L2pS2sO9bfem/t9+zjckTjeFKRmfXJgrw/KGkZ8HrgiIj4qaRPAScBHy6pbJndfv3t7HvYvnXf1qwXHti4lp2WrGi6GLl4/Umz2bdlyxa/1zMq2r27AHiopAXAQmBm/up9yohMerM4G2JmZtZ9uYO+iFgP/DVwI3AzcGdEfKmsglm3eNxT+bwVm5mZlSl30CdpH+BE4BBgKbC7pFOGnHe6pNWSVt8bD+Qvqc2EabKGbUzXuz7bMM6Gm1kXFOnePR74YUTcFhFbgQuBp80/KSJWRcTKiFi5UPlv18Sgd7P5yqrPVfMkpXI4g21ms6RIq3Uj8BRJCyUJOA5YU06xyjVueQuvaWZd1cVlW8reRrDIouPOzplZ3xQZ03cZcD5wJfCd9FqrSiqXlajqhtHZkG7pyySlNg4PMDNrUqH+qYg4MyIeHRFHRsSpEXF/WQVrg1GN4yx1g7lhtLbbtGlT00UwM5sJ7R2UNKVZmumYpwusCw2ju9Os7/yQZWZNmpmgrw0e2Li26SLYjCnyMOPxqqOVOWzBDzNm1hUO+nDjOIrH87VfFydzmJlZM3oT9LlxtDZqaimiWZ/M4W5UM7Md9Sboy2vWG0ebbdNmsWdpktIozmCbWV856Ev1sYt3XDak7IbR456qMwtZ7DyTl4osRWRm1kczFfRNGvTe9caxTw2ju+esCmU/fPhhxsy6pFNBX5vGP/WhG2wcd3/Vp+hyRKOy2B66MJzrtpnNqk4FfVXrYxfvMM5elOe6e7Y0XYTecZbYzGy4VgV9dTSQk7p4uxj45VmYuc6G0UFk9aqu102tQVnmouOewGFmfdeqoK/NutgVlmc8nxvG+tU1bGFY4NfVoQuj6vaoh5kqHjymvaYzkGbWtJkL+srYji1rVqSKxnHUZI08kzhGKdIwOuBrRpZ6nWWiUpOZ7DrqdhGu22Y262Yu6Msib+NYZrbvfu1aynVGdX9VMWs3T6Port32KzuL3VTdrjPLZ2bWRb0M+rJqy/i+NmT56syCuBssv7Zn++qWNeBzls/M+qC3QV/eNfvmZ0Xmd/EWGfA+bXBXVybEDWI7ZB26kKduT6rXVRlV58uq21k4g21mfdG5oC/LoPcyxvXNmfWsSJUBnxvG9qq7Xlc5bs9jVM2sKZKeJ+laSddLOmPI8d+TdI2kb0v6iqSHN1HOOZ0L+spUVravSk1mQtwgdleXd58pY5yqAz4zq5qknYFzgOcDRwAnSzpi3mnfAlZGxOOA84F311vK7fU66GvKsIHu02RC6ujWLdogOss3nSqWbZk28Kt7WaIysn95u3WdwTazEjwJuD4iboiILcAngRMHT4iIiyNi7gP+P4DlNZdxO70P+to88L0tS1lYe7R96ELdDzTDtHGMqicnmXXWfpJWD3ydPnBsGXDTwPfr0tdGeTXw+SoKmdWCJm9epduvv519D9u3lnttvXktuxy0opZ7Octn07jj+tvY57D9c/1slfW6qixfGwM+M6vW1q1bq2x3fhwRK0cc05DXYuiJ0inASuBZZRUsj95n+qCc2Y5FDGsAp2kU2xTwFeWMSPnKGt9X9VZsVa85WUbd9gONmQ1YBxw88P1yYIdGTNLxwFuAEyLi/prKNpSDvlQX9uQtcx9S6648Xbyj6vf8ej3Nw4wWHzz5pIrkyfKZmZXscuCRkg6RtCtwEnDR4AmSjgL+D0nAd2sDZdxOoaBP0iJJ50v6nqQ1kp5aVsH6rMtZPje8NqiqLPa0XLfNrGwRsQ14LfBFYA3wqYi4WtJZkk5IT/srYA/g05KuknTRiMvVouiYvrOBL0TES9Iod2EJZZpo/a33suyAybeadlxfkfFPWeXdospZPiuqjvqd1zT120MAzKwtIuJzwOfmvfa2gX8fX3uhxsid6ZO0F/BM4FyAiNgSEZ5uOqUiA9qd5euG6+7Zkum8KpZt6bKyZuw2PU7VzKwtinTvPgK4DfiQpG9J+oCk3eefJOn0uanO98YDBW5Xjy4vapuXG8Xs2lKf8y7dMqx+lzVeNW8Wuwv8QGNms6BI0LcAeCLwvog4CvgJsMMWJBGxKiJWRsTKher2vJE6JnMMy/wN6/oqe3/dovrSKM5SfW5C1vo9yrRdu36gMTN7UJFWax2wLiIuS78/nyQILCRrV1hWZc50zKPJGY5ZtKlR9Fit2VLF0IVh2v7A4XptZm2RO+iLiFuAmyQ9Kn3pOOCaUkpluTTVuLS90bXZ4PptZlZM0f6p1wEfl/Rt4AnAnxcvUju0ZWzfNF27o7jRmk1lbsk2i9qUxTYza4NCS7ZExFUk24pYTzmgtCK8FJGZWX08En2MSdm+oluxlbH/aBHOhFidstb3WRrPZ2bWJg76WqRIENjErF03uDbrXMfNbJZ0NuibZiFbj32yWdX07PS6eKkWM7PiOhv09UEZkzisW5ralaOONSjNzKxZDvosF3d7mZmZdYuDvp5y95eNUnSCkpmZtZODPjOzIZzNNrNZ46CvRlVtSO9tnszMzGwSB30zrEuZCgeuNq1R9dtDF8zMhnPQZ2at4dnpZmbVcdA3w5YvX950ETJbunRp00WwFliyZEnmc0fV72XLlpVVHDOzmeKgz8zMzKwHHPTVaLfYUsl1nSUzMzOzSRz0mdl2Fq5Y0XQRWqFLwyPMzLJY0HQBrBnLli0rNMtx+fLlnZodbGZms2nLli2etZ+RM31mZmZmPdCLoG/fw/Ztugi5LF68eIfXppndaN2z7ICFjdx370M949XMbNZ1NuhrqnHsmirHJXnMU/PyPNDsc9j+FZSkWtNOVvKyLWZmO+ps0FeHuhvHRYsW5f5Zz+A1K58fbMxsljjoK6DoLMciQV4ZnA2xOg2r70WHMDgoMzPLzkGfFeJG18zMrBsKB32Sdpb0LUmfLaNAbdGWcU9lTOZwYDabujpBKS+P6zMzK6aMTN8bgDUlXKcSfWoYmxrX56DShml6+EJZXL/NbFYUCvokLQd+BfhAOcWBw/fYtaxL5TYqy5dnWYvYdNNU52dtKMvK9rUpG+LJKDanrHF9bajfrtdm1hZFM33vBd4MPFBCWWyEYV28o+RpYMpoGJ0N6YZhDzR1rtGXdTLHKM5mm5nllzvok/RC4NaIuGLCeadLWi1p9b3R/tiwLWP5rJ2qrM/TrD3ZxmELu8WWUq/nhcjNzMpVJNN3DHCCpLXAJ4HnSPrY/JMiYlVErIyIlQvlycLzFRn3NKpRHJUNqbobrA/ZkDz1uQ1DFmZd27t4zczaIHcUFhF/FBHLI2IFcBLw1Yg4pbSSjZE1IzJtNqSOLF+WbEjRLjCzYdqcxXYXr5lZ9Zx6S7W5QRzH2T5rs2EPMNNkt2dpQoeZWdNKCfoi4pKIeGEZ12qrNmxIPyobMm3gN44bx/Yrc7/d+fV63C4zOy3Z/ti0M9Mn6cKEJTOzLnOmj3xZvqJbsA0qmg2ZVtXZuCLX9/IW5WtjFruJbB8UD/ycyTazLpvZoK/O2Y27HLSitnuVle1zV1g/FAn4yqrX0wR4dYzta6J++2HGzNpgZoO+rLI0imV37WZd2mLabF+Z3bxWryomJ02q21UMWSh72ZY509TtLNk4P9iYWR/1PujLo8yu3TlVZUPGqTLb526w7qmiXo8zqs5Pm80epsrAz3XbzLpqJoO+rNmQNo51mq+shtHdvP3SRJavKUUy2a7jZtYnnQv6ptm1oAxNb1FVprpnPDojUq4uP8yMqttNZ/vAgZ+Z9Ufngr6y5G0Y53eBzR/sPn9ZizJU2TBC9TMerV2GPchMqtdtlyeTXYQfaMysi3ob9GVRZZavqgHvc8rs5oV6Az9PPMmvjVm+LvCDjZn1QS+Dvrwzduse6J5FmduzVRH4OSPSLnXU62kfaJoct2pm1iczF/SVsT5f1gxfFV1g0459GqeKJVycESlflnGqWep127N8VY9ZhWKB37R1e9pg0hlsM2vagqYLULe8sxrbmOWbs3jxYjZt2jTVzyxdupQNGzbs8Pry5ctZt27d2J9dtmwZ69evn+p+1rysWb42jecbVbeXLFnCxo0bS7+f67ZZ92zZsmViu2WJVmX6Dt9j16aLYBVw91r1ql6ipYoJSnVzPTSzvmtV0Ne0Lmb5quAdDcrjB5n6uRvVzGy4TgV9RdfoK3PMU5u6wMaZdvkWa5+i41T9MPOgsh9onD00sy7pVNA3SRmTONqsjoHwbeFsTXZtn8CRRZ/qtplZU2Yq6CtilralympcYOWMiLWds9hmZtNx0DdBH7vAbHZM+zDTlWELRXjMqpn1VW+CvlnoAmsjN47F1L2X9JxZf5jx8AAzsx31Jugbp49du1m4S7Z5k8ap+mGmeX6fmFlXOOgr0SysZVYlN4718sPMaO7iNbM+mpmgr86Zu1WPe8ozk7HMPXjNrBrudjazJuUO+iQdLOliSWskXS3pDWUWrA1mfdwTlNMIOSNiVfIDjZlZOYpk+rYBb4qIxwBPAX5H0hHlFKtcfR/3VGRpC3fJzp4+PMyAs2pmZvPlDvoi4uaIuDL9993AGqCylE9Vsxw97sm6qosPM2Uvwlz1A42z2GY2S0oZ0ydpBXAUcFkZ1zMDZ2qKyPMw04c1+szM+qxw0CdpD+AC4Hcj4q4hx0+XtFrS6nvjgaK3s47rendxnfV51rcVnCVdr9dm1g+Fgj5Ju5AEfB+PiAuHnRMRqyJiZUSsXKjuTBbuy7inLNwN9qCu1mczM7Mis3cFnAusiYj3lFckq5u7UW0crz9pZjYbiqQqjgFOBZ4j6ar06wUllWsq7gYzs2EmPdC4W9bM+mRB3h+MiK8DKrEslRg3w9GD3a1JRWak552562EL01u2bBnr169vuhhmZoV5UJKZmZlZDzjoM5sxXnvSzMyGcdDXE0UWsTUzM7Puc9BnmXjZFmsrP9CYmWXjoK8kdS1rUfY2Vk3wjEkzM7P6OegboqszHBcvXpz7Z71WnzVt3ANNkbpdFz/MmFnbOegzq8jhe+zadBEy6/NSRA7WzKwvcq/TZ2bV8YLjs2vp0qVs2LCh6WKYzYytW7f6PZWRM31mNtM8dMHMLOGgz8zMzKwHWhP01T3+yQvYmhWjxQc3XYTaeDkiM5sFrQn6zMzMzKw6DvrMOmafw/bP9XPcl4yTAAAGEklEQVRdXYrIzMzK4aDPzMzMrAcc9JnNEI9VNTOzUWY66MvbDWZmZmY2a2Y66DMzMzOzhIM+y8zbVZmZmXVXJ4K+ZQcsHHnM21W1i9czMzMza6dOBH1t0edN6c3MzKzbHPSZmZmZ9YCDPjMzM7MeKBT0SXqepGslXS/pjLIKZWZmZmblyh30SdoZOAd4PnAEcLKkI8oqmJmZmZmVp0im70nA9RFxQ0RsAT4JnFhOsczMzMysTEWCvmXATQPfr0tf246k0yWtlrT63nigwO3Mmuf6bGZmXbWgwM9qyGuxwwsRq4BVACtXrow/X726wC1ny0PGHDvwwANzHeszaViVLJfrczau28XVUZ/NrF+KZPrWAQcPfL8c2FCsOGZmZmZWhSJB3+XAIyUdImlX4CTgonKKZWZmZmZlyh30RcQ24LXAF4E1wKci4uqyCmZmZmbWZpOWrpO0m6R/TI9fJmlF/aV8UJExfUTE54DPlVQWMzMzs04YWLruuSRD3i6XdFFEXDNw2quBOyLiMEknAX8J/Lf6S5vwjhxmZmZm08uydN2JwHnpv88HjlODs7QKZfrMzMzMmhQR3HfffU3cetjSdU8edU5EbJN0J7AY+HEtJZyn1qDviiuuuEfStXXeM7UfDf2BG7x33+4L8Kg6b9ZgfYb+/f/t43u41vps1lVbtmzZfOONN95c0eUfJWlwba5V6dJdkG3pukzL29Wl7kzftRGxsuZ7Iml1E/dt8t59u+/cvWu+ZSP1Gfr3/7fpetXU71z3Pc26KCL2aejWWZaumztnnaQFwN7A7fUUb0ce02dmZmY2vSxL110EnJb++yXAVyOiN5k+MzMzs85Lx+jNLV23M/DBiLha0lnA6oi4CDgX+Kik60kyfCc1V+L6g75Vk0+Zqfs2ee++3beJe/fpd+3rfZu8d5O/s5llMGzpuoh428C/7wNeWne5RlGDWUYzMzMzq4nH9JmZmZn1QC1B36RtSiq878GSLpa0RtLVkt5Q173T++8s6VuSPlvjPRdJOl/S99Lf+6k13vuN6d/5u5I+IekhFd3ng5JulfTdgdf2lfRlSd9P/1vpbK4m6nQf63N630bqdF31Ob1X43XazGZf5UHfwDYlzweOAE6WdETV901tA94UEY8BngL8To33BngDyb7EdTob+EJEPBp4fF33l7QMeD2wMiKOJBnUWtWA1Q8Dz5v32hnAVyLikcBX0u8r0WCd7mN9hgbqdM31GRqu02bWD3Vk+rJsU1KJiLg5Iq5M/303SWOxrI57S1oO/ArwgTrul95zL+CZJLOFiIgtEbG5rvuTTAx6aLoW0UJ2XK+oFBFxKTuuczS41c15wIuquHeqkTrdt/qc3rfJOl1LfYZW1Gkz64E6gr5h25TU0lANkrQCOAq4rKZbvhd4M/BATfcDeARwG/ChtBvuA5J2r+PGEbEe+GvgRuBm4M6I+FId904tiYib07LcDBxQ4b0ar9M9qc/QUJ1uQX2Geuu0mfVAHUFf41uQSNoDuAD43Yi4q4b7vRC4NSKuqPpe8ywAngi8LyKOAn5CTV1C6XijE4FDgKXA7pJOqePeDWi0TveoPkNDdbpn9dnMeqKOoC/LNiWVkbQLSQP58Yi4sKbbHgOcIGktSdffcyR9rIb7rgPWRcRc9ud8kgazDscDP4yI2yJiK3Ah8LSa7g2wUdJBAOl/b63wXo3V6Z7VZ2iuTjddn6HeOm1mPVBH0Jdlm5JKSBLJWKA1EfGeOu4JEBF/FBHLI2IFye/71YioPEsQEbcAN0ma26j9OOCaqu+buhF4iqSF6d/9OOod9D+41c1pwD9VeK9G6nTf6nN676bqdNP1Geqt02bWA5XvyDFqm5Kq75s6BjgV+I6kq9LX/jhdQXtWvQ74eBqM3AC8qo6bRsRlks4HriSZZfotKtpRQNIngGOB/SStA84E/gL4lKRXkzTYla2A3mCd7mN9hgbqdJ31GZqv02bWD96Rw8zMzKwHvCOHmZmZWQ846DMzMzPrAQd9ZmZmZj3goM/MzMysBxz0mZmZmfWAgz4zMzOzHnDQZ2ZmZtYDDvrMzMzMeuD/A1C1hDyCJJ4rAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 17 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# https://jdhao.github.io/2017/06/11/mpl_multiplot_one_colorbar/\n",
    "# https://matplotlib.org/tutorials/toolkits/axes_grid.html\n",
    "\n",
    "automated_plot = True  # set that to False if you don't want the automated 9 plots\n",
    "if automated_plot:\n",
    "    fig = plt.figure(figsize=(10, 10))\n",
    "\n",
    "    grid = AxesGrid(\n",
    "        fig,\n",
    "        111,\n",
    "        nrows_ncols=(3, 3),\n",
    "        axes_pad=0.40,\n",
    "        cbar_mode=\"single\",\n",
    "        cbar_location=\"right\",\n",
    "        cbar_pad=0.1,\n",
    "    )\n",
    "\n",
    "    for fig_num, the_ax in enumerate(grid):\n",
    "        the_ax.axis(\"equal\")\n",
    "        try:\n",
    "            time_step = plot_times[fig_num]\n",
    "            conc = save_figs[fig_num]\n",
    "            im = the_ax.contourf(x, y, conc, 20, cmap=cm)\n",
    "            the_ax.set_title(f\"timestep {time_step} days\")\n",
    "        except IndexError:\n",
    "            fig.delaxes(the_ax)\n",
    "\n",
    "    cbar = grid.cbar_axes[0].colorbar(im)\n",
    "    cbar.set_label_text(\"Concentration (mg/L)\", rotation=270, size=20, va=\"bottom\")\n",
    "    fig.suptitle(\n",
    "        \"Evolution of the concentration through time\", y=0.9, size=25, va=\"bottom\"\n",
    "    )\n",
    "    fig.savefig(\"evolution.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Conclusions\n",
    "\n",
    "What you may have noticed, is that, even for small simple 2d transient problems like the one you have just solved, the computation times are already becoming significant...\n",
    "\n",
    "This is partly because we are dealing with big matrix which are filled with zeros. It is a complete waste of time and memory to deal with all of these 0 values. There are other ways to make our calculation way faster. We will probably dedicate a lecture to understand how we can improve this."
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "all",
   "formats": "ipynb,python//py:percent",
   "notebook_metadata_filter": "all"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  },
  "nbsphinx": {
   "execute": "never"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}