"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0,
"nbgrader": {
"grade": false,
"grade_id": "cell-46d188f356ba20fd",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"source": [
"# Week 6 assignment: Introduction\n",
"\n",
"**The one-dimensional steady-state finite-volume approximation**\n",
"\n",
"We just developed the fundamental stencil for 1-D steady-state diffusion in terms of concentrations in gridblocks as:\n",
"\n",
"## 1-D steady - state diffusion stencil (aligned in the x direction)\n",
"\n",
"For every cell (gridblock) within the domain (not on the boundary), we can write that (for the steady-state), the sum of all fluxes is zero.\n",
"\n",
"\\begin{align}\n",
"&\\left(D\\theta \\frac{c_E - c_C}{\\Delta x} + D\\theta \\frac{c_W - c_E}{ \\Delta x} \\right) (\\Delta y) (\\Delta z) = 0 \\\\\n",
"\\end{align}\n",
"\n",
"For the cells corresponding to the boundary problem, the equations are simple:\n",
"\n",
"\\begin{equation}\n",
"c_{\\mathrm{left}} = \\text{Specified value for the left side of the domain}\n",
"\\end{equation}\n",
"\\begin{equation}\n",
"c_{\\mathrm{right}} = \\text{Specified value for the right side of the domain}\n",
"\\end{equation}\n",
"\n",
"## Translation into a computational problem\n",
"\n",
"We have seen that this stencil can be applied to any grid block and this lead to a system of $n$ equations, where $n$ is the number of gridblocks. These equations can be written in a matrix form\n",
"\n",
"\\begin{equation}\n",
"Ac = b\n",
"\\end{equation}\n",
"\n",
"where A is a 2D matrix ($n \\times n$), c is the solution vector (array of size $n$) and b is the right hand side (array of size $n$).\n",
"\n",
"Apart from the boundary conditions, every line of this matrix represents the equation above:\n",
"\n",
"\\begin{equation}\n",
"\\left( -2 \\frac{D\\theta}{\\Delta x} c_C + \\frac{D\\theta}{\\Delta x}c_E + \\frac{D\\theta}{\\Delta x} c_W \\right) \\Delta y \\Delta z = 0\n",
"\\end{equation}\n",
"\n",
"In a computational sense, this can be written as:\n",
"\n",
"\\begin{equation}\n",
"\\left( -2 \\frac{D\\theta}{\\Delta x} c_{\\left[i\\right]} + \\frac{D\\theta}{\\Delta x}c_{\\left[i-1\\right]} + \\frac{D\\theta}{\\Delta x} c_{\\left[i+1\\right]} \\right) \\Delta y \\Delta z = 0\n",
"\\end{equation}\n",
"\n",
"Therefore, for line $i$, representing the $i^{th}$ gridblock, the value for the $i-1$ columnn of matrix $A$ should be\n",
"\\begin{equation}\n",
"A\\left[i \\right]\\left[i-1 \\right] = + \\frac{D\\theta}{\\Delta x}\\Delta y \\Delta z\n",
"\\end{equation}\n",
"\n",
"The value representing the contribution from the cell on the right ($i+1$ column) is\n",
"\\begin{equation}\n",
"A\\left[i \\right]\\left[i+1 \\right] = + \\frac{D\\theta}{\\Delta x}\\Delta y \\Delta z\n",
"\\end{equation}\n",
"\n",
"And the term on the diagonal is\n",
"\\begin{equation}\n",
"A\\left[i \\right]\\left[i\\right] = -2 \\frac{D\\theta}{\\Delta x}\\Delta y \\Delta z\n",
"\\end{equation}\n",
"\n",
"The boundary conditions will be written as:\n",
"\n",
"\\begin{equation}\n",
"c{\\left[0\\right]} = c_{\\text{left}}\n",
"\\end{equation}\n",
"\\begin{equation}\n",
"c{\\left[n-1\\right]} = c_{\\text{right}}\n",
"\\end{equation}\n",
"So that means that the matrix A will have:\n",
"\\begin{equation}\n",
"A\\left[0 \\right]\\left[0 \\right] = 1\n",
"\\end{equation}\n",
"\\begin{equation}\n",
"A\\left[n-1 \\right]\\left[n-1 \\right] = 1\n",
"\\end{equation}\n",
"\n",
"And the boundary conditions have to be in the vector $b$ (the right-hand-side):\n",
"\\begin{equation}\n",
"b{\\left[0 \\right]} = c_{\\text{left}}\n",
"\\end{equation}\n",
"\\begin{equation}\n",
"b{\\left[n-1 \\right]} = c_{\\text{right}}\n",
"\\end{equation}\n",
"\n",
"## Building these matrixes\n",
"\n",
"In the cell below, we have defined a function which builds the matrix A and the vector b. A lot of comments are given. We strongly advise you to make the link between the code and the previous equation. We don't need you to be able to write this alone, or by memory, but to be able to understand these lines and adapt them."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-1a2af915a90d524e",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from numpy.testing import assert_almost_equal\n",
"import pdb"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"nbgrader": {
"grade": false,
"grade_id": "cell-07bfc86663956631",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"outputs": [],
"source": [
"def build_1d_matrix(c_left, c_right, n, D, width, poro):\n",
" \"\"\"\n",
" Constructs a coefficient matrix A and an array b related to the problem Ac = b.\n",
" \n",
" Parameters\n",
" ----------\n",
" \n",
" c_left: float\n",
" left boundary condition for concentration [mg/L]\n",
" c_right: float\n",
" right boundary conditions for concentration [mg/L]\n",
" n: int \n",
" number of cells/gridblocks [-]\n",
" D: float \n",
" value of the diffusion coefficient (supposed the same everywhere) [m2/s]\n",
" width float: \n",
" total physical width of the domain [m]\n",
" poro: float\n",
" porosity value, supposed the same everywhere [-]\n",
" \n",
" Returns\n",
" -------\n",
" \n",
" A, b: tuple\n",
" A: 2-d np.float array \n",
" b: 1-d np.float vector\n",
" \n",
" These solve the \n",
" discretized 1D diffusion problem Ax = b\n",
" \"\"\"\n",
" # The name of this function is build_1d_matrix. We have to know its name to use it!\n",
" # In this function, we will use the values of every parameters which are given to the function\n",
" # ( c_left, c_right, n ,..): we have to use the same names and not overwrite these!\n",
" A = np.zeros((n, n)) # This is the initialization of the square matrix \"A\".\n",
" # it is a n*n matrix which we initially fill with 0\n",
" b = np.zeros(n)\n",
" # b is the 1D array representing the RHS\n",
" dx = width / (n - 1)\n",
" # for a 10m width with 1 meter between cels, we need 11 cells\n",
" # so, dx, which is the length between the cells, would be 1m (10/(11-1))\n",
" dy = dz = dx # here it is assumed that all grid cells are cubes\n",
" coef = poro * D / dx * dy * dz\n",
" # this coefficient is the one which appears in front of every term in the conservation equation\n",
"\n",
" # Now we are going to put the coefficients in the matrix A and in vector b.\n",
" # So we are going to loop over every equation, each of them being described by one line of the matrix A\n",
" for i in range(n):\n",
" if i == 0:\n",
" # The first line represents the left boundary condition\n",
" # the \"index\" i representing the left-value is 0\n",
" # The equation for the left gridblock is: 1 * c[0] = c_left.\n",
" # So, in the matrix form, the coefficicent associated with c[0] in line 0 is 1.\n",
" A[i][i] = 1\n",
" # while the RHS associated to that is bc_left\n",
" b[i] = c_left\n",
" elif i == n - 1:\n",
" # Same for the right side. the equation there is c[n-1] = c_right\n",
" A[i][i] = 1\n",
" b[i] = c_right\n",
" else:\n",
" # For every other node, we have to write the long equation.\n",
" East = coef\n",
" West = coef\n",
" # In the future, these might not be precisely the same, so we define one for East and one for West.\n",
" A[i][i] = East + West\n",
" A[i][i + 1] = -East\n",
" A[i][i - 1] = -West\n",
" # As we indexed the unknowns from left to right, the unknow corresponding\n",
" # to the concentration in the gridblock right to the i-th gridblock is the i+1 unknown\n",
"\n",
" # And, as there is nothing in the RHS in the equation above, the right hand side should be zero.\n",
" # But, as it was initialized as 0, there is no need to be b[i] = 0 (it is alreay done !)\n",
"\n",
" # We could have written things like this below, this would be precisely the same:\n",
" # A[i][i] = 2*coef\n",
" # A[i][i+1] = - coef\n",
" # A[i][i-1] = -coef\n",
" return A, b\n",
" # This last line is very important. It says that the function has to give back the matrix A and the array b.\n",
" # Basically, when you call that function, specifying all the arguments, the function will give back these\n",
" # two arrays, so that you can work with them."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"lines_to_next_cell": 2,
"nbgrader": {
"grade": false,
"grade_id": "cell-e5cda80cdf50fd9f",
"locked": true,
"schema_version": 1,
"solution": false
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd4VHXaxvHvk1BCLxKU3ov0EhEpwUJXwQV07XVVVCywvq5l11XXXdsu2LCvfa0oiAKCWAhVCNKrAZEqREEQpOd5/5ghOxsDDCEnk3J/rmsu5pw5Z+bmMOTJac/P3B0RERGAuFgHEBGR/ENFQUREMqkoiIhIJhUFERHJpKIgIiKZVBRERCSTioKIiGRSURARkUwqCiIikqlYrAMcqypVqnjdunVjHUNEpECZO3fuj+6eeLTlClxRqFu3LqmpqbGOISJSoJjZ99Esp8NHIiKSSUVBREQyqSiIiEgmFQUREcmkoiAiIpkCKwpm9rKZbTGzxYd53czsSTNLM7OFZtYuqCwiIhKdIPcUXgV6H+H1PkCj8OM64NkAs4iISBQCKwrungJsPcIi/YHXPWQWUNHMqgWVZ1X6Tv41aQV79h8M6iNERAq8WJ5TqAGsi5heH573G2Z2nZmlmllqenp6jj7ss6WbeeqLNM5+ciqpa45Uq0REiq5YFgXLZp5nt6C7v+DuSe6elJh41Lu0szW4WwNeu7oDe/ZncP7zM7lv7BJ27T2Qo/cSESmsYlkU1gO1IqZrAhuD/MBujROZODSZK06ry2sz19BzRAopK3O25yEiUhjFsiiMBS4PX4XUEdju7puC/tCyJYtxX7/mvH/9aZQsHsflL8/m9vcX8POv+4L+aBGRfC/IS1LfBmYCTcxsvZldY2aDzWxweJHxwGogDXgRuDGoLNlJqluZ8bd05aYzGjB63ga6D09hwqLAa5KISL5m7tkexs+3kpKSPLe7pC7ZuJ07Ri1kycYd9GlxEvf3b07Vcgm5+hkiIrFkZnPdPeloy+mOZqB59QqMuakzd/RuwufLt9BjeAqj5q6noBVMEZHjpaIQVjw+jhtPb8iEW7vS+MSy3P7+Aq54ZQ7rt/0a62giInlGRSGLBollefe603igf3NS12yl54gUXpuxhowM7TWISOGnopCNuDjj8tPqMmloMkl1K/PXsUu44PmZpG3ZGetoIiKBUlE4gpqVSvPaVafwr/Nb8+2WnfR9Yiojv0xj/8GMWEcTEQmEisJRmBkD29dk8rBu9Gh2Io9NXEH/p6ezeMP2WEcTEcl1KgpRSixXkpGXtOO5S9uTvnMv/UdO55FPl6vBnogUKioKx6h3i5OYPLQbA9vV4NmvVtH3ianMUYM9ESkkVBRyoELp4jw6qDVvXnMq+w5mcP5zM7n3o8XsVIM9ESngVBSOQ5dGVZg0NJmrO9fjjVnf02tECl+t2BLrWCIiOaaicJxKlyjGvec2Y9TgTpQqEc+Vr8xh2Hvz2bZLDfZEpOBRUcgl7etUYtwtXbjlzIaMnb+RHiOmMH7RJrXKEJECRUUhF5UsFs+wnk0YO6QL1SqU4sb/fMPgN+eyZceeWEcTEYmKikIAmlUvz+gbO3FXn6Z8tSKds4ZP4b0567TXICL5nopCQIrFx3F9twZ8elsyJ1crzx0fLOSyf89m3VY12BOR/EtFIWD1qpThnWs78uB5LZi/7md6jkjh5WnfcVAN9kQkH1JRyANxccalHeswaWgyHetX5oFPljLouRl8u/mXWEcTEfkfKgp5qHrFUrx85Sk8/vs2rPlxF2c/OY0nP/+WfQfUYE9E8gcVhTxmZpzXtgafDetGrxYnMfyzlfR7ehoL1/8c62giIioKsVKlbEmeuqgtL16exLZf93HeyOk8NH6ZGuyJSEypKMRYj2Yn8tmwbvz+lFo8n7Ka3o+nMGv1T7GOJSJFlIpCPlA+oTgPDWjFW384lQyHC1+YxT2jF/HLnv2xjiYiRYyKQj7SqWEVJt6WzB+61OPt2WvpOSKFL5erwZ6I5B0VhXymVIl4/nxOMz64oRPlEopx1atzuO2deWxVgz0RyQMqCvlU29qV+OTmrtx6ViPGLdpE9+FTGLtgo1pliEigVBTysRLF4hjaozEf39yFWpVKccvb87j29bn8sF0N9kQkGCoKBUDTk8rz4Y2duafvyUxLS6fH8Cm8PXut9hpEJNepKBQQ8XHGtcn1mXhbMi1qVOCuDxdx8Ytf8/1Pu2IdTUQKERWFAqbOCWV469pTeWhASxZv2E6vx1N4aepqNdgTkVyholAAmRkXdajNZ8O60aVhFR4ct4wBz0xnxQ9qsCcixyfQomBmvc1shZmlmdmd2bxe28y+NLN5ZrbQzPoGmaewOalCAi9ensSTF7Vl/bbdnPPUVEZ8tlIN9kQkxwIrCmYWD4wE+gDNgIvMrFmWxf4MvOfubYELgWeCylNYmRn9Wlfns2HdOLtlNZ74/FvOeWoq89epwZ6IHLsg9xQ6AGnuvtrd9wHvAP2zLONA+fDzCsDGAPMUapXLlODxC9vy8pVJ/LLnAAOemc6Dnyxl9z412BOR6AVZFGoA6yKm14fnRboPuNTM1gPjgZsDzFMknNn0RCYNTeaiDrV5adp39Ho8hRlpP8Y6logUEEEWBctmXtZLZC4CXnX3mkBf4A0z+00mM7vOzFLNLDU9PT2AqIVLuYTi/P13LXnnuo7EGVz80tfc9eFCdqjBnogcRZBFYT1QK2K6Jr89PHQN8B6Au88EEoAqWd/I3V9w9yR3T0pMTAwobuHTsf4JfHpbMtd3q8+7c9bRY/gUJi/dHOtYIpKPBVkU5gCNzKyemZUgdCJ5bJZl1gJnAZjZyYSKgnYFclFC8Xju6nMyY27qTKXSJfjD66nc/PY8ftq5N9bRRCQfCqwouPsBYAgwEVhG6CqjJWb2gJn1Cy/2R+BaM1sAvA1c6erdEIhWNSsydkgX/tijMRMX/0D34VMYM2+DWmWIyP+wgvZDISkpyVNTU2Mdo0D7dvMv3PHBQuat/Zkzm1blwfNaUL1iqVjHEpEAmdlcd0862nK6o7kIanRiOUYN7sS95zRj5qqf6DkihTdnfU+GWmWIFHkqCkVUfJxxdZd6TBqaTJtaFfnzmMVc+OIsvvtRDfZEijIVhSKuVuXSvHFNBx4d2Iplm3bQ+/EUnpuyigMH1SpDpCgqdrQFzKwq0BmoDuwGFgOp7q6fGoWEmXHBKbXo1iSRv4xZzMMTljNu4SYeGdiKZtXLH/0NRKTQOOyegpmdYWYTgXGE+hdVI9TD6M/AIjO738z0E6MQObF8As9f1p5nLmnHpu276ff0NP41aQV7D6hVhkhRcaQ9hb7Ate6+NusLZlYMOAfoAXwQUDaJATOjb8tqnFb/BP42bilPfZHGhMU/8MjAVrSvUynW8UQkYDm6JNXMBrp7TIqBLknNW1+t2MI9oxezcfturuxUl//r1YTSJY561FFE8pmgL0kdkcP1pIA5vUlVJg5N5rKOdXhl+hp6jkhh2rdqsCdSWOW0KGTX7E4KqbIli/FA/xa8d/1pFI+P49J/f80doxawfbca7IkUNjktCrrLqQjqUK8yE27tyg2nN+CDbzbQY/gUJi75IdaxRCQXHfacgpktIvsf/gY0dveSQQY7HJ1TyB8Wrd/OHR8sZNmmHZzdqhr3nducxHIx+UqISBSiPadwpKJQ50gruvv3Ocx2XFQU8o/9BzN4fsoqnvw8jdIl47n3nGb8rm0NzHR0USS/yY0TzS8AA4BS7v591keuJZUCq3h8HEPObMT4W7tQv0oZhr23gKtencOGn3fHOpqI5NCRisIVwDbgPjP7xsyeNbP+ZlY2j7JJAdGwajneH9yJ+85txuzvttJz+BTemLlGDfZECqCo7lMID5F5KqE7m88i1O5ikrs/Gmy839Lho/xt3dZfuXv0IqZ++yOn1K3EwwNb0SBRv0eIxFqu3qfg7hnuPtPd73X3zoRGUdtwvCGl8KlVuTSvX92Bxwa1YsUPv9Dniak881WaGuyJFBDRNMR7it9ehbQd0K/rki0z4/ykUIO9e8cs4dFPVzB+UajBXvPqFWIdT0SOIJo9hZJAG+Db8KMVUBm4xsweDzCbFHBVyyXw3GXtefaSdvywfS/9np7OYxOXs2e/GuyJ5FfRNLFpCJwZHnMZM3sWmESoGd6iALNJIdGnZTVOa3ACD45bxsgvV/Hp4h94dFAr2tepHOtoIpJFNHsKNYAyEdNlgOrufhDYG0gqKXQqli7BP89vzetXd2DP/gwGPTeT+8YuYdfeA7GOJiIRoikKjwLzzewVM3sVmAf808zKAJODDCeFT3LjRCYNTeaK0+ry2sxQg72UlemxjiUiYdFekloN6ECoxcVsd98YdLDD0SWphUfqmq3c8cFCVqfvYlD7mvz57JOpWLpErGOJFEq53To7kdAVSHFARzMbcDzhRACS6lZm/C1duemMBoyet4Huw1OYsGhTrGOJFGlH3VMws5cJXXG0BDh0sbm7+9UBZ8uW9hQKpyUbt3PHqIUs2biDPi1O4v7+zalaLiHWsUQKjeNuiBfxRkvdvVmuJTtOKgqF14GDGbw49TtGTF5JQrE4/nJOMwa1r6kGeyK5IDcPH800s3xTFKTwKhYfxw2nN2DCrV1pclI5/m/UQi5/eTbrtv4a62giRUY0ReE1QoVhhZktNLNFZrYw6GBSdDVILMu7153G3/o355vvt9Hr8RRenf6dGuyJ5IFoDh+lAcMI3aiW2cBG4ylIXtjw827u/nARU1am075OJR4Z2JKGVcvFOpZIgZObh4/WuvtYd/9O4ylIXqtRsRSvXnUKwy9ozar0nfR9Yhojv0xjvxrsiQQimjYXy83sLeBjIu5gdvcPA0slEsHMGNCuJl0bJXLfx0t4bOIKPlm4iccGtaJFDTXYE8lN0ewplCJUDHoC54Yf5wQZSiQ7ieVKMvLidjx/WXt+2rmX/iOn8/AENdgTyU1H3VNw96ty+uZm1ht4AogHXnL3h7NZ5gLgPkI3xy1w94tz+nlSNPRqfhId653AP8Yv47kpq5i05AceHtiKDvXUYE/keB12T8HM/mxmh/1fZmZnmtlh9xjMLB4YSWi0tmbARVkvbTWzRsBdQGd3bw7cdoz5pYiqULo4jwxqxZvXnMq+gxlc8PxM/jJmMTvVYE/kuBxpT2ER8LGZ7QG+AdKBBKARofEVJgP/OML6HYA0d18NYGbvAP2BpRHLXAuMdPdtAO6+JYd/DymiujSqwqShyfxz4kpemfEdny/bzN8HtOSMJlVjHU2kQDrsnoK7fxQeenMwoRYX8cAO4E2gg7sPdfcjtbesAayLmF4fnhepMdDYzKab2azw4abfMLPrzCzVzFLT09VRU/5X6RLFuPfcZowa3IkyJYtx1StzGPbufLbt2hfraCIFTjTnFA6NuHassutNkPWmiGKE9jxOB2oCU82shbv/nCXDC8ALELpPIQdZpAhoX6cSn9zShZFfpPHMV6uYsjKd+/s35+yW1dQqQyRK0XZJzYn1QK2I6ZpA1pbb64GP3H2/u38HrCBUJERypGSxeIb1bMLHN3ehesVSDHlrHte/MZfNO/bEOppIgRBkUZgDNDKzemZWArgQGJtlmTHAGQBmVoXQ4aTVAWaSIuLkauUZfWMn7u7blCkr0+k+fArvzllLNOOHiBRlgRWF8JjOQ4CJwDLgPXdfYmYPmFm/8GITgZ/MbCnwJfB/7v5TUJmkaCkWH8d1yQ2YeFsyzaqV508fLOKSl75m7U9qsCdyONH0PkokdJVQXSLOQWg8BSlIMjKct+es5aHxyzmY4dzeqwlXdqpLfJzONUjREG3vo2jaXHwETCV0CapuHZUCKS7OuOTUOpzZtCr3jF7M3z5ZyscLNvLooFY0PlEN9kQOiWZPYb67t8mjPEelPQU5Xu7O2AUbuW/sEnbuPcCQMxpxw+kNKFEsyFNsIrGVm11SPzGzvrmQSSRfMDP6t6nB5GHd6N2iGiMmr6Tf09NYsO7no68sUshFUxRuJVQY9pjZL+HHjqCDiQTthLIleeqitrx0eRI//7qf3z0znYfGL2P3Ph0llaLrqEXB3cu5e5y7J4Sfl3P38nkRTiQvdG92IpOGJfP7U2rzfMpq+jyRwsxVughOiqaoDqKaWT8z+2f4obbZUuiUTyjOQwNa8ta1p+LARS/O4u7Ri9ixZ3+so4nkqaMWBTN7mNAhpKXhx63heSKFTqcGVfj01mT+0KUe78xeS8/hKXyxfHOsY4nkmWiuPloItHH3jPB0PDDP3VvlQb7f0NVHklfmr/uZP41ayIrNv9C/TXXuPacZJ5QtGetYIjmSm1cfAVSMeK7xD6VIaFOrIh/f3IXbujdi/KJN9BiRwkfzN6hVhhRq0RSFh4B5Zvaqmb0GzOXI4yiIFBolisVxW/fGfHJzV2pVLs2t78znD6+lsmn77lhHEwnEUQ8fAZhZNeAUQu2wv3b3H4IOdjg6fCSxcjDDeWX6d/xz0gqKx8VxV9+TufCUWsSpVYYUAMd9+MjMmob/bAdUI9Tmeh1QPTxPpEiJjzP+0LU+E29LpkWNCtw9ehEXvzSLNT/uinU0kVxz2D0FM3vB3a8zsy+zednd/cxgo2VPewqSH7g7785Zx9/HLWN/RgZ/7NGEq7vUU4M9ybei3VOI5uqjBHffc7R5eUVFQfKTH7bv4c9jFjF52RZa16zAI4Na0fQk3dsp+U9uXn00I8p5IkXOSRUSePHyJJ66qC3rtu3m3KemMeKzlew7kBHraCI5ctjW2WZ2ElADKGVmbfnvmMvlgdJ5kE2kQDAzzm1dnc4Nq/DAx0t44vNvmbB4E48Oak2bWhWP/gYi+ciRzilcAVwJJAGRx2t+AV519w8DT5cNHT6S/O6L5Zu5Z/RiNu/Yw9Wd6/HHnk0oVSI+1rGkiMvNcwoD3f2DXEt2nFQUpCD4Zc9+Hp6wnP98vZbalUvz8ICWdGpYJdaxpAjLtaIQfrOzgeZAwqF57v7AcSXMIRUFKUhmrf6JOz9YyJqffuXCU2pxV9+TqVCqeKxjSRGUayeazew54PfAzYTOK5wP1DnuhCJFQMf6JzDh1mSuT67Pe6nr6DliCp8tVYM9yb+iufqok7tfDmxz9/uB04BawcYSKTxKlYjnrr4nM+amzlQqXYJrX09lyFvf8OPOvbGOJvIb0RSFQ/cj/Gpm1YH9QL3gIokUTq1qVmTskC78sUdjJi3ZTI/hUxgzTw32JH+Jpih8bGYVgceAb4A1wNtBhhIprEoUi+Pmsxox7pYu1K1Shtvenc/Vr85h489qsCf5wxFPNJtZHNDR3WeEp0sCCe6+PY/y/YZONEthcTDDeW3GGh6buIL4OONPfZpySYfaarAngciVE83hgXX+FTG9N5YFQaQwiY8zru5Sj0lDk2lTqyJ/GbOYC1+cxer0nbGOJkVYNIePJpnZQDPTry8iAahVuTRvXNOBRwe2YtmmHfR5YirPTVnFgYNqlSF5L5qb134BygAHCJ10NkJdUmPS9UuHj6Qw27xjD38Zs5hJSzfTokZ5Hh3YmmbV1WBPjl+u3afg7uXcPc7dS7h7+fC0vqUiATixfALPX9aekRe344fte+j39DT+NWkFew8cjHU0KSKiuXnt82jmiUjuMDPOblWNz4Z2o1+b6jz1RRpnPzmNud9vi3U0KQKONPJagplVBqqYWSUzqxx+1AWq51VAkaKqUpkSDL+gDa9edQq79x1k0HMzuP/jJezaeyDW0aQQO9KewvXAXKBp+M9Dj4+AkdG8uZn1NrMVZpZmZnceYblBZuZmdtTjXSJFzelNqjJxaDKXdazDK9PX0OvxFKZ+mx7rWFJIRXOi+WZ3f+qY39gsHlgJ9CA0vvMc4CJ3X5pluXLAOKAEMMTdj3gWWSeapSib/d1W7vxgIat/3MUFSTW5p28zKpRWgz05utw80fyUmXUys4vN7PJDjygydADS3H21u+8D3gH6Z7Pc34BH+W87DRE5jA71KjP+1q7ceHoDPvhmA91HTOHTxT/EOpYUItGcaH4D+CfQBTgl/IjmME8NYF3E9PrwvMj3bgvUcvdPog0sUtQlFI/njt5N+eimziSWLcngN+dy43/mkv6LGuzJ8TvscJwRkoBmfuxdu7K72S3zPcItNEYQGt3tyG9kdh1wHUDt2rWPMYZI4dSiRgU+GtKZF1JW88Tn3zI9bQr3ntOMAe1qoHtNJaeiuaN5MXBSDt57Pf/bYrsmsDFiuhzQAvjKzNYAHYGx2Z1sdvcX3D3J3ZMSExNzEEWkcCoeH8dNZzRk/C1daVi1LH98fwFXvDKH9dt+jXU0KaCiKQpVgKVmNtHMxh56RLHeHKCRmdUzsxLAhUDmeu6+3d2ruHtdd68LzAL6He1Es4j8VsOqZXn/+tO4v19zUtdspdeIFF6fuYaMDLXllmMTzeGj+3Lyxu5+wMyGABOBeOBld19iZg8Aqe4eTWERkSjFxRlXdKrLmU2rcvfoRdz70RI+XrCRhwe2okFi2VjHkwIi2jGa6wCN3H2ymZUG4t39l8DTZUOXpIocnbszau56Hhy3jN37D3Jb90Zc27U+xeOjOTgghVFujtF8LTAKeD48qwYw5vjiiUiQzIzzk2rx2bBkzmpalUc/XcF5I6ezeIM638uRRfNrw01AZ2AHgLt/C1QNMpSI5I6q5RJ49tL2PHtJOzbv2Ev/kdN59NPl7NmvBnuSvWiKwt7wzWcAmFkxIi4tFZH8r0/Lanw+rBsD2tbgma9W0ffJqaSu2RrrWJIPRVMUppjZ3UApM+sBvA98HGwsEcltFUoX57HzW/P61R3Yuz+D85+fyV8/WsxONdiTCNEUhTuBdGARoSZ544E/BxlKRIKT3DiRSUOTueK0urw+63t6jUhhyko12JOQaBrilQH2uPvB8HQ8UNLdY3J3jK4+Esk9c7/fyh2jFrIqfRcD2tXg3nOaUbF0iVjHkgDk2tVHwOdAqYjpUsDknAYTkfyjfZ3KjLulK0POaMjY+RvpPnwK4xdtinUsiaFoikKCu+88NBF+Xjq4SCKSlxKKx3N7ryZ8NKQzJ1VI4Mb/fMPgN+ayZYcaFxdF0RSFXWbW7tCEmbUHdgcXSURioXn1Coy5sTN/6t2UL1ZsofvwKbyXuo5j74UpBVk05xROITQWwqFmdtWA37v73ICzZUvnFESCtzp9J3d+sIjZa7bStVEV/vG7ltSqrAMEBVm05xSibXNRHGhCqB32cnfff/wRc0ZFQSRvZGQ4/5m9lofHLyPD4Y7eTbj8tLrEx6ktd0GUmyeaITSwTiugLXBRlCOviUgBFhdnXNaxDpOGdePU+pW5/+OlnP/cDNK2xKTtmeSRIEdeE5FCoEbFUrxy5SmM+H1rVv+4i75PTOPpL75l/8GMWEeTAAQ58pqIFBJmxu/a1qRro0T+OnYJ/5y0kk8WbuKxQa1pWbNCrONJLgpy5DURKWSqlC3JyIvb8fxl7dm6ax/9R07joQnL1GCvEIlmT+HQyGuzgcyRwd29X2CpRCRf69X8JDrWP4F/jFvG81NWM2nJZh4e0JJT658Q62hynKK5JLVbdvPdfUogiY5CVx+J5C/T037kzg8Xsm7rbi7tWJs/9W5KuYTisY4lWeTa1UfhH/7LgXLhx7JYFQQRyX86N6zCxNuSuaZLPf7z9Vp6jkjhy+VbYh1Lciiaq48uAGYD5wMXAF+b2aCgg4lIwVG6RDH+ck4zPrihE2VLFuOqV+cw9N35bN217+grS74SzeGjBUAPd98Snk4EJrt76zzI9xs6fCSSv+09cJCRX67imS/TqFCqOPf3b87ZLathppveYik3b16LO1QQwn6Kcj0RKYJKFotnWI/GfHxzF2pUKsWQt+Zx3Rtz2awGewVCND/cPzWziWZ2pZldCYwDJgQbS0QKupOrlefDGzpxd9+mpKxMp/vwKbwze60a7OVz0fY+GkDojmYDUtx9dNDBDkeHj0QKnjU/7uJPHyzk6++20qnBCTw8oBW1T1CDvbx03A3xzKwhcKK7T88yPxnY4O6rciXpMVJRECmYMjKct2av5eEJyzmQkcHtPZtwVed6arCXR3LjnMLjQHadr34NvyYiErW4OOPSjnX4bFgynRpU4cFxyxj47AxWblaDvfzkSEWhrrsvzDrT3VOBuoElEpFCrVqFUvz7iiSeuLAN3/+0i7OfnMoTk79l3wE12MsPjlQUEo7wWqkjvCYickRmRv82NZg8rBu9W1RjxOSV9Ht6GgvW/RzraEXekYrCHDO7NutMM7sGiMmoayJSuJxQtiRPXdSWFy9PYtuv+/jdM9P5x/hl7N6nBnuxcqQTzScCo4F9/LcIJAElgN+5+w95kjALnWgWKZx27NnPQ+OX8/bstdQ5oTQPD2jFaQ3UYC+35NpwnGZ2BtAiPLnE3b/IhXw5pqIgUrjNWPUjd324iO9/+pWLT63NnX2aUl4N9o5bro7RnJ+oKIgUfrv3HWTE5JW8NHU1Vcsl8PffteCsk0+MdawCLbfHaM5piN5mtsLM0szszmxeH2ZmS81soZl9bmZ1gswjIgVDqRLx3N33ZEbf2JmKpYtzzWup3PL2PH7auffoK8txCawomFk8MBLoAzQDLjKzZlkWmwckuXsrYBTwaFB5RKTgaV2rImOHdGFo98ZMWLyJHiNS+Gj+BrXKCFCQewodgDR3X+3u+4B3gP6RC7j7l+7+a3hyFlAzwDwiUgCVKBbHrd0bMe6WrtSuXJpb35nPH15LZdP23bGOVigFWRRqAOsipteH5x3ONRym0Z6ZXWdmqWaWmp6enosRRaSgaHxiOT64oRN/Pvtkpq/6kZ7DU3jr67VkZGivITcFWRSya2iS7b+emV1K6HLXx7J73d1fcPckd09KTEzMxYgiUpDExxl/6Fqfibcl06JGBe4evYiLX5rFmh93xTpaoRFkUVgP1IqYrglszLqQmXUH7gH6ubvOIonIUdU5oQxvXXsqDw1oyZINO+j1eAovpKziwEG1yjheQRaFOUAjM6tnZiWAC4GxkQuYWVvgeUIFQYO6ikjUzIyLOtTms2Hd6NqoCv8Yv5yBz85g+Q87Yh2tQAusKLj7AWAIMBFYBrzn7kvFGSD7AAAOJ0lEQVTM7AEz6xde7DGgLPC+mc03s7GHeTsRkWydVCGBFy9P4qmL2rJ+227OeXIawz9byd4DapWRE7p5TUQKja279vG3T5Yyet4GGp9YlkcGtqJt7UqxjpUv5Iub10RE8lLlMiUY8fs2vHxlEr/sOcCAZ2fwt0+W8uu+A7GOVmCoKIhIoXNm0xOZNDSZS06tzb+nfUfvx6cyI+3HWMcqEFQURKRQKpdQnAfPa8m713UkPs64+KWvufODhWzfvT/W0fI1FQURKdROrX8CE27tyvXd6vNe6jp6DJ/CpCUx6fxfIKgoiEihl1A8nrv6nMyYmzpTuUwJrntjLkPe+oYf1WDvN1QURKTIaFWzIh/f3IU/9mjMpCWb6T58CqPnrVeDvQgqCiJSpBSPj+Pmsxox7pYu1KtShqHvLuDqV+ew8Wc12AMVBREpohqdWI5Rgzvx13ObMWv1VnqOSOGNWd8X+QZ7KgoiUmTFxxlXda7HpKHJtKlVkb+MWcyFL8xidfrOWEeLGRUFESnyalUuzRvXdODRQa1Y/sMO+jwxlWe/KpoN9lQUREQINdi7IKkWk4d14/QmiTzy6XLOe2Y6SzcWrQZ7KgoiIhGqlk/g+cuSePaSdvywfS/9np7GPyeuYM/+otFgT0VBRCQbfVpWY/KwZPq3qcHTX6Zx9pNTmfv91ljHCpyKgojIYVQsXYJ/XdCa167uwJ79GQx6bib3jV3Crr2Ft8GeioKIyFF0a5zIxKHJXNaxDq/OWEPPESmkrCyc48WrKIiIRKFsyWI80L8F7w8+jZLF47j85dnc/v4Ctv9auBrsqSiIiByDU+pWZvwtXbnx9AaMnreB7iOm8OniTbGOlWtUFEREjlFC8Xju6N2Uj27qTGLZkgx+8xtueHMuW37ZE+tox01FQUQkh1rUqMBHQzrzf72a8PnyLfQYnsKouQW7wZ6KgojIcSgeH8dNZzRk/C1daVS1LLe/v4DLX57Nuq2/xjpajqgoiIjkgoZVy/Le9adxf7/mzP1+G70eT+HV6d8VuAZ7KgoiIrkkLs64olNdJg1NJqluZe77eCkXPD+TtC0Fp8GeioKISC6rWak0r111Cv86vzXfbtlJ3yemMvLLNPYXgAZ7KgoiIgEwMwa2r8nkYd3o3qwqj01cQf+np7N4w/ZYRzsiFQURkQAllivJM5e057lL25G+cy/9R07nkU+X59sGeyoKIiJ5oHeLakwe2o0BbWvw7Fer6PvkVOasyX8N9lQURETySIXSxXns/Na8cU0H9h3I4PznZnLvR4vZmY8a7KkoiIjksa6NEpl4WzJXda7LG7O+p9eIFL5asSXWsQAVBRGRmChTshh/Pbc5owafRkLxOK58ZQ7D3pvPtl37YppLRUFEJIba16nMuFu6MuSMhoydv5EeI6YwbuGmmLXKCLQomFlvM1thZmlmdmc2r5c0s3fDr39tZnWDzCMikh8lFI/n9l5N+GhIZ06qkMBNb33D4DfnsmVH3jfYC6womFk8MBLoAzQDLjKzZlkWuwbY5u4NgRHAI0HlERHJ75pXr8CYGzvzp95N+XJFOt2HT+G91HV5utcQ5J5CByDN3Ve7+z7gHaB/lmX6A6+Fn48CzjIzCzCTiEi+Viw+jhtOb8Cnt3al6UnluWPUQi77d9412AuyKNQA1kVMrw/Py3YZdz8AbAdOCDCTiEiBUD+xLO9c15G/ndeCeWu30XNECh8v2Bj45wZZFLL7jT/rPlA0y2Bm15lZqpmlpqcXznFRRUSyioszLutYh0nDutG5YRXqVSkT/GcG+N7rgVoR0zWBrGUucxkzKwZUAH5zi5+7v+DuSe6elJiYGFBcEZH8qUbFUrx0RRItalQI/LOCLApzgEZmVs/MSgAXAmOzLDMWuCL8fBDwhRfkIYtERAq4YkG9sbsfMLMhwEQgHnjZ3ZeY2QNAqruPBf4NvGFmaYT2EC4MKo+IiBxdYEUBwN3HA+OzzLs34vke4PwgM4iISPR0R7OIiGRSURARkUwqCiIikklFQUREMqkoiIhIJitotwWYWTrwfQ5XrwL8mItxcotyHRvlOnb5NZtyHZvjyVXH3Y9692+BKwrHw8xS3T0p1jmyUq5jo1zHLr9mU65jkxe5dPhIREQyqSiIiEimolYUXoh1gMNQrmOjXMcuv2ZTrmMTeK4idU5BRESOrKjtKYiIyBEUmqJgZr3NbIWZpZnZndm8XtLM3g2//rWZ1Y147a7w/BVm1iuPcw0zs6VmttDMPjezOhGvHTSz+eFH1rbjQee60szSIz7/DxGvXWFm34YfV2RdN+BcIyIyrTSznyNeC3J7vWxmW8xs8WFeNzN7Mpx7oZm1i3gtkO0VRaZLwlkWmtkMM2sd8doaM1sU3lapuZXpGLKdbmbbI/697o147YjfgYBz/V9EpsXh71Tl8GuBbDMzq2VmX5rZMjNbYma3ZrNM3n2/3L3APwi15l4F1AdKAAuAZlmWuRF4Lvz8QuDd8PNm4eVLAvXC7xOfh7nOAEqHn99wKFd4emcMt9eVwNPZrFsZWB3+s1L4eaW8ypVl+ZsJtWQPdHuF3zsZaAcsPszrfYEJhEYT7Ah8nQfb62iZOh36LKDPoUzh6TVAlRhur9OBT473O5DbubIsey6hMV4C3WZANaBd+Hk5YGU2/x/z7PtVWPYUOgBp7r7a3fcB7wD9syzTH3gt/HwUcJaZWXj+O+6+192/A9LC75cnudz9S3c/NCL3LEIj1AUtmu11OL2Az9x9q7tvAz4Desco10XA27n02Ufk7ilkMypghP7A6x4yC6hoZtUIcHsdLZO7zwh/JuTdd+vQZx9tex3O8Xw3cztXnny/3H2Tu38Tfv4LsIzfjmefZ9+vwlIUagDrIqbX89uNmrmMux8AtgMnRLlukLkiXUPot4FDEiw0NvUsMzsvlzIdS66B4V3VUWZ2aGjVfLG9wofZ6gFfRMwOantF43DZg9xexyLrd8uBSWY218yui0EegNPMbIGZTTCz5uF5+WJ7mVlpQj9cP4iYHfg2s9Bh7bbA11leyrPvV6CD7OQhy2Ze1suqDrdMNOvmVNTvbWaXAklAt4jZtd19o5nVB74ws0XuviqPcn0MvO3ue81sMKG9rDOjXDfIXIdcCIxy94MR84LaXtGIxfcrKmZ2BqGi0CVidufwtqoKfGZmy8O/ReeVbwi1XdhpZn2BMUAj8sH2CjsXmO7ukXsVgW4zMytLqAjd5u47sr6czSqBfL8Ky57CeqBWxHRNYOPhljGzYkAFQruR0awbZC7MrDtwD9DP3fcemu/uG8N/rga+IvQbRJ7kcvefIrK8CLSPdt0gc0W4kCy79gFur2gcLnuQ2+uozKwV8BLQ391/OjQ/YlttAUaTe4dMo+LuO9x9Z/j5eKC4mVUhxtsrwpG+X7m+zcysOKGC8B93/zCbRfLu+5XbJ01i8SC0x7Oa0OGEQyenmmdZ5ib+90Tze+HnzfnfE82ryb0TzdHkakvoxFqjLPMrASXDz6sA35JLJ9yizFUt4vnvgFn+3xNb34XzVQo/r5xXucLLNSF00s/yYntFfEZdDn/i9Gz+90Tg7KC3VxSZahM6R9Ypy/wyQLmI5zOA3rm5raLIdtKhfz9CP1zXhrddVN+BoHKFXz/0C2OZvNhm4b/368DjR1gmz75fufoliOWD0Nn5lYR+wN4TnvcAod++ARKA98P/SWYD9SPWvSe83gqgTx7nmgxsBuaHH2PD8zsBi8L/KRYB1+RxroeAJeHP/xJoGrHu1eHtmAZclZe5wtP3AQ9nWS/o7fU2sAnYT+i3s2uAwcDg8OsGjAznXgQkBb29osj0ErAt4ruVGp5fP7ydFoT/je/JzW0VZbYhEd+vWUQUruy+A3mVK7zMlYQuPolcL7BtRuiwngMLI/6t+sbq+6U7mkVEJFNhOacgIiK5QEVBREQyqSiIiEgmFQUREcmkoiAiIplUFEREJJOKgshxMLPBZnb5MSzf1sxeOsbPeMfMGh17OpFjp/sURPKQmb0PPOjuC45hnW7Ape5+bXDJREK0pyBFgpmdEu74mmBmZcKDmbTIZrlzLTQI0zwzm2xmJ4bnP3loIBgz62VmKWYWZ2b3mdnt4fm32H8HTHonm/cuB7Q6VBDC675mZpPCA7gMMLNHwwO5fBruhwMwFege7tklEih9yaRIcPc5FhqN7UGgFPCmu2c3+tY0oKO7u4VGm7sD+CNwJzDHzKYCTwJ93T0jNCRHpjuBeh7qLFsxm/dOArJ+ZgNCAy01A2YCA939DjMbTajfzZjw56QBrYG5OdoAIlFSUZCi5AFgDrAHuOUwy9QE3g0PYFKCUIMx3P1XM7sWSAGGevYtuRcC/zGzMYRaQWdVDUjPMm+Cu+83s0WERh37NDx/EaHGbYdsAaqjoiAB0+EjKUoqA2UJDXmYAGBmfz80Jm94macIDUPaErj+0HJhLYGfCP1wzs7ZhJqWtQfmZnO4Z3eW9wPYC+DuGcB+/+9Jvgz+95e2hPD6IoFSUZCi5AXgL8B/gEcA3P0ed2/j7m3Cy1QANoSfZw6CHh7p7Y+EWp33MbNTI9/YzOKAWu7+JaFDThUJFaBIy4CGOczemFB3TpFA6fCRFAnhy0YPuPtbZhYPzDCzM939iyyL3ge8b2YbCLV0rhcey/vfwO0eGnnrGuBVMzslYr144E0zq0CozfEId/858o3dfbmZVTCzch4aizfa7CcCu9190zH+tUWOmS5JFclDZjYU+MXdo75XIbzODnf/d3DJREJ0+Egkbz1L+DzCMfiZ0BjZIoHTnoKIiGTSnoKIiGRSURARkUwqCiIikklFQUREMqkoiIhIpv8HsYEpE54m+SsAAAAASUVORK5CYII=\n",
"text/plain": [
"