{ "cells": [ { "cell_type": "markdown", "id": "d562e52a-a449-41a3-a42e-6fd1d30da8ce", "metadata": {}, "source": [ "# 1D Convection-diffusion example\n", "\n", "An example adapted from the FiPy convection diffusion 1D example\n", "\n", "see: https://pages.nist.gov/fipy/en/latest/generated/examples.convection.exponential1D.mesh1D.html\n", "\n", "Written by Ali A. Eftekhari\n", "\n", "Last checked: June 2021\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 1, "id": "6c8552a0-52f2-48c2-95e7-b4f43152f05b", "metadata": {}, "outputs": [], "source": [ "import pyfvtool as pf\n", "import matplotlib.pyplot as plt\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "id": "ce01d9b3-03cb-444e-8347-626bd39dc03e", "metadata": {}, "outputs": [], "source": [ "## define the domain\n", "L = 1.0 # domain length\n", "Nx = 25 # number of cells" ] }, { "cell_type": "code", "execution_count": 3, "id": "2a37d6ab-7bf9-4c6c-9715-c35993553798", "metadata": {}, "outputs": [], "source": [ "meshstruct = pf.Grid1D(Nx, L)" ] }, { "cell_type": "code", "execution_count": 4, "id": "a4ce839a-124f-4c17-a71d-91f4ea6fa346", "metadata": {}, "outputs": [], "source": [ "c_central = pf.CellVariable(meshstruct, 0.0)" ] }, { "cell_type": "code", "execution_count": 5, "id": "8e52b20d-6f7d-4cd3-b406-7fa248afa0b4", "metadata": {}, "outputs": [], "source": [ "# Set boundary conditions\n", "c_central.BCs.left.a = 0.0 \n", "c_central.BCs.left.b = 1.0 \n", "c_central.BCs.left.c = 0.0 # left boundary\n", "c_central.BCs.right.a = 0.0 \n", "c_central.BCs.right.b = 1.0 \n", "c_central.BCs.right.c = 1.0 # right boundary" ] }, { "cell_type": "code", "execution_count": 6, "id": "d00aba31-44ff-4dac-bd55-33af030eb419", "metadata": {}, "outputs": [], "source": [ "c_upwind = c_central.copy() # separate copy for comparison" ] }, { "cell_type": "code", "execution_count": 7, "id": "c769b389-f284-4b71-8141-25de33800360", "metadata": {}, "outputs": [], "source": [ "x = meshstruct.cellcenters.x" ] }, { "cell_type": "code", "execution_count": 8, "id": "08083985-fe24-4126-b11a-5c0850a32f6d", "metadata": {}, "outputs": [], "source": [ "## define the transfer coeffs\n", "D_val = 1.0\n", "D = pf.CellVariable(meshstruct, D_val)\n", "Dave = pf.harmonicMean(D) # convert a cell variable to face variable" ] }, { "cell_type": "code", "execution_count": 9, "id": "36427ec8-f484-4bc8-988d-5da9abee7ed8", "metadata": {}, "outputs": [], "source": [ "alfa = pf.CellVariable(meshstruct, 1.0)" ] }, { "cell_type": "code", "execution_count": 10, "id": "2e895c26-b2ab-4a9c-b847-771300658771", "metadata": {}, "outputs": [], "source": [ "u = 10.0\n", "u_face = pf.FaceVariable(meshstruct, u)" ] }, { "cell_type": "code", "execution_count": 11, "id": "0eca751e-b0d3-4d37-b538-b98668afd107", "metadata": {}, "outputs": [], "source": [ "## solve\n", "pf.solvePDE(c_central, [ pf.convectionTerm(u_face),\n", " -pf.diffusionTerm(Dave)])\n", "\n", "pf.solvePDE(c_upwind, [ pf.convectionUpwindTerm(u_face),\n", " -pf.diffusionTerm(Dave)]);" ] }, { "cell_type": "code", "execution_count": 12, "id": "05b0dc00-043c-41da-92bd-700925b9a54f", "metadata": {}, "outputs": [], "source": [ "c_analytic = (1-np.exp(u*x/D_val))/(1-np.exp(u*L/D_val))" ] }, { "cell_type": "code", "execution_count": 13, "id": "7cac7a1b-1702-4450-bc81-0737f529407b", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABYT0lEQVR4nO3deXxU1cH/8c+dmUw2SEICWYBA2PdFQDQguGMVcamKW1VcWilaRao+8tinLj/7YK0LlApqC6KWItWqtZVHRCtrcGFTTFgTQoAkxCSQhOyZub8/JhMSkkAmJJnM5Pt+veY1d+5y5uQGvd+cc+65hmmaJiIiIiJeYvF2BURERKRjUxgRERERr1IYEREREa9SGBERERGvUhgRERERr1IYEREREa9SGBERERGvUhgRERERr7J5uwJN4XQ6yczMpHPnzhiG4e3qiIiISBOYpklRURHdu3fHYmm8/cMnwkhmZibx8fHeroaIiIg0w6FDh+jZs2ej230ijHTu3Blw/TBhYWFero2IiIg0RWFhIfHx8TXX8cb4RBhxd82EhYUpjIiIiPiYMw2x0ABWERER8SqFEREREfEqhRERERHxKp8YM9IUpmlSVVWFw+HwdlXkFFarFZvNptuyRUSkQX4RRioqKsjKyqKkpMTbVZFGhISEEBcXh91u93ZVRESknfH5MOJ0Ojlw4ABWq5Xu3btjt9v1F3g7YpomFRUV/Pjjjxw4cIABAwacduIbERHpeHw+jFRUVOB0OomPjyckJMTb1ZEGBAcHExAQwMGDB6moqCAoKMjbVRIRkXbEb/5E1V/b7Zt+PyIi0hhdIURERMSrFEZERETEqxRGxGMJCQnMnz/f29UQERE/oTDiZdnZ2fzqV7+ib9++BAYGEh8fz7Rp0/jiiy9a7DsuuugiZs+e3WLliYiItCSfv5vGl6WnpzNx4kQiIiJ44YUXGDlyJJWVlaxevZoHHniA3bt3t1ldTNPE4XBgs+mfhIhIh/LNn+FYOoy4Ebqf45Uq+F3LiGmalFRUeeVlmqZHdZ01axaGYfDNN99w4403MnDgQIYNG8acOXP46quvACgoKOAXv/gF0dHRhIWFcckll/Ddd9/VlPH0008zevRo3nnnHRISEggPD+eWW26hqKgIgBkzZrBu3ToWLFiAYRgYhkF6ejpr167FMAxWr17NuHHjCAwMZMOGDaSmpnLttdcSExNDp06dOPfcc/n8889b7hckIiLty873YfOfIHef16rgd38Gl1Y6GPrb1V757pRnryDE3rRTmp+fz6effsrvfvc7QkND622PiIjANE2mTp1KZGQkq1atIjw8nNdff51LL72UvXv3EhkZCUBqaiofffQR//73vzl27BjTp0/n+eef53e/+x0LFixg7969DB8+nGeffRaAbt26kZ6eDsDjjz/Oiy++SN++fYmIiODw4cNcddVVPPfccwQFBfHWW28xbdo09uzZQ69evVrmRImISPvgdED2Ttdy3CivVcPvWkZ8xf79+zFNk8GDBze6z5dffsnOnTt57733GDduHAMGDODFF18kIiKC999/v2Y/p9PJsmXLGD58OJMmTeKOO+6oGXMSHh6O3W4nJCSE2NhYYmNjsVqtNcc+++yzXH755fTr14+oqChGjRrF/fffz4gRIxgwYADPPfccffv25eOPP269kyEiIt6Rl8pmaxXX9uzO5opcr1XD71pGggOspDx7hde+u6ncXTqnm7p+69atnDhxgqioqDrrS0tLSU1NrfmckJBA586daz7HxcWRk5PTpHqMGzeuzufi4mKeeeYZ/v3vf5OZmUlVVRWlpaVkZGQ0qTwREfEdZuYOFkRGkBZgY8H2hZzffYJXHqnid2HEMIwmd5V404ABAzAMg127dnHdddc1uI/T6SQuLo61a9fW2xYREVGzHBAQUGebYRg4nc4m1ePULqLHHnuM1atX8+KLL9K/f3+Cg4O58cYbqaioaFJ5IiLiO5IOriE5MBCA5LxkkjKTmNhjYpvXo/1ftf1UZGQkV1xxBa+++ioPPfRQvVBw/PhxxowZQ3Z2NjabjYSEhGZ/l91ux+FwNGnfDRs2MGPGDK6//noATpw4UTO+RERE/IdpmizM34YFE6dhYDEsLNy+kAleaB3RmBEvWrRoEQ6Hg/Hjx/OPf/yDffv2sWvXLv74xz+SmJjIZZddRmJiItdddx2rV68mPT2dpKQkfvOb37Bly5Ymf09CQgJff/016enp5ObmnrbVpH///nzwwQfs2LGD7777jttuu63JrSwiIuI7ko5sItmowFkdPJyms6Z1pK0pjHhRnz592LZtGxdffDG//vWvGT58OJdffjlffPEFixcvxjAMVq1axeTJk7nnnnsYOHAgt9xyC+np6cTExDT5ex599FGsVitDhw6lW7dupx3/8corr9ClSxcmTJjAtGnTuOKKKxgzZkxL/LgiItJOmKbJwi0vYTllSgp364inU1WcLcNs629shsLCQsLDwykoKCAsLKzOtrKyMg4cOECfPn30aPp2TL8nEZH2Y9ORTcz8fGaj21+77LUWGTtyuut3bWoZERER6UBM02Th9oU0NirEwGjz1pFmhZFFixbV/IU7duxYNmzYcNr9ly9fzqhRowgJCSEuLo67776bvLy8ZlVYREREmq/SWUl2cTaNRQ0Tk+zibCqdlW1WJ4/vplm5ciWzZ89m0aJFTJw4kddff50rr7ySlJSUBmfo3LhxI3feeSevvPIK06ZN48iRI8ycOZP77ruPDz/8sEV+CBEREWkau9XOu1NXkP/GJCg7Dte/AdF1J+CMDIrEbrW3WZ08bhl5+eWXuffee7nvvvsYMmQI8+fPJz4+nsWLFze4/1dffUVCQgIPPfQQffr04YILLuD+++/36G4QERERaTmxDgdDC39kaJVJcMRkwq19GBI5hKFRQxkaNZTY0Ng2rY9HYaSiooKtW7cyZcqUOuunTJlCUlLDtwJNmDCBw4cPs2rVKkzT5OjRo7z//vtMnTq10e8pLy+nsLCwzktERERaSFb1A1e7DeF3n6Uz8fn/8Pbmg16rjkdhJDc3F4fDUe+20piYGLKzsxs8ZsKECSxfvpybb74Zu91ObGwsERERLFy4sNHvmTdvHuHh4TWv+Ph4T6opIiIip+MOI3Gj2HvU9ZT3ATGdvFadZg1gPXVmNtM0G52tLSUlhYceeojf/va3bN26lU8//ZQDBw4wc2bjtxTNnTuXgoKCmtehQ4eaU00RERFpSHUYKe82nIz8EgAGxXQ+3RGtyqMBrF27dsVqtdZrBcnJyWl0Eq558+YxceJEHnvsMQBGjhxJaGgokyZN4rnnniMuLq7eMYGBgQRWz5UvIiIiLaw6jBy0DwDK6dopkKhO3rvuetQyYrfbGTt2LGvWrKmzfs2aNUyYMKHBY0pKSrBY6n6N+xH2PjDfml9Zu3YthmFw/PjxsypnxowZjT7cT0RE2rmio1CUBRh8X9UTgMGx3msVgWZ008yZM4e//OUvLF26lF27dvHII4+QkZFR0+0yd+5c7rzzzpr9p02bxgcffMDixYtJS0tj06ZNPPTQQ4wfP57u3bu33E/iYy666CJmz55db/1HH33Uag8omjBhAllZWYSHh7dK+SIi4gPc40W6DiQ51/UQ1UFeDiMezzNy8803k5eXx7PPPktWVhbDhw9n1apV9O7dG4CsrKw6zz6ZMWMGRUVF/OlPf+LXv/41ERERXHLJJfz+979vuZ9CmsQ9gFhERDqwWoNX92S7Bq96c7wINHMA66xZs0hPT6e8vJytW7cyefLkmm3Lli1j7dq1dfb/1a9+RXJyMiUlJWRmZvLXv/6VHj16nFXFO4Knn36a0aNH8/rrrxMfH09ISAg33XRTTTfLzp07sVgs5ObmAnDs2DEsFgs33XRTTRnz5s0jMTERqN9Ns2zZMiIiIli9ejVDhgyhU6dO/OQnPyErK6vmeIfDwZw5c4iIiCAqKorHH39c3WsiIr4sa4frvfvok2HE17pp2j3ThIpi77xa4SK9f/9+/v73v/Ovf/2LTz/9lB07dvDAAw8AMHz4cKKioli3bh0A69evJyoqivXr19ccv3btWi688MJGyy8pKeHFF1/knXfeYf369WRkZPDoo4/WbH/ppZdYunQpS5YsYePGjeTn52vmXBERX5b1PQDHI4aSV1yBYXj3tl5oRjdNu1dZAv/rpbEo/50J9tAWLbKsrIy33nqLnj1dg4wWLlzI1KlTeemll4iNjWXy5MmsXbuWG264gbVr13LXXXfx1ltvkZKSwsCBA0lKSuKRRx5ptPzKykpee+01+vXrB8CDDz7Is88+W7N9/vz5zJ07lxtuuAGA1157jdWrV7fozygiIm2kJB8KXEMpdpsJQAq9IkMIsXs3Dvhfy4if6dWrV00QAUhMTMTpdLJnzx7ANRDW3S22bt06Lr74YiZPnsy6dev49ttvKS0tZeLExh8DHRISUhNEAOLi4sjJyQGgoKCArKysmm4eAJvNxrhx41ryRxQRkbbiHi8S2ZfkfNeit8eLgD+2jASEuFoovPXdTRQWFkZBQUG99cePHycsLKzR49x32rjfL7roIh5++GH279/PDz/8wKRJk0hNTWXdunUcP36csWPH0rlz4//QAgIC6pWvMSEiIn6q9syr1eNFvH1bL/hjy4hhuLpKvPHy4JbcwYMHN/iwwG+//ZZBgwbVfM7IyCAz82S42rx5MxaLhYEDBwInx40899xzjBo1irCwMC688ELWrVt3xvEiZxIeHk5cXBxfffVVzbqqqiq2bt3a7DJFRMSLaoWR3Ufdg1cb/wO4rfhfGPERs2bNIjU1lQceeIDvvvuOvXv38uqrr7JkyZKa2WoBgoKCuOuuu/juu+/YsGEDDz30ENOnT6+5RdcwDCZPnsxf//pXLrroIsA1y21FRQVffPFFzbrmevjhh3n++ef58MMP2b17N7NmzTrrSdNERMRLqsOIM2YU+2rCiHcHr4LCiNckJCSwYcMGUlNTmTJlCueeey7Lli1j2bJldW7N7d+/Pz/96U+56qqrmDJlCsOHD2fRokV1yrr44otxOBw1wcMwDCZNmgTABRdccFb1/PWvf82dd97JjBkzSExMpHPnzlx//fVnVaaIiHhBWSHkpwJwJLg/JRUO7DYLCVEte+NFcximDwwQKCwsJDw8nIKCgnrjKcrKyjhw4AB9+vQhKCjISzVsHU8//TQfffQRO3bs8HZVzpo//55ERHxC+iZYdhWEx/PZlM/5xTtbGRoXxqqHJ7XaV57u+l2bWkZEREQ6ggZmXm0Pg1dBYURERKRjcM+8Wmvw6kCFETmTp59+2i+6aEREpB1o6Jk0CiMiIiLSJiqKIXcvAOXRIziQWwyom0ZERETaytFkMJ3QKZbUkk44nCZhQTZiw9rHDQUKIyIiIv6udhfN0ULA1UVjeDBZZ2tSGBEREfF3tQav7sk+AbSf8SKgMCIiIuL/6gxedbeMeH8aeDeFEREREX9WVQ45u1zLte+kaQdP63VTGPFzCQkJzJ8//6zKWLt2LYZh6Jk0IiK+KCcFnFUQHEmBPYbMgjJAYURqSUpKwmq18pOf/MTbVQHgoosuYvbs2XXWTZgwgaysLMLDw71TKRERab5aXTT7clzjReLCgwgPCfBipepSGKllc+Zmrv3oWjZnbm6z71y6dCm/+tWv2LhxIxkZGW32vZ6w2+3Exsa2m1HXIiLigVphZHc7m+zMTWGkmmmaLNi2gLSCNBZsW0BbPD+wuLiYv//97/zyl7/k6quvZtmyZTXb3F0jX3zxBePGjSMkJIQJEyawZ8+emn1SU1O59tpriYmJoVOnTpx77rl8/vnnjX7fPffcw9VXX11nXVVVFbGxsSxdupQZM2awbt06FixYgGEYGIZBenp6g900mzZt4sILLyQkJIQuXbpwxRVXcOzYsRY7NyIi0kIyd7je2+l4EVAYqZGUmURyXjIAyXnJJGUmtfp3rly5kkGDBjFo0CB+9rOf8eabb9YLQU8++SQvvfQSW7ZswWazcc8999RsO3HiBFdddRWff/4527dv54orrmDatGmNtrDcd999fPrpp2RlZdWsW7VqFSdOnGD69OksWLCAxMREfv7zn5OVlUVWVhbx8fH1ytmxYweXXnopw4YNY/PmzWzcuJFp06bhcDha6MyIiEiLcFS6JjyD6jlG1DLSbpmmycLtC7EYrtNhMSws3L6w1VtHlixZws9+9jMAfvKTn3DixAm++OKLOvv87ne/48ILL2To0KE88cQTJCUlUVbmGnw0atQo7r//fkaMGMGAAQN47rnn6Nu3Lx9//HGD3zdhwgQGDRrEO++8U7PuzTff5KabbqJTp06Eh4djt9sJCQkhNjaW2NhYrFZrvXJeeOEFxo0bx6JFixg1ahTDhg3jwQcfpGvXri11akREpCX8uAcc5RAYjtmlT7t7Jo2bwggnW0WcphMAp+ls9daRPXv28M0333DLLbcAYLPZuPnmm1m6dGmd/UaOHFmzHBcXB0BOTg7g6uZ5/PHHGTp0KBEREXTq1Indu3efduzJfffdx5tvvllTzieffFKntaUp3C0jIiLSztWMFxnJ0aIKCkorsVoM+nXr5N16ncLm7Qp4W+1WEXcYgZOtIxO6T2iVgZtLliyhqqqKHj161KlLQEBAnbEXAQEnRzu76+F0uur52GOPsXr1al588UX69+9PcHAwN954IxUVFY1+75133skTTzzB5s2b2bx5MwkJCUyaNMmjugcHB3u0v4iIeEmdwauuyc4SokIICqjf6u1NHb5l5NRWEbfWbB2pqqri7bff5qWXXmLHjh01r++++47evXuzfPnyJpWzYcMGZsyYwfXXX8+IESOIjY0lPT39tMdERUVx3XXX8eabb/Lmm29y991319lut9vPOPZj5MiR9bqTRESkHaoVRvZWjxcZ3I5mXnXr0GHE3Spi0HDLh4HRKmNH/v3vf3Ps2DHuvfdehg8fXud14403smTJkiaV079/fz744IOaIHPbbbfVtJqczn333cdbb73Frl27uOuuu+psS0hI4OuvvyY9PZ3c3NwGy5s7dy7ffvsts2bN4vvvv2f37t0sXryY3Nzcpp0AERFpfU4HZO90Lbfj23qhg4eRSmcl2cXZmDQcNkxMsouzqXRWtuj3LlmyhMsuu6zBScRuuOEGduzYwbZt285YziuvvEKXLl2YMGEC06ZN44orrmDMmDFnPO6yyy4jLi6OK664gu7du9fZ9uijj2K1Whk6dCjdunVrcPzJwIED+eyzz/juu+8YP348iYmJ/POf/8Rm6/C9fiIi7UdeKlQWQ0AIRPVvt4NXAQyzLSbUOEuFhYWEh4dTUFBAWFjd5qWysjIOHDhAnz59CAoK8rjs7OJs8svyG90eGRRJbGisx+W2ZyUlJXTv3p2lS5fy05/+tE2+82x/TyIi4qHv34MP7oP486ia8SlDn1pNRZWTtY9eRELX0Dapwumu37V1+D9lY0Nj/S5sNMbpdJKdnc1LL71EeHg411xzjberJCIirSVrh+s9bhQH80uoqHISHGClV2SIV6vVkA4fRjqSjIwM+vTpQ8+ePVm2bJm6VURE/FmtwavuLpqBMZ2wWNrfoz10NepAEhIS2mSaexER8TLThKzvXctxo9i90x1G2t94EejgA1hFRET80rEDUF4A1kDoNpi97XjwKiiMiIiI+B93F03MMLAG1DyTpj3OMQJ+FEbU/dC+6fcjItKGao0XKa1wkJ5XDKhlpNW4p0svKSnxck3kdNy/n9rT24uISCupFUb25RRhmhAZaqdrJ7t369UInx/AarVaiYiIqHl4XEhISKs8S0aaxzRNSkpKyMnJISIiosGnAIuISAsyzbp30mRVjxeJ6dxur48+H0YAYmNd84S4A4m0PxERETW/JxERaUWFR6AkDyw2iB7Knu1pQPvtogE/CSOGYRAXF0d0dDSVlS07dbucvYCAALWIiIi0FXerSLchEBBUa/CqwkibsFqtuuiJiEjHVquLBjg54Vk7DiM+P4BVREREaqkVRo4VV5BTVA603wnPQGFERETEv9QKI7urW0XiI4PpFNh+O0MURkRERPxF0VEoygIMiB3OnuxCwHUnTXumMCIiIuIv3K0iXQeCPZQ9R08A7ftOGlAYERER8R/uMNJ9NMDJlpF2Og28m8KIiIiIv8ja4XqPG4Vpmuytbhlpz7f1gsKIiIiI/8j63vUeN4rDx0o5UV5FgNWgT9dQ79brDBRGRERE/EFJPhRkuJZjR7C3erKzft06EWBt35f79l07ERERaRr3eJHIvhAUXnNbb3sfvAoKIyIiIv6hkZlXFUZERESkbZwSRtzdNO19jhFQGBEREfEPtcJIpcNJ6o++MccIKIyIiIj4vrJCyE91LceOIu3HYiodJp0DbfSICPZu3ZpAYURERMTXZe90vYfHQ2gUu6snOxsY2xnDMLxYsaZRGBEREfF1jYwXac9P6q1NYURERMTX1cy8Oho4eSdNe5951U1hRERExNed0jLiS3OMgMKIiIiIb6sohty9ruW4UZwor+LwsVLAN27rBYURERER33Y0GUwndIqFzjE140WiOwfSJdTu5co1jcKIiIiIL/PhmVfdFEZERER8Wc3g1bphxFcGr4LCiIiIiG9rpGXEV27rBYURERER31VVDjm7XMtxozBNkz1H3S0jYV6smGcURkRERHxVTgo4qyA4EsJ78uOJcvKLK7AYMCCmk7dr12QKIyIiIr6qdheNYdR00SREhRIUYPVixTyjMCIiIuKr/GC8CCiMiIiI+K7MHa737qMB37ytFxRGREREfJOj0jXhGZxsGTnqe7f1gsKIiIiIb/pxDzjKITAcuvTB6TRrZl/tEC0jixYtok+fPgQFBTF27Fg2bNhw2v3Ly8t58skn6d27N4GBgfTr14+lS5c2q8IiIiJCrfEiI8EwyMgvoazSSaDNQu+oUO/WzUM2Tw9YuXIls2fPZtGiRUycOJHXX3+dK6+8kpSUFHr16tXgMdOnT+fo0aMsWbKE/v37k5OTQ1VV1VlXXkREpMNq5Em9A2I6YbUY3qpVs3gcRl5++WXuvfde7rvvPgDmz5/P6tWrWbx4MfPmzau3/6effsq6detIS0sjMjISgISEhLOrtYiISEfX2DNpYnxnsjM3j7ppKioq2Lp1K1OmTKmzfsqUKSQlJTV4zMcff8y4ceN44YUX6NGjBwMHDuTRRx+ltLS0+bUWERHpyJwOyN7pWq4OIyfHi/jOZGduHrWM5Obm4nA4iImJqbM+JiaG7OzsBo9JS0tj48aNBAUF8eGHH5Kbm8usWbPIz89vdNxIeXk55eXlNZ8LCws9qaaIiIh/y0uFymIICIGo/gDsznZdKwf50DTwbs0awGoYdfuiTNOst87N6XRiGAbLly9n/PjxXHXVVbz88sssW7as0daRefPmER4eXvOKj49vTjVFRET8k7uLJnYEWKyUVTpIzysBfO+2XvAwjHTt2hWr1VqvFSQnJ6dea4lbXFwcPXr0IDw8vGbdkCFDME2Tw4cPN3jM3LlzKSgoqHkdOnTIk2qKiIj4t8Pfut6ru2j255zA4TSJCAkgunOgFyvWPB6FEbvdztixY1mzZk2d9WvWrGHChAkNHjNx4kQyMzM5ceJEzbq9e/disVjo2bNng8cEBgYSFhZW5yUiIiLV0r50vfeZDJwcLzIwpnOjPRXtmcfdNHPmzOEvf/kLS5cuZdeuXTzyyCNkZGQwc+ZMwNWqceedd9bsf9tttxEVFcXdd99NSkoK69ev57HHHuOee+4hODi45X4SERGRjuD4IcjdC4YVEiYBJ++k8cUuGmjGrb0333wzeXl5PPvss2RlZTF8+HBWrVpF7969AcjKyiIjI6Nm/06dOrFmzRp+9atfMW7cOKKiopg+fTrPPfdcy/0UIiIiHYW7VaTnOAiOAE7OMeJrM6+6eRxGAGbNmsWsWbMa3LZs2bJ66wYPHlyva0dERESaYf8Xrvd+l9Ss2uujz6Rx07NpREREfIXTAWlrXcvVYaSgpJKsgjIABsQojIiIiEhrytwBZcddD8frPgY4+aTeHhHBhAUFeK9uZ0FhRERExFek/sf13ncyWF0jLfbUTHbmm60ioDAiIiLiO9xhpNZ4kT1HfXvwKiiMiIiI+IayQjj8jWu5dhipeUCewoiIiIi0pvSN4KyCyH7QJQFwPY7F12/rBYURERER39BAF012YRlFZVXYLAb9uvne03rdFEZERER8QWr9+UXcrSJ9uoZit/nuJd13ay4iItJR5B+A/DSw2CDhgprVe/ygiwYURkRERNq/mingx0PQyYfH+vozadwURkRERNq7BsaLbM7czNrix7CG7GNQrG8/3V5hREREpD1zVEHaetdydRgxTZMF2xZQZc0mMHo1A6N9d/AqKIyIiIi0b5nboLwAgrtA99EAJGUmkZyXDIA1+DAZZdu9WMGzpzAiIiLSntVMAX8RWKyYpsnC7Qsx3Jdw0+DVHX/CNE2vVfFsKYyIiIi0Z6eMF3G3ipg4XesNk+S8ZJIyk7xUwbOnMCIiItJelR6Hw1tcy30vrmkVsRh1L98Ww8LC7Qt9tnVEYURERKS9OrAeTAd0HQgR8TWtIk7TWWc3p+n06dYRhREREZH2qlYXzcmxIkaDuxoYPts6ojAiIiLSHplmnSngK52VZBdnY9Jw2DAxyS7OptJZ2YaVbBk2b1dAREREGpCfBsczwBIAvSdit9p59+p3yS/LZ03yURb8Zx+jeobzu+tH1BwSGRSJ3Wr3YqWbR2FERESkPXJ30fQ6HwJdk5rFhsYSGxrLinwnzrISxsT2YWjUUC9WsmWom0ZERKQ9Sq1+Hk2tKeDdNqXmAjC+T2Rb1qjVKIyIiIi0N45K1500UC+MZBeUkfZjMRYDzu8b5YXKtTyFERERkfbm8LdQUQQhURA7ss6mTftdrSIjekYQHhzgjdq1OIURERGR9qZmCviLwVL3Uu3uopnYzz9aRUBhREREpP05ZQp4N9M0a1pGJvbv2ta1ajUKIyIiIu1JST4c2eZa7ndxnU2pPxZztLAcu83C2N5dvFC51qEwIiIi0p4cWAeY0G0IhHWvsympuotmXO8uBAVYvVC51qEwIiIi0p400kUD+GUXDSiMiIiItB+meXJ+kf51w4jDabI5NQ+ACX40eBUURkRERNqPvP1QcAisgdBrQp1NyZkFFJZV0TnQxoge4V6qYOtQGBEREWkv9lc/GK93IthD6mzatN/VKnJe3yhsVv+6fPvXTyMiIuLLTjNexD14dWJ//+qiAYURERGR9qGqHNI3uJZPCSPlVQ6+Tc8H/G/wKiiMiIiItA+HvoHKEgiNhuhhdTZtO3icskon3ToHMiC6k5cq2HoURkRERNqDmi6a+lPAJ9WaAt4wjLauWatTGBEREWkPTjNeZGP1/CIT/LCLBhRGREREvK84F7K+cy33rTsFfFFZJd8fLgD8c7wIKIyIiIh4X9pawISYEdA5ps6mr9PycThNEqJC6BER7JXqtTaFEREREW+rPV7kFJtS/buLBhRGREREvMs0Tz+/SPVkZxP7KYyIiIhIa/hxNxRlgS0IeiXW3VRUzp6jRQAk+tnzaGpTGBEREfEmd6tI74kQEFRnk/uW3qFxYUSG2tu6Zm1GYURERMSbmtJF44dTwNemMCIiIuItlWWQvsm13EAY6QiDV0FhRERExHsOfQVVpdA5DqKH1NmUkVfC4WOlBFgNxidEeqmCbUNhRERExFtqd9GcMs27e9bVc+K7EBpoa+uatSmFEREREW/Z3/h4kZNdNP49XgQURkRERLyj6Cgc3ela7ntRnU1Op8nmVPfgVf8eLwIKIyIiIt6Rttb1HjcKQusGjt3ZReQXVxBitzKqZ0SbV62tKYyIiIh4w+lu6a3uohnfJxK7zf8v1f7/E4qIiLQ3Z5gCflP14FV/ngK+NoURERGRtnY0GYpzICAE4s+rs6nS4eSbA/lAxxi8CgojIiIibc/dKpJwAdgC62z67tBxiiscRIbaGRIb5oXKtT2FERERkbaW+oXrvd+l9TZtqp4CPrFvFBaLUW+7P1IYERERaUsVJXBws2u5g88v4qYwIiIi0pYyksBRDmE9oeuAOptKKqrYnnEM6DiDV0FhREREpG2lful673dxvSngvzmQT6XDpEdEML2jQrxQOe9QGBEREWlLp51fxD3rahSG0THGi4DCiIiISNspzIKcFMCoNwU81JpfpANMAV+bwoiIiEhbSavuoul+DoRE1tl0rLiClKxCABL7dZzBq6AwIiIi0nZO00WzOS0P04SBMZ2I7hzUxhXzLoURERGRtuB0ngwj/RuaX6T6lt4OdBeNm8KIiIhIW8j+DkrywN4Jep5bb/PJwasKIyIiItIadr7veu93MVgD6mzKPF7KgdxiLAac1zeygYP9m8KIiIhIa3NUwnfvupZH315vs7uLZmTPCMKCAupt93cKIyIiIq1t76dQkgudYqD/5fU2155fpCNSGBEREWlt2//qeh91C1htdTaZpslG9/wiHXDwKiiMiIiItK7CLNj3mWv5nDvqbd6fc4Ifi8oJtFkY07tLG1eufVAYERERaU3fvwumE+LPr/dgPDg5XuTchEiCAqxtXbt2oVlhZNGiRfTp04egoCDGjh3Lhg0bmnTcpk2bsNlsjB49ujlfKyIi4ltM82QXzTk/a3CXTdXjRSZ00PEi0IwwsnLlSmbPns2TTz7J9u3bmTRpEldeeSUZGRmnPa6goIA777yTSy+tP9GLiIiIX8r4CvL2Q0AoDLuu3uYqh5Ov0qoHr3bQ8SLQjDDy8ssvc++993LfffcxZMgQ5s+fT3x8PIsXLz7tcffffz+33XYbiYmJza6siIiIT3G3igy7HgI719v8Q2YhRWVVhAXZGN4jvI0r1354FEYqKirYunUrU6ZMqbN+ypQpJCUlNXrcm2++SWpqKk899VSTvqe8vJzCwsI6LxEREZ9SXgTJH7qWG+uiqR4vcn7fKKwWo61q1u54FEZyc3NxOBzExMTUWR8TE0N2dnaDx+zbt48nnniC5cuXY7PZGtznVPPmzSM8PLzmFR8f70k1RUREvC/5I6gshqj+0Ov8BndJSq2+pbcDTgFfW7MGsBpG3fRmmma9dQAOh4PbbruNZ555hoEDBza5/Llz51JQUFDzOnToUHOqKSIi4j21B642cI0sq3SwJf0Y0HEnO3NrWlNFta5du2K1Wuu1guTk5NRrLQEoKipiy5YtbN++nQcffBAAp9OJaZrYbDY+++wzLrmk/mOUAwMDCQwM9KRqIiIi7UfuPjj0FRhWGHVrg7tsO3iM8ion0Z0D6detUxtXsH3xqGXEbrczduxY1qxZU2f9mjVrmDBhQr39w8LC2LlzJzt27Kh5zZw5k0GDBrFjxw7OO++8s6u9iIhIe7T9Hdf7gMuhc2yDu2yq1UXTUO9CR+JRywjAnDlzuOOOOxg3bhyJiYm88cYbZGRkMHPmTMDVxXLkyBHefvttLBYLw4cPr3N8dHQ0QUFB9daLiIj4BUcl7FjhWm5k4CrAxv3V84v069hdNNCMMHLzzTeTl5fHs88+S1ZWFsOHD2fVqlX07t0bgKysrDPOOSIiIuK39n8OxTkQ0hUGXNHgLgWllew8fBzQ4FUAwzRN09uVOJPCwkLCw8MpKCggLCzM29URERFp3Lu3w+5/Q+KDcMXvGtzls+RsfvHOVvp2DeU/j17UtvVrQ029fuvZNCIiIi3lRA7s/dS1fJoumiRNAV+HwoiIiEhL+X4lOKugxziIHtLobu7JzjryFPC1KYyIiIi0BNOEbdV30TTSKrI5czNTP7iGtBPbMQxI1OBVQGFERESkZRzeArl7wBYMw39ab7NpmizYtoCMogMERq9maPfORITYvVDR9sfju2lERESkAe65RYZeC0H1H3qXlJlEcl4yANbgw/SNPdKWtWvX1DIiIiJytiqK4YcPXMtj7qi32TRNFm5fiMWwVH82SKv6Bz5wQ2ubUBgRERE5WykfQ0URdOkDvSfW2+xuFXGaTgAMwySjeA9JmY0/8b4jURgRERE5WzUPxbu93kPxTm0VcbMYFhZuX6jWERRGREREzk5eKhzcCBgw6rZ6m09tFXFzmk6S85LVOoLCiIiIyNnZsdz13v9SCO9RZ5O7VcSg4QfhGRhqHUFhREREpPmcDtjxN9dyA3OLVDoryS7OxqThsGFikl2cTaWzsjVr2e7p1l4REZHmSv0PFGVBcCQMuqreZrvVzrtXv0t+WT5LNx7gg+1HOL9PJL+5emjNPpFBkditHXu+EYURERGR5nLPLTLyZrAFNrhLbGgskYHRfPn9EZxlPbj73HEMjYppw0q2f+qmERERaY7iXNi9yrV8zu2n3fXzXUfJK64gJiyQiwd1a4PK+RaFERERkeb4/u/grIS40RA74rS7rvgmA4CbxsZjs+rSeyqdEREREU+ZZq25RRp+KJ7bofwSNlY/pXf6uPjWrplPUhgRERHxVOZ2yEkGWxCMuOm0u7635RCmCRf070qvqJA2qqBvURgRERHxlLtVZMg0CI5odLcqh5O/bzkMwC3j1SrSGIURERERT1SWws73Xctn6KJZt/dHsgvL6BISwOVDdQdNYxRGREREPLHrX1BeAOG9IGHyaXd999tDANwwpieBNmtb1M4nKYyIiIh4wj23yDm3g6Xxy2hOYRn/2Z0DqIvmTBRGREREmupYOhxYDxgwuv5D8Wp7b+thHE6Tcb270D+6c5tUz1cpjIiIiDSV+zk0fS+CiF6N7uZ0mqys7qK5ZXzj+4mLwoiIiEhTOB2wvfoJvWcYuLo5LY+M/BI6B9q4akRsG1TOtymMiIiINMWBdVB4GILCYfDVp93VPXD12nO6E2LXY+DORGFERESkKbZVD1wdMR0CghrdLb+4gtU/ZANwy7nqomkKhREREZEzKcmH3f92LZ+hi+aDbYepcDgZ3iOM4T3C26Byvk9hRERE5Ex2vg+OCogZAXGjGt3NNGsNXFWrSJMpjIiIiJyJe26RMXeAYTS627aMY+zLOUFwgJVrR3dvo8r5PoURERGR08n6DrK/B6v9jA/FW/GNq1Vk6sg4OgcFtEXt/ILCiIiIyOm4b+cdPBVCIhvdrbCskk++zwLgVs246hGFERERkcZUlMD3K13LZxi4+vGOTEorHQyI7sSYXl3aoHL+Q2FERESkMV+/BmXHIaI39L34tLu++20GADefG49xmnElUp/CiIiISENK8mHjfNfyxU+CpfGn7v5wpIAfjhRit1r46ZiebVM/P6IwIiIi0pCNL0N5AcQMP+PAVXeryBXDY4kMtbdF7fyKwoiIiMipjh+Cr99wLV/2NFgav1yWVjj45/ZMAG45VwNXm0NhRERE5FRrnwdHOfS+APpfdtpdP9mZRVF5Fb0iQ0jsG9VGFfQvCiMiIiK15eyC7/7mWr78mdNOcgbw7jcnB65aLBq42hwKIyIiIrV98SyYThgyDXqOO+2u+3OK2HLwGFaLwU1jNXC1uRRGRERE3DK+gj2rwLDCpU+dcfd3q2dcvWRwNNFhjT/JV05PYURERATANGFNdQA552fQdcBpdy+vcvCPbYcBDVw9WwojIiIiAHs/hUNfgS0ILnrijLuvSTnKsZJKYsOCuHBgtzaooP9SGBEREXE64PNnXMvn/xLCzvzEXXcXzfRxPbFZdTk9Gzp7IiIi370LP+6CoAiYOPuMu2fklbBxfy6GATeNUxfN2VIYERGRjq2yDL78X9fypF9DcMQZD/n7FleryAX9uxIfGdKKlesYFEZERKRj+/bPUHgYwnrA+F+ccfcqh5P3trrCyK3je7V27ToEhREREem4So/Dhpdcyxf/NwSc+fbctXt+5GhhOVGhdi4bEtO69esgFEZERKTj2rQASo9Bt8Ew6tYmHeJ+KN4NY3tit+ky2hJ0FkVEpGMqzIKvFruWL30KLNYzHpJdUMZ/ducArunfpWUojIiISMe07nmoKoX482HQlU065P2th3CaMD4hkn7dOrVyBTsOhREREel4cvfBtndcy5c9fcaH4QE4nSYrq++iuWW8WkVaksKIiIh0PF88C6YDBl4JvRNPu+vmzM1c+9G1LNm6mkP5pXQOsnHl8Lg2qmjHoDAiIiIdy+GtsOtjwIBLf3vaXU3TZMG2BaQVpLE0ZRFgcv05PQi2n3l8iTSdwoiIiHQcpgmfVz8Mb/RtEDP0tLsnZSaRnJcMwAkOYA3dxy3nam6RlqYwIiIiHcf+LyB9A1gD4aK5p93VNE0Wbl+IxbBUfzaI6PEFQ+I6t0VNOxSFERER6RicTvj8adfy+J9DxOkHobpbRZymEwDDMKmwHiQpM6mVK9rxKIyIiEjH8MP7cHQnBIa5nkFzGqe2irhZDAsLty/ENM3WrGmHozAiIiL+r6oc/vP/XMsXzIaQyNPufmqriJvTdJKcl6zWkRamMCIiIv5vy5twPAM6xcJ5vzztru5WEYOG5x4xMNQ60sIURkRExL+VF8H6P7iWL/ovsIecdvdKZyXZxdmYNBw2TEyyi7OpdFa2dE07LJu3KyAiItKqkv4EJbkQ1R/OueOMu9utdt69+l02p6fz6/e+B+DFm0YyODasZp/IoEjsVnurVbmjURgRERH/dSIHkha6li/5H7AGNOmw6OAYln25H2dZD24Y05OfDhvVipUUddOIiIj/Wv8HqCyGHmNh6LVNPuy9rYf47nABnQJt/NeVg1qxggIKIyIi4q/y02DLUtdyEx+GB1BQWskLn+4BYPZlA4juHNRKFRQ3hREREfFP//kdOKug/2XQZ3KTD3tlzV7yiivoH92JuyYktF79pIbCiIiI+J/MHa5JzgAufarJh+3JLuKdrw4C8PS0YQRYdZlsCzrLIiLiX5xO+Ow3ruUR0yFuZJMOM02Tpz7+AYfT5CfDYrlgQNdWrKTUpjAiIiL+ZdN818PwbEFw8X83+bBPdmbxVVo+gTYLT04d0nr1k3oURkRExH8cTIL/POdavuoPENmnSYeVVFTxu092AfDLi/oRH3n6idGkZSmMiIiIfyjOhffvAdMBI29u0gRnbou+TCWroIyeXYKZeWG/VqykNKRZYWTRokX06dOHoKAgxo4dy4YNGxrd94MPPuDyyy+nW7duhIWFkZiYyOrVq5tdYRERkXqcTvjgF1CUBV0HwtSXm3wr78G8Yt5YnwbAb6YOJSjA2po1lQZ4HEZWrlzJ7NmzefLJJ9m+fTuTJk3iyiuvJCMjo8H9169fz+WXX86qVavYunUrF198MdOmTWP79u1nXXkREREANr4MqV+ALRhuegsCOzX50P/37xQqHE4mDejKFcNiWrGS0hjD9PCxg+eddx5jxoxh8eLFNeuGDBnCddddx7x585pUxrBhw7j55pv57W9/26T9CwsLCQ8Pp6CggLCwsDMfICIiHUf6JnjrajCdcO2rcM7Pmnzol3tyuPvNb7FZDD6dPZn+0U0PMXJmTb1+e9QyUlFRwdatW5kyZUqd9VOmTCEpKalJZTidToqKioiMjGx0n/LycgoLC+u8RERE6jnxY/U4ESeMuhVG397kQ8urHDz7rxQA7p6YoCDiRR6FkdzcXBwOBzExdZuxYmJiyM7OblIZL730EsXFxUyfPr3RfebNm0d4eHjNKz4+3pNqiohIR+B0wgc/hxPZ0HUQTH2pyeNEAJZuTOdAbjHdOgfy0KUDWrGicibNGsBqnPLLNk2z3rqGrFixgqeffpqVK1cSHR3d6H5z586loKCg5nXo0KHmVFNERPzZhpcg7UvXOJHpb4E9tMmHZheUsfA/+wB44ieD6RzUtKf5SuuwebJz165dsVqt9VpBcnJy6rWWnGrlypXce++9vPfee1x22WWn3TcwMJDAwEBPqiYiIh3JgQ2w9n9dy1NfgmjPJimb93+7KKlwMKZXBNef06MVKiie8KhlxG63M3bsWNasWVNn/Zo1a5gwYUKjx61YsYIZM2bwt7/9jalTpzavpiIiIgAncuAf97rGiYy+Hc5p+jgRgG8O5PPPHZkYBjx77XAslqZ37Ujr8KhlBGDOnDnccccdjBs3jsTERN544w0yMjKYOXMm4OpiOXLkCG+//TbgCiJ33nknCxYs4Pzzz69pVQkODiY8PLwFfxQREfF7Tkf1OJGj0G2wa5ZVDzicJk99nAzALef2YngPXYfaA4/DyM0330xeXh7PPvssWVlZDB8+nFWrVtG7d28AsrKy6sw58vrrr1NVVcUDDzzAAw88ULP+rrvuYtmyZWf/E4iISMex4SVIWwsBIa75RDwYJwLwt68PsiurkPDgAB67YlDr1FE85vE8I96geUZERIQD6+Hta13dM9e9BqNv9ejw/OIKLn5xLQWllTx77TDuTExonXpKjVaZZ0RERMQrio7C++5xIj/zOIgAvPjZHgpKKxkc25nbxvdqhUpKcymMiIhI++Z0wAf3QXEOdBvi8TgRgB+OFLDiG9cQgmeuGYbNqstfe6LfhoiItG/r/+DqogkIrZ5PJMSjw03TNWjVNOGaUd05r29UK1VUmkthRERE2q+0dbD2edfy1a9At6YNOt2cuZlrP7qWzZmb+XD7EbYePEaI3cp/X+XZfCTSNjy+m0ZERKRNFB2Ff9wHmHDOHTDq5iYdZpomC7YtIK0gjVe2LiDt+3sBePCS/sSGB7VihaW51DIiIiLtj9PhmtisOAeih3k0TiQpM4nkPNdcIrvykznm/IE+XUO594I+rVVbOUsKIyIi0v6s+z2kb3CNE7lpGQQEN+kw0zRZuH0hFsNS/dkgsNtn/M/UIQTarK1YYTkbCiMiItK+pH4J615wLU+bD90GNvlQd6uI03QCYBgm1uDDBIbtb4WKSktRGBERkfajKNs13TsmjLkLRk5v8qGntoq4WbCwcPtCfGCOzw5LYURERNoHR5VrwGrxjxAzHK78vUeHn9oq4ubESXJeMkmZSS1ZW2lBCiMiIuJ9pgmfPekaJ2Lv5NE4EdfhrlYRg4afwGtgqHWkHVMYERER73I64ZM58PVrrs/TFkDXAR4VUemsJKs4C5OGw4aJSXZxNpXOyrOtrbQCzTMiIiLe46iCf86C71cCBlyzEEbc6HExNiOAgY7/4YsDaYQEWHn+hhH07dapzj6RQZHYrfYWqri0JIURERHxjqpyeP8e2P1vsNjg+tebFUQAfv/pbtZ8X47N0pPFt5/LpAHdWriy0poURkREpO1VlMDKn0HqF2ANdD1zZtCVzSpq2aYDvL4+DYAXbhypIOKDFEZERKRtlRXCilvg4CYICIFbV0Dfi5pV1Kc/ZPPMv1MAeOyKQfx0TM8WrKi0FYURERFpOyX58NcbIHMbBIbB7e9Br/ObVdTWg/k8/O52TBNuP68Xsy7q18KVlbaiMCIiIm2j6Ci8cx3kpEBwJNzxIXQf3ayiUn88wb1vbaG8ysllQ6J55pphGEbDt/VK+6cwIiIire/4IXj7WshPhU6xcOc/IXpws4rKKSrjrqXfcLykklHxEfzx1nOwWTVThS9TGBERkdaVl+oKIgWHILwX3PVPiOzbrKKKy6u4Z9m3HD5WSkJUCEvvGkeIXZcyX6ffoIiItJ6jKa6umRNHIaq/q0UkvHmDTCsdTmYt38YPRwqJCrWz7O7xRHUKbNn6ilcojIiISOvI3A7vXA+lx1zPmrnjQ+gU3ayiTNPkyQ93sm7vjwQFWFgy41wSuoa2cIXFWxRGRESk5R3cDH+bDuWF0GMs3P4+hEQ2u7j5n+/j71sOYzHgT7eOYXR8RMvVVbxOYURERFpW6n9gxW1QVQq9L4Db3oXAzs0ubuW3GSz4Yh8A/++64Vw2NKalairthMKIiIi0nN2fwHszwFEB/S+D6e+APaTZxX25J4f//vAHAB68uD+3n9e7hSoq7YnuhRIRkZbx/Xuw8g5XEBkyDW75m0dBZHPmZq796Fo2Z252FXf4OA8s34bDafLTMT349ZSBrVVz8TKFEREROXtbl8EHPwfTASNvgRuXga3pd7qYpsmCbQtIK0hjwbYFHMwt5p5l31JS4eCC/l15/qcjNamZH1MYERGRs7P5VfjXw4AJ4+6F6xaD1bNRAEmZSSTnJQOQnJfM7cvfJvdEBUPiwlj8szHYbbpc+TP9dkVEpHnKi+Djh2D1f7s+T3gIpr4EFs8uLaZpsnD7QixG9XGmQX7gx3SPCGLZ3efSOSighSsu7Y0GsIqIiOcObIB/zoLjGa7PF/8GJj8KzehKqd0qAoBhYg0+zOwLnMSEBbVQhaU9U8uIiIg0XUUJ/N9/wVtXu4JIRC+4699w4WPNCiL1WkWqGVj4IH0Jpmm2VM2lHVMYERGRpjn0Dbx2AXz9muvz2Lvhl0nQZ1Kzi3S3ijhNZ531Jk6S85JJykw6mxqLj1A3jYiInF5VOXz5v5D0RzCd0Lk7XLMQBlx2VsU6nU6e2fgSpmlgGPVbQAwMFm5fyITuE3QnjZ9TGBERkcZl7oCPfgk5Ka7PI2+BK5+H4C5nVWxJRRX//eF2jpzIwmJruCvGxCS7OJtKZyV2q/2svk/aN4URERGpz1EJG16G9S+AswpCusK0+a7JzM7SvqNFzFq+jX05J7AGPMgdF3TlpjE9MSz1Wz8igyIVRDoAhREREakrZxd8OBOydrg+D7kGrn4FQrueddH/2HqY33z0A6WVDqI7B/LHW3/C+X2jzrpc8W0KIyIi4uJ0wOY/wX+ec03pHhThmjdk+A3NulOmttIKB099/AN/33IYgAv6d2X+LaPp2qnps7SK/1IYERERyEt1jQ059LXr84ArYNoCCIs766L355zggeXb2HO0CMOA2ZcO5MFL+mNtoFtGOiaFERGRjszphC1LYM1vobIE7J3hJ/PgnJ+ddWsIwD93HGHuBzspqXDQtVMgf7xlNBP6n313j/gXhRERkY7q+CH45wNwYJ3rc5/JcO2rronMzlJZpYNn/pXCim9cM7Qm9o1iwa2jie6sGVWlPk16JiLS0TgdsO1tWJToCiK2YLjyD3DHP5sVRDZnbubaj65lc+ZmAA7kFnP9oiRWfJOBYcBDlw7gr/edpyAijVLLiIhIR+Gogp3vwYaXIG+fa138ea6n7Eb1a1aRpmmyYNsC0grSWLBtAbk/9mLuBz9woryKqFA7828ZzaQB3VrwhxB/pDAiIuLvqsrhuxWueUOOH3StC4pwPdju/FlgsTa76NoPuUvOS2b29r/jKB/I+D6RLLz1HD3oTppEYURExF9Vlrq6YzYtgMIjrnWh3SDxQTj3XgjsfFbF1zzkDgtOnJimQWC3z5hx7k+Yc/kgbFaNBJCmURgREfE35UWwZSkk/QmKc1zrOsfBxIdhzF1gD2mRr/kyY0NNqwiAYZhYgw8zcUS+goh4RGFERMRflB6Hb96ArxZB6THXuohecMEjMPp2sLXMBGOlFQ7e2ZzOn/Y8j2mv+5A7i2HRw+3EYwojIiK+rjjPFUC+eQPKC13rovrDpF/DiJvAGtAiX1Na4WD51wd5bV0ax8ydhPQ6xKlxw2k6Sc5LJikziYk9JrbI94r/UxgREfFVRdmQtNDVJVNZ4loXPdQ1MHXodWc1MLW2skoHy7/O4LV1qfxYVA6YRPRfgwMDqP/EXQNDrSPiEYURERFfc/wQJP0Rtr4FjnLXurjRMPkxGHQVWFpmvEZZpYO/fZ3B4poQAj27BPPLi3rzRnoJ+WX1gwiAiUl2cTaVzko9cVeaRGFERMRX5O6HpAWwYwU4K13r4s+DyY9D/0s9mr59c+Zmnv/meZ4Y/wSJ3RPrbCurdLDimwwWr00lpzqE9IgI5leX9OenY3pit1m4dPi75JflN1p+ZFCkgog0mcKIiEh7VpgJP3wAP/wDMredXN/nQldLSMIFHj9D5tSJys6POx/DMCirdLDy20MsWrufo4UnQ8iDl/TnhuoQ4hYbGktsaGyL/IgiCiMiIu1NcR6kfOQKIQc3UTMuw7DCgMtdA1Pjxze7+FMnKlubsZFDmfEs+jKV7MIyALqHB/HAJf25aWx8nRAi0hoURkRE2oOyQtj9CfzwPqR+Cabj5LZeiTD8Bteg1E5nN7V6zURlhgWn6cTAwuzP5lGQOgswiAsP4oGL+3PTuJ4E2lpmAKzImSiMiIh4S2Up7F3tCiB7Pzs5GBUgbhQMvxGG/xTCe7bYV9ZuFQEwcWLaD9EtOp2HJkxjukKIeIHCiIhIW3JUulo+fnjf1RJSceLktq4DqwPIDdC1f4t+7fGSClb/kM1LKc9jUneiMgMLfQds5GfnPaBbccUrFEZERFqb0wEHk1wBJOWfJ2dHBQjv5Wr9GHEjxAz3eDDq6RSUVLI6JZtPvs9i0/5czOA9hPRKrzdRmYmTlHxNVCbeozAiItIaCjMhfROkb3B1xZzIPrktNBqGXe8KID3PbXIAOd3tuG4FJZV8lpLNqp1ZbNyfS6XD3QJiEtn9cyo1UZm0QwojIiIt4fgh150v6Rtdr2MH6m4PCoch17gCSMIkj2dHbex2XICC0ko+TznKJzuz2LDvx1oBBAbHdmbqiDguHxbFzHUvkqeJyqQdUhgREWmOYwddocMdQI4frLvdsEDsSNc8IH0mQ9+LzupBdafejvt5+nqKjvXjk++zWH9KABkU05mpI+O4akQc/aM71ax/92pNVCbtk8KIiMiZmKarpSN908kAUnCo7j6GFbqPht4TXQGk1/mu1pAW+fq6t+OCwezPnqf4wANQPQJkYEwnrhoRx9QRcQyI6dxgOZqoTNorhRERkVOZJuSlwsGNJwNIUWbdfSw26H6OK3j0vgB6nQeB9UNAU8Z5NFwFk9Qfi9mWcYzVaetILk2uvRVL0GHiexzip4MvZerIOAY2EkBEfIHCiIh0bKXH4GgK5KTA0WTXK2cXVBTV3c8SAD3GQkJ1y0fP8RDYqeEyq51unMepTpRX8d2h42w7eIxtGcfYlnGcgtJKwCQkYTmWoLq341qw0KPPOmZfNlMDTsXnKYyISMdQVQ65e6uDR7Lr/Why/RYPN6vddadL74muANJzPNhDPPrKU8d5uG+dNU2T9LwStrqDx8Fj7D1ahPOUsaWBNgv9eh3hUODhemU7cdYpU8SXKYyIiH8xTdd4jppWjhRX8MjbB86qho8J7wUxQyF6KMQMY7PVwfN7lvPEeXM96lqpW436067/z/oXSSgPYMehAvKLK+od0yMimLG9uzCmVwRjendhcGxn7vz0dow8A1O344ofUxgREd/jqHQFjuMZrtexg9Xv6fDjbigvbPi4wHBX6IgZVhM8iB5SZ6CpaZos+ORW0goPnLFrpTbTNPmxqJzUH4s5kFtMUuZGkgvrTrv+Y8V+Dh7ahKN4IHabhZE9whnjDh+9uhAdFlSnzApHBdnF2Q0GEVeZuh1X/IPCiIi0P44qKDxSHTYO1g8dRZlgOusdtjkokOe7duGJY1UkdupTq7VjuGs5rMcZJxhrrGvFrbCskgPVgSMt1/V+IPcEB34sprjC/XA7k5CEd+qN8wALfQas54XEuxjeI+KMT8O1W+26HVc6BIUREWlbjkoo/hFO5FS/H4WCI7VCx0HX59pPrW2ILQgietW8zPBeLMj+lLSSTBYMmcz5V7/rcddFY10rY22hpOeVcCC3mNwT9btX3CwGxEeGENX1APss9cd5gJOj5fspte3CbmvaOA/djisdgcKIiJy9qvKTAeNEDhTn1Aobtd9z6j6XherWjKguPJF3jMSyWk+ttdrrhA3Xq3f1qxd0iq7TypF0ZBPJaUsBSM5PqdeiYZomhWVV5BSWkVNUztHq95zCco4WlfFjYTmHyrdzIqJ+18oH+7/EUTywZn1050D6dA2lb7dQ+nQNpU/XTvTpGkqvyBACrAa3fvK6xnmIeEBhRKSDaNJ8F5VlUHYcygqgtPq9rNZ79brNxYd4viKDJ0oMEgtzXds9YVghtBtmp64sCC4njXIW9B7K+UN/idEloTpsxIDl9N0YDqdJUVklx0sqeOHr+RhYMHG1aMz98gWG81tX4KgOH+VV9bt2TjIJSfgXFrN+10qvfut4eMhN9OvWmYSuIXQOCmi0FI3zEPFcs8LIokWL+MMf/kBWVhbDhg1j/vz5TJo0qdH9161bx5w5c0hOTqZ79+48/vjjzJw5s9mVbinNnYxIZZ1dWe2xTu21rCaV43RAZQlUFFe/TkBFSa3lYsyKYhYcWEFaeS4L/vMo5wcNxigrqBs2So+Do7zh76jFBBZ0jyEtMJAFtnLOLytwzQFqCXC1VoR2q36Phk7dqt9PWR/cBSwWV2vG567/FyRX5PG2M4LepX0oyK+ksOwQhaWVFJRWUlBS/V7rVVhaSVG56+4Ya+heQnrtrlVHJ8ccaazOWF+nRQMgLMhGdFgQMWGBRHcOIrr6/bi5k2VpDXet5Fam0i36ICOacAutxnmIeM4wTbPh+N6IlStXcscdd7Bo0SImTpzI66+/zl/+8hdSUlLo1atXvf0PHDjA8OHD+fnPf87999/Ppk2bmDVrFitWrOCGG25o0ncWFhYSHh5OQUEBYWFhnlS3UaZpcusnt5Kcl8ywqGGsmLqi2U2mKsu369SqZXUZxIrJr2A4KqCqzNXyUFUGVaWuro0668rqfDYry7j12CaSHUUMM4JZYcZiVJZUB40TJ8NHVekZ67QpOIiZsdE1n1/LzmFiaVnDOxsW190lQeGYQRE4A8Nw2MOpsodRZQ9nk/MEjxesrdn9we6/okfYRArMUEoqnJRUOCiprKKk3EFJhYPSyiqKyx2UutdXOCgpd1BcUYkzbgFG4BEMw8Q0DZxlPShJPznFeRPPOJ36vIoReARqtWgYWIgO7MvsIa8SGx5UEzyCAuo/oM79O0vJS2m0a2Vo1NCz+nch0hE19frtccvIyy+/zL333st9990HwPz581m9ejWLFy9m3rx59fZ/7bXX6NWrF/PnzwdgyJAhbNmyhRdffLHJYaQ1nGnEvMpq4bJME0wnSYc31C0n/TMmRo9zDVZ0Ok55d5787Kyqty4pb2fdsr6ez8Swfq59T305GlhX65VUkkny8VplvXczE40Q12BLR0X1q/Fls9a6TQEGybHdXGUd20PS6+Mav/Cf7pwHB5FcHSCSzVKSsr8+bTlOLFRag6mwhFBhDabCCKLcEkyZJYjfR+RjmOWYhut6/XTXAVxfdBkFZijHnaEcc4aQ7wgi3xHCsSo75cVQXuDEceosXJiEJLxac5eIaRosSP2QkvTueBYgqlszgk62RBiGiTX4MH3ijxATMJLw4IB6r7AGPv+Q/w2/+rJ+i4ZZPVi0S9cDjDvDv9VKZ6W6VkS8yKOWkYqKCkJCQnjvvfe4/vrra9Y//PDD7Nixg3Xr1tU7ZvLkyZxzzjksWLCgZt2HH37I9OnTKSkpISCgft9reXk55eUnm4sLCwuJj49vsZaRU/8KMkzo7Qzhf0oGY2Bg1PwPyfVumGbN8sl1rmUTk2dCUzloLcVpgMWE3o4gnjmRUKsss2Z/oGada33t0+/kybBM0q0VNWX1cQTwvwXRWEzDtX+t42qWMcE0qy8FrmVw8mhkEWk2R01Z/SotvJIf6OoTx1mnHNdnMExn9efa9XPyy2gr+wIMnIaBxTQZWOHkL0dLsODEYjoxqD7OdGKpeXets1SXbQK3do9hl91eU86QigpWZB718DLmW2UNLq/kjcwCKgigHDvlZgBl2CkngHLTTln1+rLqz+UEUEoAK+PTyA0sofpXT+eycIZmXEYpQZSYQRQTRAlBFJuBlBBEOQE0FAhcXRhL660vybinXhfG6VgtBoGd92Ht/pd62/pUzSY2YBQhdivBdhuhdmvNckj1ckj1crDdSkiAlae2/oK0gj04OTmOw2JYGBI5pMktEC3ZopFdnH3GrhXd1SLimVZpGcnNzcXhcBATE1NnfUxMDNnZ2Q0ek52d3eD+VVVV5ObmEhcXV++YefPm8cwzz3hSNY/U/isfwDQg3VpCZfGnHv8Fuyk4iAPhJ5u/nQYcsJVRUra+WWWl2eqWlWqr5Fjl9maVtT+gbln77E4Oc5CJZZ6XtcdeuyyD3YFWfrCXelSvpOAgkgNPPkLdaRgkBwaSFBzExNIyqkwLDle8qfNeb9k0+DbYSnKgvV5ZywL7MaTUTpVpoQobVVhwYKUSKw6sVGGhyrThwEIVVqqwkhpSQnJger2yZgZcTWRJDJXYqDBtVHLyVYGNStNa9zM2SkMOURr4jzplpQTZGWP7LygdhNViYDUM13v1y2IY2NzLFrBZLFTad1EQ9MPJk2dAUXAB2X160sUygi5WC90sBgFWCzargc1iIcBqnLJswWbAv/P+Ql6la2DnyeIsDBi0kUeG3kKQ3YbdaiEwwILdaiEowILdaq357H63Wgxu/eRtduW7nxzrYjEshMSs4fWp9zS5C2PTkU3sL9hVb73T9GyK85Zs0dAttCLe06wBrKf+D8c0zdP+T6ih/Rta7zZ37lzmzJlT89ndMtIS6j+Ku7qOJsyL7s8TjvMxcI3gN0+tn+FqM6j+gInJC5aNGGaB66/XWmX9b8xAHjMvwqg5xqjex8DdxuFaZ9TU6yU+x+BYnf+xGhj8LnYIj1ivwjAsruOM6mNrLWMYNfUzTVhQ8Q8MM7deWc91H8EDwbdjsbjLstTUBcPANCyApbo8C07TZHHhnzEcmfXKerbnWH4R9QgWixXTsLrGGhjuMi2YFotrfXXZfznyNEZFer2L4jMJE7m373wshoFhuM6P4f6RMGru3nRvA5O/7HkIo3R/vbLe69+XXw//ExbDgmFAgGEQWF2OpfrHtFSXY7EYYJp8tH0mliJL3b/QsXB0aAVzEp/GanHVy2K4goRhuI61GFR/di0bwP3/uZO9x+tfrM8952tWTJ3twV/7r1HUwEW/U9znvD31Xo8u+ss+31//O3BypHQvnbqkNbm7bdORTXVCvJunAcL936DB2d/6qsGiIv7BozDStWtXrFZrvVaQnJyceq0fbrGxsQ3ub7PZiIqKavCYwMBAAmv9Bd2STm0VcTMNOGgUYkz5qUf/cz74+ScNlpXBcQIu/4lHZR36/O/1y8LkEHmEXDzRo7IOf/5ag2Uddh6lS+JAj8o68vmRBsvKrDpE7IiQJpW16cgmMg+kNVCOk6yyffTsfsizn2/H3gbLOlSyh6CwVI/KOlC0u956J05SC3eT69jJxNiml7X7WEr9sjy8WDf2b9SbF/2WLKulx2eoRUPE93kURux2O2PHjmXNmjV1xoysWbOGa6+9tsFjEhMT+de//lVn3Weffca4ceMaHC/Smtrr/5z9vaz2WKf2WlZ7vei3ZFlqzRCRU3ncTTNnzhzuuOMOxo0bR2JiIm+88QYZGRk184bMnTuXI0eO8PbbbwMwc+ZM/vSnPzFnzhx+/vOfs3nzZpYsWcKKFSta9idpgvb6P2d/L6s91qm9ltVeL/otHSDUmiEitXk8zwi4Jj174YUXyMrKYvjw4bzyyitMnjwZgBkzZpCens7atWtr9l+3bh2PPPJIzaRn//Vf/+XRpGctOc9IS46YV1lNL6s91qm9lqW7OkTEXzT1+t2sMNLWWmPSMxEREWldTb1+n/7BDyIiIiKtTGFEREREvEphRERERLxKYURERES8SmFEREREvEphRERERLxKYURERES8SmFEREREvEphRERERLzK42fTeIN7ktjCwkIv10RERESayn3dPtNk7z4RRoqKigCIj4/3ck1ERETEU0VFRYSHhze63SeeTeN0OsnMzKRz585nfGy6W2FhIfHx8Rw6dEjPs2lDOu/eofPuHTrv3qHz7h3NOe+maVJUVET37t2xWBofGeITLSMWi4WePXs269iwsDD9Y/UCnXfv0Hn3Dp1379B59w5Pz/vpWkTcNIBVREREvEphRERERLzKb8NIYGAgTz31FIGBgd6uSoei8+4dOu/eofPuHTrv3tGa590nBrCKiIiI//LblhERERHxDQojIiIi4lUKIyIiIuJVCiMiIiLiVT4dRhYtWkSfPn0ICgpi7NixbNiw4bT7r1u3jrFjxxIUFETfvn157bXX2qim/sWT8/7BBx9w+eWX061bN8LCwkhMTGT16tVtWFv/4em/d7dNmzZhs9kYPXp061bQT3l63svLy3nyySfp3bs3gYGB9OvXj6VLl7ZRbf2Hp+d9+fLljBo1ipCQEOLi4rj77rvJy8tro9r6vvXr1zNt2jS6d++OYRh89NFHZzymRa+ppo969913zYCAAPPPf/6zmZKSYj788MNmaGioefDgwQb3T0tLM0NCQsyHH37YTElJMf/85z+bAQEB5vvvv9/GNfdtnp73hx9+2Pz9739vfvPNN+bevXvNuXPnmgEBAea2bdvauOa+zdPz7nb8+HGzb9++5pQpU8xRo0a1TWX9SHPO+zXXXGOed9555po1a8wDBw6YX3/9tblp06Y2rLXv8/S8b9iwwbRYLOaCBQvMtLQ0c8OGDeawYcPM6667ro1r7rtWrVplPvnkk+Y//vEPEzA//PDD0+7f0tdUnw0j48ePN2fOnFln3eDBg80nnniiwf0ff/xxc/DgwXXW3X///eb555/fanX0R56e94YMHTrUfOaZZ1q6an6tuef95ptvNn/zm9+YTz31lMJIM3h63v/v//7PDA8PN/Py8tqien7L0/P+hz/8wezbt2+ddX/84x/Nnj17tlod/VlTwkhLX1N9spumoqKCrVu3MmXKlDrrp0yZQlJSUoPHbN68ud7+V1xxBVu2bKGysrLV6upPmnPeT+V0OikqKiIyMrI1quiXmnve33zzTVJTU3nqqadau4p+qTnn/eOPP2bcuHG88MIL9OjRg4EDB/Loo49SWlraFlX2C8057xMmTODw4cOsWrUK0zQ5evQo77//PlOnTm2LKndILX1N9YkH5Z0qNzcXh8NBTExMnfUxMTFkZ2c3eEx2dnaD+1dVVZGbm0tcXFyr1ddfNOe8n+qll16iuLiY6dOnt0YV/VJzzvu+fft44okn2LBhAzabT/5n7nXNOe9paWls3LiRoKAgPvzwQ3Jzc5k1axb5+fkaN9JEzTnvEyZMYPny5dx8882UlZVRVVXFNddcw8KFC9uiyh1SS19TfbJlxM0wjDqfTdOst+5M+ze0Xk7P0/PutmLFCp5++mlWrlxJdHR0a1XPbzX1vDscDm677TaeeeYZBg4c2FbV81ue/Ht3Op0YhsHy5csZP348V111FS+//DLLli1T64iHPDnvKSkpPPTQQ/z2t79l69atfPrppxw4cICZM2e2RVU7rJa8pvrkn0xdu3bFarXWS8k5OTn1kppbbGxsg/vbbDaioqJara7+pDnn3W3lypXce++9vPfee1x22WWtWU2/4+l5LyoqYsuWLWzfvp0HH3wQcF0kTdPEZrPx2Wefcckll7RJ3X1Zc/69x8XF0aNHjzqPTB8yZAimaXL48GEGDBjQqnX2B8057/PmzWPixIk89thjAIwcOZLQ0FAmTZrEc889p5bvVtDS11SfbBmx2+2MHTuWNWvW1Fm/Zs0aJkyY0OAxiYmJ9fb/7LPPGDduHAEBAa1WV3/SnPMOrhaRGTNm8Le//U19uM3g6XkPCwtj586d7Nixo+Y1c+ZMBg0axI4dOzjvvPPaquo+rTn/3idOnEhmZiYnTpyoWbd3714sFgs9e/Zs1fr6i+ac95KSEiyWupczq9UKnPxrXVpWi19TmzXstR1w3/q1ZMkSMyUlxZw9e7YZGhpqpqenm6Zpmk888YR5xx131Ozvvg3pkUceMVNSUswlS5bo1t5m8PS8/+1vfzNtNpv56quvmllZWTWv48ePe+tH8EmenvdT6W6a5vH0vBcVFZk9e/Y0b7zxRjM5Odlct26dOWDAAPO+++7z1o/gkzw972+++aZps9nMRYsWmampqebGjRvNcePGmePHj/fWj+BzioqKzO3bt5vbt283AfPll182t2/fXnM7dWtfU302jJimab766qtm7969Tbvdbo4ZM8Zct25dzba77rrLvPDCC+vsv3btWvOcc84x7Xa7mZCQYC5evLiNa+wfPDnvF154oQnUe911111tX3Ef5+m/99oURprP0/O+a9cu87LLLjODg4PNnj17mnPmzDFLSkrauNa+z9Pz/sc//tEcOnSoGRwcbMbFxZm33367efjw4Taute/68ssvT/v/6ta+phqmqTYsERER8R6fHDMiIiIi/kNhRERERLxKYURERES8SmFEREREvEphRERERLxKYURERES8SmFEREREvEphRERERLxKYURERES8SmFEREREvEphRERERLxKYURERES86v8DtzfAV74qPqkAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(x, c_central.value, label=\"Central\")\n", "plt.plot(x, c_upwind.value, label=\"Upwind\")\n", "plt.plot(x, c_analytic, '^', label = \"Analytic\")\n", "plt.legend();" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.9" } }, "nbformat": 4, "nbformat_minor": 5 }