{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Step 1: 1D Linear Convection\n", "\n", "The first step in this journey begins with a simple equation called the 1-D Linear Convection equation. It is the most basic model that we can use to learn something useful about CFD.\n", "Here it is:\n", "\n", "$$\\frac{\\partial u}{\\partial t} + c \\frac{\\partial u}{\\partial x} = 0$$\n", "\n", "This equation can be uderstood as a representation of a *wave* that propagates whatever initial shape it had at a constant speed c without changing it. Assuming initial condutions to be $u(x,0) = u_0(x)$ we can find an exact solution to this equation of the type $u(x,t) = u_0(x-ct)$. Which represents the behavior of a wave being propagated with no change in shape. \n", "\n", "But we must discretize the above equation before we can work with it on our computer. To do so we perform a discretization in both space and time using the Forward Difference scheme for the tiem derivative and the backward difference scheme for the spatial derivative. \n", "\n", "This means that from now on we will consider space to have been discretized into **x** points that we index with i= 0 to **N** and step time in discrete intervals of size $\\Delta t$.\n", "\n", "Performing a taylor series approximation with a forward difference in time and eliminating the higher order terms we can obtain a discretization for the time derivative:\n", "\n", "$$\\frac{\\partial u}{\\partial t} \\approx \\frac{u_i^{n+1} - u^n_i}{\\Delta t}$$\n", "\n", "And similarly for the spatial with a backwards difference:\n", "\n", "$$\\frac{\\partial u}{\\partial x} \\approx \\frac{u_{i}^{n} - u^n_{i-1}}{\\Delta t}$$\n", "\n", "Combining and substituting into our equation:\n", "\n", "$$\\frac{u_i^{n+1} - u^n_i}{\\Delta t} + c \\frac{u_{i}^{n} - u^n_{i-1}}{\\Delta t} = 0$$\n", "\n", "Our initial condition information will be enough to solve for everything in the above equation except $u_i^{n+1}$ so we can rearrange to solve as such:\n", "\n", "$$u_i^{n+1} = u^n_i - c \\frac{\\Delta t}{\\Delta x}(u_i^n - u^n_{i-1})$$\n", "\n", "Now that we have fully defined the equation, discretized it and found an expression that we can easily solve for we are going to implement it in python. \n", "\n", "### Libraries" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Adding inline command to make plots appear under comments\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import time, sys\n", "%matplotlib inline " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Key Variables\n", "\n", "Before beginning to calculate we must first define the grid of points that will characterize our spatial domain. We begin by giving it a grid_length parameter that defines the length of our spatial domain. Next, grid_points will define the number of points in the grid and dx will be the resulting distance between any pair of adjacent grid points. nt will represent the number of timesteps we want to calculate and dt represents the amount of time each timestep will cover. Finally c will represent the wave speed of our system. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "grid_length = 10\n", "grid_points = 41\n", "dx = grid_length / (grid_points - 1) \n", "nt = 400\n", "dt = 0.025\n", "c = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initial Conditions\n", "\n", "Next we set up our initial conditions for the wave. $u = 1$ everywhere except:\n", "\n", "$$u = 2 \\Rightarrow 0.5 <= x <= 1$$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1. 1. 2. 2. 2. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.\n", " 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n" ] } ], "source": [ "u = np.ones(grid_points)\n", "u[int(.5/ dx):int(1 / dx + 1)] = 2\n", "print(u)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting initial conditions\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAHGpJREFUeJzt3X2UZHV95/H3p58bBgGdkeV58ARNMMkEHDckuskoZgWDkqwkEY0YDId9UFezObs+rEo2MYkm6olZEzhzyGQ06hiDGAjHB1wTnWRdPA5EeUwIaoDhaUZQhKGqu7rqu3/ce2vKsru6pvveur+e/rzO6cNM1e17v1XMqU/9ft97708RgZmZGcBY3QWYmVk6HApmZtblUDAzsy6HgpmZdTkUzMysy6FgZmZdDgVLgqS3Sbqq7joOF5I+I+k1dddha49DwZYk6fWS9kiak7Sz77ltkjqSnsh/9kr6hKTnDtjfZkkhaaL/uYj4vYi4tIKXsSKSXixpt6THJe2X9CVJL6u7rsVI+i1JH+l9LCLOi4gPVXCsnZLetcLffaqkT0k6IOkeSa8suz5bPYeCDfIA8C5gx1LPR8QG4CjgbOCfgL+XdM6I6ls1SeOLPHYh8FfAh4GTgOOAdwIvHW11h50/AebJ3s9XAVdIena9JdkPiAj/+GfgD1kw7Ox7bBuwd5FtPwjsWWI/m4EAJhZ57reAj/Rt9xrgXuDbwP/s2XYMeAvwDeAR4BPAU3ue/yvgIeAxYDfw7J7ndgJXAJ8GDgAv6qtD+TH/+4D3Ywx4O3APsI8sPI5ernbgBKDRV+uZ+TaT+d9fC9wJfAf4HHBqz7bPBj4PPAo8DLwNOJfsg7YFPAF8Pd/2i8Clq6l3kdd9WX6c+fxYf3MI/4aOzH/vmT2P/QXw7rr/ffvn+388UrCyXQOcJenIEvb1fOBZwDnAOyX9SP74G4BfAH6W7IP2O2TfQgufAU4Hng7cDHy0b7+vBH6XbITzD33PPQs4Gbh6QF2/lv+8AHgGsIEsDAfWHhEPAP8PeHlfLVdHREvSBWQf9P8B2AT8PbALQNJRwP8BPpu/5h8CvhARnwV+D/jLiNgQEVvKqrd/JxGxney9/IP8WC/Na7te0neX+Lk+//VnAgsRcVfPLr9OFnSWEIeCle0Bsm/bx5Swr/8VEY2I+DrZB0jxgfefyL7N7o2IObJRxoVFryIidkTE4z3PbZF0dM9+r42I/xsRnYho9h3zafl/HxxQ16uA90fENyPiCeCtwCv6eiVL1f4x4CIASQJekT9WvK7fj4g7I2KB7MP+JySdCpwPPBQR74uIZv76vjLozSup3mVFxPkRccwSP+fnm20Avtf3q4+RBbMlxKFgZTuRbDriuyXs66GePz9J9sECcCrwqeLbKNl0Sxs4TtK4pHdL+oak7wH/mv/Oxp593TfgmI/k/z1+wDYnkE3FFO4BJsjmyper/ZPAT0k6HvgZoEM2Iihe1wd6XtejZAF7Itno5RsDahpkNfWW5QngKX2PPQV4vOTj2Co5FKxsvwjcHBEHKjzGfcB5fd9IZyLifrLpmAuAFwFHk82ZQ/bhWhh0a+B/zvf/8gHbPED2AV44BVggm+cfKCK+A9wA/Epe68cjoqjnPuA/9r2u2Yj4cv7cM5ba7TKHXXG9wxwrP/31iSV+PpNvdhcwIen0nl/dAty+ghqsQg4FW5KkCUkzwDgwLmlmsdNJlTlR0uXApWTz4oNM5/sqfg713+GVwO/m0ypI2pTPx0M2HTFH9o3/CLIpmKHlH9D/DXiHpEskPUXSmKTnS9qeb7YL+A1Jp0nawME5/YUhD/Mx4GLgQg5OHRWv663FGTmSjpb0S/lz1wPHS3qTpGlJR0n6yfy5h4HNA97H1dbb62H6wimy0183LPFzXr7NAbJ+029LOlLS88jC+y9WUINVyKFgg7yd7GyZtwC/mv/57T3PnyDpCbKpga8CPwZsi4gbltnvE/m+ip8XHmJdHwCuA26Q9DhwI1B8QH6YbHrkfuCO/LlDEhFXk32Tfy3Zt+yHyc7AujbfZAfZh9lu4FtAk6z5PazryBrhD+Vz+MVxPwW8B/h4PvV1G1B8qD4O/BzZabEPAf9C1jiG7GwrgEck3bzI8VZbb68/A87Ip7j++hB/978As2RnQO0C/nNEeKSQGB0cuZqZ2XrnkYKZmXU5FMzMrMuhYGZmXQ4FMzPr+oHTC1O3cePG2Lx5c91lmJmtKTfddNO3I2LTctutuVDYvHkze/bsqbsMM7M1RdI9y2/l6SMzM+vhUDAzsy6HgpmZdTkUzMysy6FgZmZdlYWCpJMl/Z2kOyTdLumNi2wjSX8s6W5Jt0g6q6p6zMxseVWekroA/GZE3JwvJXiTpM9HxB0925xHdrfI08nucnkFB+92aWZmI1ZZKETEg+RLGkbE45LuJFtBqjcULgA+nN/D/kZJx0g6Pv/dZEQEH7nxHvY/PreyHUi8bMvx/NDTvfKgmaVtJBevSdoMnAn0ryl7It+/NOLe/LHvCwVJlwGXAZxyyilVlbmkBx5r8o5rb89rOfTfj4CHH2vyngt/vOTKzMzKVXko5Cs9fRJ4U0T0L9w9lIjYDmwH2Lp168gXgHhyLlug6n9fdCYv3XLCIf/+C977RZ5stcsuy8ysdJWefSRpkiwQPhoR1yyyyf1kC5IXTsofS0oj/0CfnRxf0e/PTI7TmHcomFn6qjz7SGRL990ZEe9fYrPrgIvzs5DOBh5LrZ8AdD/QZ6dWFgqzk2M0PVIwszWgyumj5wGvBm6V9LX8sbcBpwBExJXAp4GXAHcDTwKXVFjPihUjhZkVjhRmp8a7+zAzS1mVZx/9AzCwLZufdfS6qmooS3OV00ezk+N850CrzJLMzCrhK5qH0O0prHD6aGZynOaCRwpmlj6HwhCarQ6wupFC041mM1sDHApDKBrNM5Mre7tmJt1TMLO1waEwBDeazWy9cCgModlqI8H0xMpHCs1Wh05n5NfdmZkdEofCEBrzbWYnx9FK7nHBwV7E3EKnzLLMzErnUBhCo9VecZMZsovXiv2YmaXMoTCEZquz4n4CHDyV1Vc1m1nqHApDaLbaK75GAQ42qD1SMLPUORSG0Gi1V3w6KvSEgq9VMLPEORSGUDSaV6r4XU8fmVnqHApDyEYKq+8pePrIzFLnUBhCc9VnH3n6yMzWBofCEBpuNJvZOuFQGMKqRwp5oMy1fPGamaXNoTCExvwqewoeKZjZGuFQGMJqL16b8RXNZrZGOBSWsdDuMN/urGr6aGbCjWYzWxscCsto5jexm51a+Vs1NiamJ8Z8nYKZJc+hsIzi2/1qRgrgNRXMbG1wKCyjucoFdgqzk+MeKZhZ8hwKyyg+yFdznQJkodDwKalmljiHwjKKKZ/VTh/NTI670WxmyXMoLKP4IF/t9NHMpBvNZpY+h8IyGmX1FNxoNrM1wKGwjGZJ00eznj4yszXAobCMRkmN5hmffWRma4BDYRmN+fziNZ+SambrgENhGaVNH7mnYGZrgENhGd1G8ypucwHFdQoOBTNLm0NhGc1WGwmmxlf3Vk1PjtNsdeh0oqTKzMzKV1koSNohaZ+k25Z4/mhJfyPp65Jul3RJVbWsRmM+W2BH0qr2U0w/zS34qmYzS1eVI4WdwLkDnn8dcEdEbAG2Ae+TNFVhPSvSWOWqa4VZr6lgZmtAZaEQEbuBRwdtAhyl7Cv4hnzbharqWalGa3WrrhWKU1odCmaWsjp7Ch8EfgR4ALgVeGNELDq3IukySXsk7dm/f/8oa8zWZ17lNQpw8Ipon5ZqZimrMxReDHwNOAH4CeCDkp6y2IYRsT0itkbE1k2bNo2yRpqt1a26Vuiu0+yrms0sYXWGwiXANZG5G/gW8MM11rOoxny7u8byanikYGZrQZ2hcC9wDoCk44BnAd+ssZ5FuadgZuvJRFU7lrSL7KyijZL2ApcDkwARcSXwO8BOSbcCAt4cEd+uqp6VarbaPP2o6VXvx9NHZrYWVBYKEXHRMs8/APz7qo5flkbJjWaPFMwsZb6ieRnFxWurVQTLnJfkNLOEORSWUVpPwSMFM1sDHArLmGt1Spk+ciiY2VrgUBhgod1hvt1hZmL1oTA9kd/mwo1mM0uYQ2GAZn7zutlV3jYbYGxMTE+M+ToFM0uaQ2GA4lt9GY1m8EI7ZpY+h8IAxbf6MhrNkC+04+kjM0uYQ2GA4lt9GY1myNdp9noKZpYwh8IAZU8fzXikYGaJcygMUEwfldlTcKPZzFLmUBigmD6aLm2kMOZGs5klzaEwQOkjBU8fmVniHAoDlN1onpn09JGZpc2hMEBjPr94rcyRgkPBzBLmUBig4Uazma0zDoUBuhevlXCbC/BIwczS51AYoNlqMyaYGi/nbcp6Ch06nShlf2ZmZXMoDNCYz9ZSkFTK/orbZcz5qmYzS5RDYYBGq5xV1wqzk2Pd/ZqZpcihMEBZq64VilNbHQpmliqHwgDNVru0axTg4PSRz0Ays1Q5FAZozJc9fTTe3a+ZWYocCgOU3lOY8kjBzNLmUBig2eowU+L0UXek4FAws0Q5FAZottrMTJT3Fs14+sjMEudQGKBRUaPZIwUzS5VDYYDSG83uKZhZ4hwKA5R+nUL3lFRf0WxmaXIoDFD2dQpuNJtZ6hwKS1hod2i1o9Tpo+m8ae1Gs5mlqrJQkLRD0j5Jtw3YZpukr0m6XdKXqqplJZoL5S6wAzA2JmYmx9xTMLNkVTlS2Amcu9STko4B/hR4WUQ8G/ilCms5ZMW3+ZnJct+iGa+pYGYJqywUImI38OiATV4JXBMR9+bb76uqlpXoLrBT4kgB8oV2PH1kZomqs6fwTOBYSV+UdJOki5faUNJlkvZI2rN///6RFNddirPERjN49TUzS1udoTABPAf4eeDFwDskPXOxDSNie0RsjYitmzZtGklxxbf5MnsKcHD1NTOzFE3UeOy9wCMRcQA4IGk3sAW4q8aaurojhbKnj6bG3Wg2s2TVOVK4Fni+pAlJRwA/CdxZYz3fpwiFMm+IB54+MrO0VTZSkLQL2AZslLQXuByYBIiIKyPiTkmfBW4BOsBVEbHk6aujNlfRSGFmcpxHD8yXuk8zs7JUFgoRcdEQ2/wh8IdV1bAajYrOPvJ1CmaWMl/RvITGfPkXrxX78/SRmaXKobAEN5rNbD1yKCyhe/HaVLlvkUcKZpYyh8ISGvNtxgRT4+Xf5qLZ6tDpRKn7NTMrg0NhCc1WtsCOpFL3W1whPbfgC9jMLD0OhSWUvRRnwWsqmFnKHApLaLTaTE+UHwrFXVcdCmaWIofCEspeda1QXPfgO6WaWYocCktozLdLPx0VetdpdiiYWXocCktotCoKhSmHgpmly6GwhEarU/rN8MCNZjNLm0NhCXOtNrMlL8UJ7imYWdocCkuoevrIIwUzS9FQd0mV9M7FHo+I3y63nHQ05tul3yEVDo4U3FMwsxQNe+vsAz1/ngHOJ6EFcarQaFUTCrOePjKzhA0VChHxvt6/S3ov8LlKKkpEVdcpdE9J9W0uzCxBK+0pHAGcVGYhKWm1O7TaUUlPYXoiv6LZIwUzS9CwPYVbgeK2nuPAJuCw7Sc0K1pLAWBsTF59zcySNWxP4fyePy8AD0fEQgX1JKHZyqZ2qrhOAbymgpmla9iewj1VF5KSKkcKxX49fWRmKfJ1CosovsXPVHDxWrZfjxTMLE0OhUUU3+KrGikUq6+ZmaXGobCIRtXTR1PjbjSbWZIcCovoTh+50Wxm64xDYRHNEUwfudFsZilyKCyiueDpIzNbnxwKi2jM59cpVDVSmBjz9JGZJcmhsIhRNJodCmaWIofCIprdRnM1b8/spKePzCxNDoVFNObbjAmmxqu7eK3Z6tDpxPIbm5mNkENhEcWqa5Iq2X9xS+453z7bzBJTWShI2iFpn6TbltnuuZIWJF1YVS2HqlHRWgqF7kI7nkIys8RUOVLYCZw7aANJ48B7gBsqrOOQNStada3gUDCzVFUWChGxG3h0mc3eAHwS2FdVHStRdShMT3qhHTNLU209BUknAr8IXDHEtpdJ2iNpz/79+yuvrTHfrux0VOhZktMjBTNLTJ2N5j8C3hwRy3ZbI2J7RGyNiK2bNm2qvLCi0VyVol/hUDCz1Ay78loVtgIfz8/w2Qi8RNJCRPx1jTUB0Gh1OHp2srL9u6dgZqmqLRQi4rTiz5J2AtenEAiQ3RDv3zxlurL9F/0K9xTMLDWVhYKkXcA2YKOkvcDlwCRARFxZ1XHLMKrpI48UzCw1lYVCRFx0CNv+WlV1rERzRNcpuKdgZqnxFc2LaLTaTE9UFwqePjKzVDkUFjGykYJvc2FmiXEo9Gm1O7TaUWlPYXrCF6+ZWZocCn2aFa+lADA2JmYmx9xTMLPkOBT6NLprKVQXCpCFjs8+MrPUOBT6NPOlOKscKRT79/SRmaXGodCnuVD99BFkIxGPFMwsNQ6FPsW395nJat+amQkvyWlm6XEo9GmMoNEM2VXNzZZPSTWztDgU+rjRbGbrmUOhT3N+RD0FN5rNLEEOhT6jnT5yKJhZWhwKfbqhUPn00Zinj8wsOQ6FPkXzt8o1msE9BTNLk0OhTzGlU/kpqe4pmFmCHAp9GvNtxgRT49WHwtxCh04nKj2OmdmhcCj0KVZdy9eOrkzRs5jz7bPNLCEOhT6NitdSKBRnN7mvYGYpcSj0ac63K28yg0PBzNLkUOhTTB9Vrbhi2s1mM0uJQ6FP1UtxFrpLcnqkYGYJcSj0abTazEyMYKSQn/LqUDCzlDgU+jRancpvhgfuKZhZmhwKfZrzbWYrvnANDl4x7Z6CmaXEodBnVI3mom/hkYKZpcSh0GfU1ym4p2BmKXEo9Bn5dQqePjKzhDgU+jQXRj195NtcmFk6HAo9Wu0OrXaMZKQwPeFTUs0sPQ6FHs0RrboGIImZyTGHgpklpbJQkLRD0j5Jty3x/Ksk3SLpVklflrSlqlqGVZwJNIrrFMAL7ZhZeqocKewEzh3w/LeAn42IHwN+B9heYS1Dac5n8/ujGCkUx3Gj2cxSMlHVjiNit6TNA57/cs9fbwROqqqWYTVGOH0E2YjEIwUzS0kqPYVfBz6z1JOSLpO0R9Ke/fv3V1ZEt6cwNZq3ZXZy3D0FM0tK7aEg6QVkofDmpbaJiO0RsTUitm7atKmyWro9hVFOHzkUzCwhlU0fDUPSjwNXAedFxCN11gKjD4UZh4KZJaa2kYKkU4BrgFdHxF111dGrOT/inoIbzWaWmMpGCpJ2AduAjZL2ApcDkwARcSXwTuBpwJ9KAliIiK1V1TOMUTeaZ6fcUzCztFR59tFFyzx/KXBpVcdfiW4ojOw6hTFPH5lZUmpvNKekmMpxo9nM1iuHQo+5hdFevDYz5Z6CmaXFodCjMd9mfExMjmskx5udHGduoUOnEyM5npnZchwKPRqtNjMTY+SN78oV01TFCMXMrG4OhR6jWnWt0F1ox30FM0uEQ6HHqFZdKzgUzCw1DoUejdZoVl0rFLfodrPZzFLhUOhR1/SRL2Azs1Q4FHo0W54+MrP1zaHQo9HqjHT6qLhFt0cKZpYKh0KPrNE8urdkesI9BTNLi0Ohx6gbzUX/wtNHZpYKh0IPN5rNbL1zKPSo7ToFTx+ZWSIcCj3qmz7ybS7MLA0OhVyr3WGhEyMNhemJ7O13T8HMUuFQyDVHvMAOgKTsTqkOBTNLhEMhV3xbnx7hSAFgxquvmVlCHAq55vxoF9gpzE56oR0zS4dDIdddn3nUI4UpL8lpZulwKOS6oTA12rdkdnLc1ymYWTIcCrliCmeU1ylAPn3kUDCzRDgUcs2FeqaPZqfGafo6BTNLhEMh16xppDA94UazmaXDoZCrq9GcjRQcCmaWBodCrlHDxWsAs75OwcwS4lDIudFsZuZQ6GrWeZ2CewpmlgiHQq7RajM+JibHNdLjzk6OM7fQodOJkR7XzGwxDoVcM1+fWRp9KADMLfi0VDOrn0Mh12iNdn3mQtHDcF/BzFJQ2aegpB2S9km6bYnnJemPJd0t6RZJZ1VVyzBGvepaYdahYGYJqfKr8U7g3AHPnwecnv9cBlxRYS3LGvWqa4WZKS/JaWbpmKhqxxGxW9LmAZtcAHw4IgK4UdIxko6PiAerqOdLd+3nXdffseTzD3y3wTM2baji0AMVQfTanV/trsRmZraYX3nuyVz6755R6TEqC4UhnAjc1/P3vfljPxAKki4jG01wyimnrOhgG6YnOP24pT/0Tz9uAz93xnEr2vdqbD31WF5+1kk0WgsjP7aZrS0bN0xXfow6Q2FoEbEd2A6wdevWFZ27+ZxTj+U5pz6n1LrKcOyRU7zvl7fUXYaZGVDv2Uf3Ayf3/P2k/DEzM6tJnaFwHXBxfhbS2cBjVfUTzMxsOJVNH0naBWwDNkraC1wOTAJExJXAp4GXAHcDTwKXVFWLmZkNp8qzjy5a5vkAXlfV8c3M7ND5HEgzM+tyKJiZWZdDwczMuhwKZmbWpazfu3ZI2g/cs8Jf3wh8u8Ry1gK/5vXBr3l9WM1rPjUiNi230ZoLhdWQtCcittZdxyj5Na8Pfs3rwyhes6ePzMysy6FgZmZd6y0UttddQA38mtcHv+b1ofLXvK56CmZmNth6GymYmdkADgUzM+taN6Eg6VxJ/yzpbklvqbueqkk6WdLfSbpD0u2S3lh3TaMgaVzSP0q6vu5aRiVfyvZqSf8k6U5JP1V3TVWS9Bv5v+nbJO2SNFN3TVWQtEPSPkm39Tz2VEmfl/Qv+X+PLfu46yIUJI0DfwKcB5wBXCTpjHqrqtwC8JsRcQZwNvC6dfCaAd4I3Fl3ESP2AeCzEfHDwBYO49cv6UTgvwJbI+JHgXHgFfVWVZmdwLl9j70F+EJEnA58If97qdZFKAD/Frg7Ir4ZEfPAx4ELaq6pUhHxYETcnP/5cbIPihPrrapakk4Cfh64qu5aRkXS0cDPAH8GEBHzEfHdequq3AQwK2kCOAJ4oOZ6KhERu4FH+x6+APhQ/ucPAb9Q9nHXSyicCNzX8/e9HOYfkL0kbQbOBL5SbyWV+yPgfwCdugsZodOA/cCf59NmV0k6su6iqhIR9wPvBe4FHiRbsfGGeqsaqeN6Vqh8CDiu7AOsl1BYtyRtAD4JvCkivld3PVWRdD6wLyJuqruWEZsAzgKuiIgzgQNUMKWQinwO/QKyMDwBOFLSr9ZbVT3yhcpKv6ZgvYTC/cDJPX8/KX/ssCZpkiwQPhoR19RdT8WeB7xM0r+STQ++UNJH6i1pJPYCeyOiGAVeTRYSh6sXAd+KiP0R0QKuAX665ppG6WFJxwPk/91X9gHWSyh8FThd0mmSpsgaU9fVXFOlJIlsnvnOiHh/3fVULSLeGhEnRcRmsv+/fxsRh/03yIh4CLhP0rPyh84B7qixpKrdC5wt6Yj83/g5HMaN9UVcB7wm//NrgGvLPkBlazSnJCIWJL0e+BzZ2Qo7IuL2msuq2vOAVwO3Svpa/tjbIuLTNdZk1XgD8NH8C883gUtqrqcyEfEVSVcDN5OdYfePHKa3u5C0C9gGbJS0F7gceDfwCUm/TraEwC+Xflzf5sLMzArrZfrIzMyG4FAwM7Muh4KZmXU5FMzMrMuhYGZmXQ4FMzPrciiYmVmXQ8FslSQ9V9ItkmYkHZnf6/9H667LbCV88ZpZCSS9C5gBZsnuRfT7NZdktiIOBbMS5LeY+CrQBH46Ito1l2S2Ip4+MivH04ANwFFkIwazNckjBbMSSLqO7JbdpwHHR8Tray7JbEXWxV1Szaok6WKgFREfy9cD/7KkF0bE39Zdm9mh8kjBzMy63FMwM7Muh4KZmXU5FMzMrMuhYGZmXQ4FMzPrciiYmVmXQ8HMzLr+Pxaki4oDzG32AAAAAElFTkSuQmCC\n", "text/plain": [ "