ConvectionDiffusionFV.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import fipy\n",
    "from fipy import *\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from diffsolver import diffsolver\n",
    "######################################################################\n",
    "# Plot configurations\n",
    "# mpl.style.use('classic')\n",
    "mpl.rcParams['lines.linewidth'] = 1.0\n",
    "mpl.rcParams['lines.dashed_pattern'] = [6, 6]\n",
    "mpl.rcParams['lines.dashdot_pattern'] = [3, 5, 1, 5]\n",
    "mpl.rcParams['lines.dotted_pattern'] = [1, 3]\n",
    "mpl.rcParams['lines.scale_dashes'] = False\n",
    "\n",
    "mpl.rcParams['grid.color'] = 'k'\n",
    "mpl.rcParams['grid.linestyle'] = ':'\n",
    "mpl.rcParams['grid.linewidth'] = 0.5\n",
    "\n",
    "mpl.rcParams['axes.autolimit_mode'] = 'round_numbers'\n",
    "mpl.rcParams['axes.xmargin'] = 0\n",
    "mpl.rcParams['axes.ymargin'] = 0\n",
    "######################################################################\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "L = 1.\n",
    "nx = 1000\n",
    "dx = L / nx\n",
    "\n",
    "peclet = 142.\n",
    "D = 1/peclet \n",
    "convCoeff = (1,)\n",
    "\n",
    "tend = 3.5\n",
    "alfa1   = 10\n",
    "dt      = alfa1*dx\n",
    "steps = int(tend/dt)\n",
    "\n",
    "valueLeft = 1.\n",
    "valueRight = 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# initial and boundary conditions\n",
    "mesh = Grid1D(nx=nx, dx=dx)\n",
    "phi = CellVariable(name=\"solution variable\", \n",
    "                   mesh=mesh,\n",
    "                   value=0.)\n",
    "\n",
    "phi.faceGrad.constrain(valueRight, where=mesh.facesRight)\n",
    "# phi.constrain(valueRight, where=mesh.facesRight)\n",
    "phi.constrain(valueLeft, where=mesh.facesLeft)\n",
    "\n",
    "x = mesh.cellCenters[1]\n",
    "\n",
    "cout = 0 * np.ndarray(steps+1)\n",
    "pv = np.arange(0,tend+dt,dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# differential equation\n",
    "eqI = (TransientTerm() - DiffusionTerm(coeff=D) + \n",
    "       ExponentialConvectionTerm(coeff=convCoeff) == 0)\n",
    "\n",
    "# solving the equation\n",
    "for step in range(steps):\n",
    "        eqI.solve(var=phi, dt=dt)\n",
    "        cout[step+1] = phi._value[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtgAAAF3CAYAAACfcTyjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmYXGWZ/vHvQxaykIUQSCAJ6xAg\nIBhANkUjyqqCu4CioIjLgDrqzA9mGGRQB8VRUWQGFRFERdRxQYZVpNkUSdjCGghhSUIWIEl3tg5Z\nnt8fdVqLpruruqvO8tS5P9dVV/qcOlV1v6dOvf3m9FPvMXdHRERERESaY7O8A4iIiIiItBINsEVE\nREREmkgDbBERERGRJtIAW0RERESkiTTAFhERERFpIg2wRURERESaKLUBtpldZmZLzezhXu7/oJnN\nTm5/NrN90soiIiK11dFvm5l918zmJn33vllnFBGJIM0z2JcDR/Vx/9PAm9x9b+DLwA9SzCIiIrVd\nTt/99tHArsntNOB/MsgkIhJOagNsd78dWNbH/X929+XJ4t3A5LSyiIhIbbX6beA44CdecTcw1sy2\nzSadiEgcRanB/hhwfd4hRESkT5OA+VXLC5J1IiJSZXDeAczszVQG2G/oY5vTqPw5kmHDhu23/fbb\nM2rUKABWrlwJoGUta1nLIZafeOKJF919a2KyHtb5qzZSn13a5faOlazd4GzabAhr129k7csbMYPN\nBw9iyKDN8I0b2Mxg+LDNGbSZsW5dJwaMGDEcw+hcuxYzGDliBBisWb0GgJEjRwCwes2aZPuRGLB6\nzerk/pGV+1dXlrfottz9fi1rud7leXMH1mfnOsA2s72BS4Gj3f2l3rZz9x+Q1GjvtttuPmfOnIwS\nFkNbWxszZszIO0am1OZyKGObzezZvDM0YAEwpWp5MvB8940i99mRjsmiZO1cv5E/PPg8v3tgIQ/O\nb+eNO43jgJ3GsffksUzbbjRjhg8BipO3HpGyQqy8kbLCwPvs3AbYZrY98BvgJHd/Iq8cIiJSt2uA\n083sF8CBQLu7L8o5k+Rk3YaNXPmXZ7nktqfYc7sxnHTQDvzww1szYmjufxwXyZ25v+qve815YrOr\ngBnAeGAJ8CVgCIC7X2JmlwLvAbr+Z7DB3fev9bz777+/z5o1K5XMIiJpM7N76+nr8lBHv23A96jM\nNLIGOMXd++yQ1We3pnufXc4XfvkAO40fyZlH78FuE0flHUkkFQPts1P7b6a7n1Dj/lOBU9N6fRER\n6Z86+m0H/jGjOFJA7s7/3PYUl935NF95514ctZcmkRHpSVFmEalb1xcqyqStrS3vCJlTm8uhjG0u\nm2h9dqRjMuusGzc5//77h7n2wUVce8ah/R5ca9+mJ1LeSFkboUIpERER6ZO7c9ZvZjN/2Vqu/sRB\njBo2JO9IIoWWWg12WlTPJyKRFbkGOw3qs1vDt29+glvnLOUXpx2kLzFKqRSuBltERETiu+Hhxfzv\nfQv47adfr8G1SJ1Ugx1AWeqVqqnN5VDGNpdNtD470jGZRdYlHZ2c/buH+e4J09l61OYNPZf2bXoi\n5Y2UtRHhBtgiIiKSPnfni796kA8dtD37br9l3nFEQlENtohIhlSDLVFc/9AivnPLk1x7xhsYPEjn\n46ScBtpn6xMjIiIir9C5fiP/ef1jnPP2aRpciwxAuE9NtHq+ZihLvVI1tbkcytjmsonWZ0c6JtPM\netldT7PHxNEc8g/jm/ac2rfpiZQ3UtZGhBtgi4iISHrWvLyBH93xNP9y1O55RxEJSzXYIiIZUg22\nFN1ldz7NPU8v45KT9ss7ikjuVIMtIiIiDXl5wyYuvWMen5qxS95RREILN8COVs/XDGWpV6qmNpdD\nGdtcNtH67EjHZBpZ/++h59lx/Ej2mTK26c9d9n2bpkh5I2VtRLgBtoiIiKTjqnvm8+GDd8g7hkh4\nqsEWEcmQarClqOa9sIr3f/9u/nzmYQwdrPNvIqAabBEREWnA1bPm8559J2lwLdIE4T5F0er5mqEs\n9UrV1OZyKGObyyZanx3pmGxm1vUbN/G/9y7k/a+b0rTn7K6s+zYLkfJGytqIcANsERERaa4/P/US\nk7cczi5bb5F3FJGWoBpsEZEMqQZbiujM/53NzluP5LQ3ano+kWqqwRYREZF+27BxEzc9uoSj99o2\n7ygiLSPcADtaPV8zlKVeqZraXA5lbHPZROuzIx2Tzcp6z9PLmDR2OFPGjWjK8/WmjPs2K5HyRsra\niHADbBEREWme6x9ezFF7Tcw7hkhLUQ22iEiGVIMtReLuHHz+n/jZxw/UFxxFeqAabBEREemXJ5as\nYshgY+fxI/OOItJSwg2wo9XzNUNZ6pWqqc3lUMY2l020PjvSMdmMrLc/8QKH7ro1ZtZ4oBrKtm+z\nFClvpKyNCDfAFhERkea4/ckXeOOuW+cdQ6TlqAZbRCRDqsGWouhcv5H9vnwzf/nXtzB62JC844gU\nkmqwRUREpG73PL2MPbYdrcG1SArCDbCj1fM1Q1nqlaqpzeVQxjaXTbQ+O9Ix2WjWO558gTfsOr45\nYepQpn2btUh5I2VtRLgBtoiIiDTunqeXcfDOW+UdQ6QlqQZbRCRDqsGWIljz8gb2+/Ifuf+cwxk2\nZFDecUQKSzXYIiIiUpf7n1vBtO1Ga3AtkpJwA+xo9XzNUJZ6pWpqczmUsc1lE63PjnRMNpJ15jPL\n2H/HLZsXpg5l2bd5iJQ3UtZGhBtgi4iISGNmPbOcA3Ycl3cMkZalGmwRkQypBlvytmHjJl573s3c\n+f/ezNgRQ/OOI1JoqsEWERGRmh5d1MGkscM1uBZJUbgBdrR6vmYoS71SNbW5HMrY5rKJ1mdHOiYH\nmvX+51aw7w5jmxumDmXYt3mJlDdS1kaEG2CLiIjIwM1e0M7ek7MfYIuUiWqwRUQypBpsydvh37qN\nb3/gtew1aUzeUUQKTzXYIiIi0qfV6zYwf/kapk4YlXcUkZYWboAdrZ6vGcpSr1RNbS6HMra5bKL1\n2ZGOyYFkfeT5DnabOJqhg7P/9d/q+zZPkfJGytqIcANsERERGZjZC1awt0pDRFKnGmwRkQypBlvy\n9Jmr7ufQXcfzvv2n5B1FJATVYIuIiEifHlqoGUREshBugB2tnq8ZylKvVE1tLocytrlsovXZkY7J\n/mbt6FzPko5Odtl6ZDqBamjlfZu3SHkjZW1EuAG2iIiI9N+cxSvZdcIoBg/Sr36RtKVWg21mlwFv\nB5a6+1493G/Ad4BjgDXAye5+X63nVT2fiESmGmzJy5V3P8sjC9v52nv2zjuKSBhFrMG+HDiqj/uP\nBnZNbqcB/5NiFhERkVJ7fFEHu0/U/NciWUhtgO3utwPL+tjkOOAnXnE3MNbMtq31vNHq+ZqhLPVK\n1dTmcihjm8smWp8d6Zjsb9bHF69kt4mj0wlTh1bet3mLlDdS1kYMzvG1JwHzq5YXJOsWdd/QzE6j\ncpaboRP/gR3P/L9MAhbKDWpzKajNIpICd+eJxSt1BlskI3kOsK2HdT0WhLv7D4AfQFLP97W3pZlL\nRCQ19vW8E2Rr1KhYA7oZM2bkHaFu/cm6YPlaRm4+mC1HDk0vUA2tum+LIFLeSFkbkedXiRcA1TPd\nTwaezymLiIhIy3p88Up23zbWf3ZEIstzgH0N8GGrOAhod/dXlYd0F62erxnKUq9UTW0uhzK2uWyi\n9dmRjsn+ZH18UQe75Vwe0qr7tggi5Y2UtRGplYiY2VXADGC8mS0AvgQMAXD3S4DrqEzRN5fKNH2n\npJVFRESkzB5fspLD95iQdwyR0khtHuy0aE5VEYlM82BLHt76rdv47vHTmbZdfrOIiERUxHmwRURE\nJGfrN25i/rI17JzTJdJFyijcADtaPV8zlKVeqZraXA5lbHORmdlRZjbHzOaa2Zk93L+9md1qZveb\n2WwzO6bWc0brsyMdk/VmfW7ZGiaOGcawIYPSDVRDK+7booiUN1LWRoQbYIuISPOZ2SDgYipX2Z0G\nnGBm07ptdjbwS3efDhwP/He2KWUg5i5dxS5bb5F3DJFSUQ22iEiGilqDbWYHA+e6+5HJ8lkA7n5+\n1TbfB+a5+9eT7b/p7of09bzqs/P3321zWbFmPf96zB55RxEJRzXYIiLSiN6urlvtXOBDycxQ1wFn\nZBNNGlE5g636a5EshRtgR6vna4ay1CtVU5vLoYxtLrB6rq57AnC5u0+mMs3qlWb2qt8jZnaamc0y\ns1nz5s17xfvc1tZW6OULL7ywUHn6Wu76udb29z+1iH/YZosweYuw3D1z3nlaKe+FF15YqDy1lgfM\n3UPdpk6d6mVz66235h0hc2pzOZSxzcAsL0Bf2v0GHAzcWLV8FnBWt20eAaZULc8DtunreaP12ZGO\nyXqybtq0yfc85wZfvnpd+oFqaLV9WySR8kbK6j7wPls12CIiGSpwDfZg4AngLcBCYCZwors/UrXN\n9cDV7n65me0B3AJM8j5+kajPztfi9k7eftEdzDr78LyjiISkGmwRERkwd98AnA7cCDxGZbaQR8zs\nPDM7NtnsC8DHzexB4Crg5L4G15K/p17QDCIieQg3wFYNdjmozeVQxjYXmbtf5+5T3X0Xd/9qsu4c\nd78m+flRd3+9u+/j7q9195tqPWe0PjvSMVlP1rlLV7HLNsUYYLfavi2SSHkjZW1EuAG2iIiI1Gee\nzmCL5EI12CIiGSpqDXZa1Gfn68OX3cPJh+zAYbtPyDuKSEiqwRYREZFXeObF1ey4lebAFslauAF2\ntHq+ZihLvVI1tbkcytjmsonWZ0c6JmtlfXnDJha3dzJ5yxHZBKqhlfZt0UTKGylrI8INsEVERKS2\nBcvXMHHMMIYO1q96kaypBltEJEOqwZas3Pr4Ui6762mu/NiBeUcRCUs12CIiIvI3T7+4mp3Gq/5a\nJA/hBtjR6vmaoSz1StXU5nIoY5vLJlqfHemYrJX12ZdWs0OBvuDYSvu2aCLljZS1EeEG2CIiIlLb\n0y+tYafxxfiCo0jZqAZbRCRDqsGWrLzpG7dy2cmv04VmRBqgGmwREREBYP3GTSxa0cmUgkzRJ1I2\n4QbY0er5mqEs9UrV1OZyKGObyyZanx3pmOwr64Lla5kwZvNCTdHXKvu2iCLljZS1EcX55ImIiEhT\nPPOSruAokifVYIuIZEg12JKFK//yDI8uWsn5735N3lFEQlMNtoiIiAAwf/lapowbnncMkdIKN8CO\nVs/XDGWpV6qmNpdDGdtcNtH67EjHZN812GuYXLAvOLbKvi2iSHkjZW1EuAG2iIiI9G3+srVM2VJn\nsEXyohpsEZEMqQZbsjD9vJu46Z/exNajNs87ikhoqsEWERERVq3bQOf6TYzfYmjeUURKK9wAO1o9\nXzOUpV6pmtpcDmVsc9lE67MjHZO9ZZ2/bA2TtxyOmWUbqIZW2LdFFSlvpKyNCDfAFhERkd4tWL6W\nyaq/FsmVarBFRDKkGmxJ22V3Ps0zL63mvOP2yjuKSHiqwRYRERGdwRYpgHAD7Gj1fM1Qlnqlampz\nOZSxzWUTrc+OdEz2WoO9fA1TCjYHNrTGvi2qSHkjZW1EuAG2iIiI9K5yBrt4A2yRMlENtohIhlSD\nLWlyd/Y+9ybu+H9vZuwITdMn0ijVYIuIiJRc+9r1ODBm+JC8o4iUWrgBdrR6vmYoS71SNbW5HMrY\n5rKJ1mdHOiZ7ytr1BceizYEN8fdtkUXKGylrI8INsEVERKRnlYvMqP5aJG+qwRYRyZBqsCVNP7x9\nHs+3r+VL79gz7ygiLUE12CIiIiU3f7nOYIsUQbgBdrR6vmYoS71SNbW5HMrY5rKJ1mdHOiZ7q8Ge\nUtCLzETft0UWKW+krI0IN8AWERGRns1ftoYp43QGWyRvqsEWEcmQarAlLe7Onl+6kb+c9RZN0yfS\nJKrBFhERKbGOtRvYzEyDa5ECSHWAbWZHmdkcM5trZmf2cP/2Znarmd1vZrPN7Jhazxmtnq8ZylKv\nVE1tLocytrlsovXZkY7J7lmfb1/LtmOG5ROmDpH3bdFFyhspayMGp/XEZjYIuBg4HFgAzDSza9z9\n0arNzgZ+6e7/Y2bTgOuAHdPKJCLSysxs3zo2W+/uD6UeRjK3qH0t244t5hccRcomtRpsMzsYONfd\nj0yWzwJw9/Ortvk+MM/dv55s/013P6Sv51U9n4hElmYNtpmtBGYCfV3Gbyd33zGN1++J+uzsXHn3\nszz6fDvnv3vvvKOItIyB9tmpncEGJgHzq5YXAAd22+Zc4CYzOwMYCby1pycys9OA0wAmTJhAW1sb\nM2bMAP7+pwYta1nLWo6ynKKZ7n5YXxuY2Z/SDiH5WLRiLduO0RlskSJI8wz2+4Aj3f3UZPkk4AB3\nP6Nqm88nGb6ZnMH+EbCXu2/q7Xl32203nzNnTiqZi6qt6j8UZaE2l0MZ21y2WUSi9dmRjsnuWf/p\n6gc4ZJeteN/+U/IL1YfI+7boIuWNlBWKOYvIAqD6Uz4ZeL7bNh8Dfgng7n8BhgHjU8wkIlIqZnZu\n3hkkGwtXrGU71WCLFEKaZ7AHA08AbwEWUqkLPNHdH6na5nrgane/3Mz2AG4BJnkfoVTPJyKRZX0G\n28zuc/d6vvyYCvXZ2Tn0gj/xk48eyE7jR+YdRaRlFO4MtrtvAE4HbgQeozJbyCNmdp6ZHZts9gXg\n42b2IHAVcHJfg2sREem3vr7wKC1i0yZnSfu6Qk/TJ1Imqc6D7e7XuftUd9/F3b+arDvH3a9Jfn7U\n3V/v7vu4+2vd/aZazxltTtVmyOCLUYWjNpdDGducg/3yfPFofXakY7I664ur17HFsMEMGzIov0A1\nRN23EUTKGylrI3QlRxGRFpHMuPQK3b803tM2Et/zKzrZbqzOXosURWo12GlRPZ+IRJbyPNjzgC/2\ntQlwnrvvmcbr90R9djauf2gRv7l/IT/8cGkmqBHJRBHnwRYRkWzdBryjxjY3ZxFEsvV8eyfbqf5a\npDDClYhEq+drhrLUK1VTm8uhjG1Ok7ufUsftc1lmitZnRzomq7MuWlH8y6RH3bcRRMobKWsjwg2w\nRURE5JWeb9cc2CJFohpsEZEMle1Kjuqzs/HOi+/i7Lftwf47jss7ikhLKdw82CIiIpKNRe3FLxER\nKZNwA+xo9XzNUJZ6pWpqczmUsc1ZMLPNzexEM/tXMzun65ZHlmh9dqRjsivr+o2bWLb6ZSaM2jzf\nQDVE3LdRRMobKWsjNIuIiEjr+T3QDtwLrMs5i6RscXsn47fYnMGDwp0zE2lZfdZgm1lHrccDi9x9\nalNT9UH1fCISWRY12Gb2sLvvleZr1Et9dvrueXoZX7v+MX7z6dfnHUWk5aRVg/2Uu4/u4zYKWD2w\nyCIikpI/m9lr8g4h2VD9tUjx1Bpgv6eO56hnm6aJVs/XDGWpV6qmNpdDGduckTcA95rZHDObbWYP\nmdnsPIJE67MjHZNdWZ9fEeMiMxH3bRSR8kbK2og+a7DdfV73dWb2dne/tq9tREQkV0fnHUCys7h9\nLdtvNTLvGCJSpd/zYJvZfe6+b0p5alI9n4hEltU82Ga2D3BosniHuz9Yx2OOAr4DDAIudfev9bDN\n+4FzAQcedPcT+3pO9dnp+8SVszh2n0m8be9t844i0nKynAfbBvAYERHJiJl9FvgZsE1y+6mZnVHj\nMYOAi6mc/Z4GnGBm07ptsytwFvB6d98TyPSy69KzxR3rmDim2FP0iZTNQAbYn2h6in6IVs/XDGWp\nV6qmNpdDGduckY8BB7r7Oe5+DnAQ8PEajzkAmOvu89z9ZeAXwHHdtvk4cLG7Lwdw96W1gkTrsyMd\nk11Zl7R3MmG0arCbKVJWiJU3UtZG9DnANrNXlYK4+z21thERkVwZsLFqeSO1//o4CZhftbwgWVdt\nKjDVzO4ys7uTkpJXv7jZaWY2y8xmvfDCC6/4hdrW1lbo5QceeKBQeWot/+nWW3lhZSfbjBpWiDxa\n1nKt5QceeKBQeWotD1StebAfBGbQd8d8i7tPH3CCflI9n4hEltE82J8HPgL8Nln1TuByd7+wj8e8\nDzjS3U9Nlk8CDnD3M6q2uRZYD7wfmAzcAezl7it6e1712ela0tHJ2757B7POPjzvKCItaaB9dq0r\nOY6hciWwvgbYL/T3RUVEJD3u/i0za6MyXZ8Bp7j7/TUetgCYUrU8GXi+h23udvf1wNNmNgfYFZjZ\nlODSb4uDlIeIlE2fJSLuvqO77+zuO/VxOyCrsBCvnq8ZGvkTRVRqczmUsc1pMrPRyb/jgGeAnwJX\nAs8m6/oyE9jVzHYys6HA8cA13bb5HfDm5DXGUykZ6XOq1mh9dqRjsq2tjcUdnUwMMsCOtm8jiZQ3\nUtZG1DqDLSIicfwceDuVvzxW1/9Zsrxzbw909w1mdjpwI5Vp+i5z90fM7Dxglrtfk9x3hJk9SqWu\n+5/d/aV0miL1WNLRyYQAF5kRKZt+z4OdN9XziUhkWc2DXRTqs9N1wQ2PM2zIID7zll3zjiLSkrKc\nB1tERArMzG6pZ53EF6lERKRM6h5gm9kkMzvEzN7YdUszWG+i1fM1Q1nqlaqpzeVQxjanycyGJbXW\n481sSzMbl9x2BLbLI1O0PjvSMdnW1haqRCTavo0kUt5IWRtRVw22mX0d+ADQVXcHlXq+21PKJSIi\n/fcJKldX3I5XzgDVQeUqjdJiFrfrDLZIEdVVg51MxbS3u69LP1LfVM8nIpFlNA/2Ge5+UZqvUS/1\n2ena60s3ctf/O4wxI4bkHUWkJaU1D3aXecAQIPcBtoiI9M3dLzKzvYBpwLCq9T/JL5U026p1G9iw\naROjh2tCMJGiqbcGew3wgJl938y+23VLM1hvotXzNUNZ6pWqqc3lUMY2Z8HMvgRclNzeDFwAHJtH\nlmh9dqRj8pqbb2fi6GGY9XUtuOKItG8jZYVYeSNlbUS9/+29hldfcEBERIrpvcA+wP3ufoqZTQAu\nzTmTNNmKdc6E0cPzjiEiPah7Huzkyl5Tk8U5yaVyM6d6PhGJLKMa7Hvc/QAzu5fKGeyVwMPuvmea\nr9sT9dnp+c19C7jtiRf4zvHT844i0rJSrcE2sxnAFVQuvWvAFDP7iLtrFhERkeKZZWZjgR9SmU1k\nFXBPvpGk2TQHtkhx1VuD/U3gCHd/k7u/ETgS+HZ6sXoXrZ6vGcpSr1RNbS6HMrY5bVYpyD3f3Ve4\n+yXA4cBH3P2UPPJE67MjHZP3PvoUEwINsCPt20hZIVbeSFkbUe8Ae4i7z+lacPcnqMwqIiIiBeKV\nur/fVS0/4+6zc4wkKVmxzpkY5CIzImVT7zzYl1G5sMyVyaoPAoPzOCOiej4RiSyjGuyLgcvdfWaa\nr1MP9dnpOe7iuzjn7dPYb4ct844i0rLSngf7U8A/Ap+hUoN9O/Df/X0xERHJxJuBT5jZs8BqKv22\nu/ve+caSZlrS3qkz2CIFVVeJiLuvc/dvufu73f1d7v7tvK7qGK2erxnKUq9UTW0uhzK2OSNHA7sA\nhwHvAN6e/Ju5aH12lGNy4ybnhZWdbDNq87yj1C3KvoVYWSFW3khZG9HnGWwz+6W7v9/MHqJSIvIK\nOhsiIlJIX3H3k6pXmNmVwEm9bC/BvLhqHSOHGEMG1ftVKhHJUp812Ga2rbsvMrMderrf3Z9NLVkv\nVM8nIpFlVIN9n7vvW7U8CHjI3ael+bo9UZ+djgfnr+DffvcQ155xaN5RRFraQPvsPv/r6+6Lkh8/\n7e7PVt+ATw8kqIiIpMPMzjKzlcDeZtaR3FYCS4Hf5xxPmkhzYIsUW71/Wzq8h3VHNzNIvaLV8zVD\nWeqVqqnN5VDGNqfJ3c9391HAN9x9dHIb5e5buftZeWSK1mdHOSaXdHSycdWyvGP0S5R9C7GyQqy8\nkbI2olYN9qeonKne2cyq51EdBdyVZjARERkYdz/LzCYBO1DVz+vqu61jcXsnWw6zvGOISC9q1WCP\nAbYEzgfOrLprpbvn8l9n1fOJSGQZ1WB/DTgeeBTYmKx2dz82zdftifrsdHz+lw9w0M5b8f79p+Qd\nRaSlpTIPtru3A+3ACcmLbAMMA7Ywsy3c/bmBhBURkVS9C9gtr+lUJX1LVIMtUmh11WCb2TvM7Eng\naeA24Bng+hRz9SpaPV8zlKVeqZraXA5lbHNG5gFD8g4B8frsKMfk4vZO5s95KO8Y/RJl30KsrBAr\nb6Ssjaj3So5fAQ4C/uju083szSRntUVEpHDWAA+Y2S3A385iu/tn8oskzbSkYx1jhw3NO4aI9KLP\nGuy/bWQ2y933N7MHgenuvsnM7nH3A2o87ijgO8Ag4FJ3/1oP27wfOJfKhWwedPcT+3pO1fOJSGQZ\n1WB/pKf17n5Fmq/bE/XZzbdq3QZe95U/8uh5R2KmLzqKpCmVGuwqK8xsC+B24GdmthTYUCPQIOBi\nKlP8LQBmmtk17v5o1Ta7AmcBr3f35UmNt4iINMDdrzCz4cD27j4n7zzSXEs7OpkwenMNrkUKrN55\nsI+j8ifHfwJuAJ4C3lHjMQcAc919nru/DPwieZ5qHwcudvflAO6+tFaQaPV8zVCWeqVqanM5lLHN\nWTCzdwAPUOmvMbPXmtk1eWSJ1mdHOCYXd3QyYfSwEFmrRcobKSvEyhspayNqDrCTM9G/d/dN7r7B\n3a9w9++6+0s1HjoJmF+1vCBZV20qMNXM7jKzu5OSEhERacy5VE5yrABw9weAnfIMJM2ztGMdEzSD\niEih1Rxgu/tGYE0yJ3Z/9PS3q+4F34OBXYEZVL40eamZjX3VE5mdZmazzGxWe3v7K/7309bW1vLL\n1YqQJ4vlGTNmFCpPFsszZswoVJ4slrvWFSVPFssZ2ZBMs1qt9hduUjBq1Kg8XnbAuvqeIluSlIhE\nyFotUt5IWSFW3khZG1Hvlxx/SWUWkZuB1V3r+/pGupkdDJzr7kcmy2cljzm/aptLgLvd/fJk+Rbg\nTHef2dvz6gszIhJZRl9y/BFwC5ULhL0H+AwwxN0/mebr9kR9dvOd94dH2W7sME49dOe8o4i0vIH2\n2fXWYP8f8O9UvuR4b3Kr1WPOBHY1s53MbCiVq4p1rwH8HfBmADMbT6VkZF5fTxqtnq8ZMjzrVRhq\nczmUsc0ZOQPYk8oUfT+ncsHWQi9FAAAgAElEQVSwz+URJFqfHeGYXKIa7NRFygqx8kbK2oh6ZxEZ\n6+7fqV5hZp/t6wHuvsHMTgdupDJN32Xu/oiZnQfMcvdrkvuOMLOuy/n+cx213SIi0gd3XwP8W3KT\nFtM1wF6zPO8kItKbektE7nP3fbutu9/dp6eWrBf6c6OIRJZRicjNwPvcfUWyvCXwi66SvSypz26+\nQy/4Ez/92IHssNXIvKOItLxU5sE2sxOAE4Gduk3xNArQmWYRkWIa3zW4BtB1BlqHu7NEs4iIFF6t\nGuw/A98EHk/+7bp9AchlSr1o9XzNUJZ6pWpqczmUsc0Z2WRm23ctmNkO5DSLSLQ+u+jH5Io16xk+\nZBDDhgwqfNbuIuWNlBVi5Y2UtRF9nsF292eBZ4GDs4kjIiJN8G/AnWZ2W7L8RuC0HPNIkyxZWZmi\nT0SKrd4a7HcDXwe2oTK/tQHu7qPTjfdqqucTkciyqMFOXmc8lelVDfiLu7+Y9mv2RH12c932xAtc\nesc8rvzYgXlHESmFVGqwq1wAvMPdH+vvC4iISC42B5ZR6eenmRnufnvOmaRBSzo62WaU6q9Fiq7e\nebCXFGVwHa2erxnKUq9UTW0uhzK2OQtm9nXgLiqlIv+c3L6YR5ZofXbRj8kl7Z1MHFMpESl61u4i\n5Y2UFWLljZS1EfWewZ5lZldTuTDMuq6V7v6bVFKJiEgj3gns5u7ram4poSxZ2cnUCbEuPy9SRvXW\nYP+4h9Xu7h9tfqS+qZ5PRCLLaB7s66nMg70qzdeph/rs5vr4T2bxnn0nc9ReE/OOIlIKqdZgu/sp\n/Y8kIiI5WQM8YGa38Mq/On4mv0jSDJWrOGoWEZGiq6sG28ymmtktZvZwsry3mZ2dbrSeRavna4ay\n1CtVU5vLoYxtzsg1wJepXMvg3qpb5qL12UU/Jrsukw7Fz9pdpLyRskKsvJGyNqLeGuwfUvmSzPcB\n3H22mf0c+EpawUREZGDc/QozGwpMTVbNcff1eWaSxm3c5Cxb/TJbj9IZbJGiq7cGe6a7v87M7nf3\n6cm6B9z9takn7Eb1fCISWUY12DOAK4BnqMyDPQX4SB7T9KnPbp6lHZ0c8907mXX2W/OOIlIaac+D\n/aKZ7UJyqV0zey+wqL8vJiIimfgmcIS7z4FKmR9wFbBfrqmkIYtVfy0SRr3zYP8jlfKQ3c1sIfA5\n4FOppepDtHq+ZihLvVI1tbkcytjmjAzpGlwDuPsTwJA8gkTrs4t8TC7pWPe3+msodtaeRMobKSvE\nyhspayPqnUVkHvBWMxsJbObusXpMEZFymWVmPwKuTJY/SE5fcpTmqf6Co4gUW7012P8JXODuK5Ll\nLYEvuHvmM4monk9EIsuoBntzKn95fAOVGuzbgf/O48Iz6rOb51s3zWGzzYzPvXVq7Y1FpCkG2mfX\nWyJydNfgGsDdlwPH9PfFREQkE4OB77j7u939XcB3gUE5Z5IGLdYZbJEw6h1gD0rOiABgZsOBXL5p\nEa2erxnKUq9UTW0uhzK2OSO3AMOrlocDf8wjSLQ+u8jHZKUG+++/eouctSeR8kbKCrHyRsraiHpn\nEfkpcEtyyXQHPkplCigRESmeYdWXSXf3VWY2Is9A0jjVYIvEUVcNNoCZHQ28hUo9303ufmOawXqj\nej4RiSyjGuy7gDPc/b5keT/ge+5+cJqv2xP12c2z75dv5qZ/eiPjt9BUfSJZSXsebNz9euD6/r6A\niIhk7nPAr8zs+WR5W+ADOeaRBq3bsJGVnesZN2Jo3lFEpA511WCb2bvN7EkzazezDjNbaWYdaYfr\nSbR6vmYoS71SNbW5HMrY5iy4+0xgdyrXK/g0sIe75zJNX7Q+u6jH5NKOdWy9xeZstpn9bV1Rs/Ym\nUt5IWSFW3khZG1HvGewLgHe4+2NphhERkeZw9/XAw3nnkOZYurKTCWNUfy0SRb3zYN/l7q/PIE9N\nqucTkciyqMEuEvXZzfF/sxfxhwef55KTdLV7kSylPQ/2LDO72sxOSMpF3m1m7+7vi4mISHGZ2VFm\nNsfM5prZmX1s914zczMrzX8U8laZQURfbhSJot4B9mhgDXAE8I7k9va0QvUlWj1fM5SlXqma2lwO\nZWxzFszslnrWdbt/EHAxcDQwDTjBzKb1sN0o4DPAX+vJEq3PLuoxuWRlJ9t0m6KvqFl7EylvpKwQ\nK2+krI2oqwbb3U9JO4iIiDTGzIYBI4DxZrYllWlVoXKSZLsaDz8AmOvu85Ln+gVwHPBot+2+TOV7\nOV9sVm6pbWnHOqZuMyrvGCJSp7oG2GY2GbgIeD2VC83cCXzW3RekmK1Ho0aVr4OZMWNG3hEypzaX\nQxnbnLJPUJmibzvgvqr1HVTOTvdlEjC/ankBcGD1BmY2HZji7teaWa8DbDM7DTgNYMKECbS1tf3t\nve46e1XU5a51RcnTtby4fRgTRg97xf0zZswoTL56lqPl1XJ6y12KkqfevP1R75ccbwZ+DlyZrPoQ\n8EF3P3zArzxA+sKMiESW0YVmznD3i/r5mPcBR7r7qcnyScAB7n5GsrwZ8CfgZHd/xszagC+6e58d\nsvrs5jjsm218/0P7seuE8p1kEslT2l9y3Nrdf+zuG5Lb5cDW/X2xZohWz9cMjfwPKiq1uRzK2OaM\nXGZmZ5vZDwDMbFczq/W9mQXAlKrlycDzVcujgL2ANjN7BjgIuKbWFx2j9dlFPSaXdqxTDXaGImWF\nWHkjZW1EvQPsF83sQ2Y2KLl9CHgpzWAiIjJglwEvA4ckywuAr9R4zExgVzPbycyGAscD13Td6e7t\n7j7e3Xd09x2Bu4Fja53BlsatWreBjZuc0cPqvviyiOSs3hKR7YHvAQdTqcH+M5Ua7GfTjfdq+nOj\niESWUYnILHff38zud/fpyboH3X2fGo87BrgQGARc5u5fNbPzgFnufk23bdtQiUgmnnphFR+7fCZt\n//zmvKOIlM5A++x6ZxF5Dji236lERCQPL5vZcConRDCzXYB1tR7k7tcB13Vbd04v285oPKbUY0nH\nq6foE5Fiq6tExMyuMLOxVctbmtll6cXqXbR6vmYoS71SNbW5HMrY5ox8CbgBmGJmPwNuAf4ljyDR\n+uwiHpNLO9YxoYcBdhGz9iVS3khZIVbeSFkbUW9B197uvqJrwd2XJ9M1iYhIwbj7zWZ2H5UvIhqV\nkr4Xc44lA7Sko5OJuoqjSCj11mA/CMxw9+XJ8jjgNnd/Tcr5XkX1fCISWZo12Ga2b1/3u/t9fd2f\nBvXZjfuPPzzCpLHDOfXQnfOOIlI6qdZgA98E/mxmv6ZS0/d+4Kv9fTEREUnVN/u4z4HDsgoizbO0\nYx3Tt98y7xgi0g911WC7+0+A9wBLgBeAd7v7lX0/Kh3R6vmaoSz1StXU5nIoY5vT5O5vdvc3A0d3\n/Vy17pg8MkXrs4t4TC7p6GTCqFeXiBQxa18i5Y2UFWLljZS1EXVPqunujwKPpphFRESa489A93KR\nntZJAIs7Opk4RrOIiERSVw12kaieT0QiS7kGeyIwCfgpcCKVLzgCjAYucffd03jdvqjPboy7s9vZ\nN/Dgl45g+NBBeccRKZ20a7BFRKT4jgROpnKZ829VrV8J/GsegaQxK9asZ9iQzTS4Fgmm3kulF0a0\ner5mKEu9UjW1uRzK2OY0ufsVSb31yd1qsI9199/kkSlan120Y3JxR2ePc2BD8bLWEilvpKwQK2+k\nrI3QGWwRkdazl5nt2X2lu5+XRxgZONVfi8SkGmwRkQylWYNd9RpfqFocBrwdeMzdP5rm6/ZEfXZj\nrrrnOe5/bjkXvHefvKOIlJJqsEVEBAB3f8V82Gb2X8A1OcWRBixq72TimOF5xxCRflINdgBlqVeq\npjaXQxnbnJMRQC6XAYzWZxftmFzcvpZteykRKVrWWiLljZQVYuWNlLURqQ6wzewoM5tjZnPN7Mw+\ntnuvmbmZpfpnUxGRMjCzh8xsdnJ7BJgDfCfvXNJ/lTPYqsEWiSa1GmwzGwQ8ARwOLABmAickF6yp\n3m4U8H/AUOB0d++zWE/1fCISWUY12DtULW4Alrj7hjRfszfqsxtz+Ldu46ITp7P7xNF5RxEppYH2\n2WmewT4AmOvu89z9ZeAXwHE9bPdl4AKgM8UsIiKl4e7PAmOBdwDvAqblm0gGanF7J9uOVg22SDRp\nDrAnAfOrlhck6/7GzKYDU9z92r6eyMxOM7NZZjZr3rx5r6jfaWtra/nlCy+8sFB5sljuWleUPFks\nd2973nmyWL7wwgsLlSeL5SyY2WeBnwHbJLefmdkZmbx4N6rBHrhV6zawftMmRg/veT6CImWtR6S8\nkbJCrLyRsjbE3VO5Ae8DLq1aPgm4qGp5M6AN2DFZbgP2r/W8U6dO9bK59dZb846QObW5HMrYZmCW\np9Tvdt2A2cDIquWRwOy0X7enW7Q+u0jH5JNLVvqMb9za6/1FylqPSHkjZXWPlTdSVveB99lp1mAf\nDJzr7kcmy2clA/rzk+UxwFPAquQhE4FlwLHeRx226vlEJLKMarAfAl7n7p3J8jBgpru/Js3X7Yn6\n7IG788kXufjWuVx12kF5RxEprSLOgz0T2NXMdgIWAscDJ3bd6e7twPiuZTNrA77Y1+BaRETq8mPg\nr2b222T5ncCPcswjA7Cojyn6RKTYUqvB9so31k8HbgQeA37p7o+Y2XlmduxAnzdaPV8zlKZeqYra\nXA5lbHMW3P1bwClU/iq4HDjF3S/s+1HpiNZnF+mYXFxjir4iZa1HpLyRskKsvJGyNiLVKzm6+3XA\ndd3WndPLtjPSzCIiUibufh9wX945ZOAWdXSyx8RReccQkQFIrQY7LarnE5HIsqjBLhL12QP30ctn\ncsIB23P4tAl5RxEprSLOgy0iIiIDtKi9UzXYIkGFG2BHq+drhrLUK1VTm8uhjG0um2h9dpGOycXt\na1WDnZNIWSFW3khZGxFugC0iItLqOtdvZPW6jYwbMTTvKCIyAKrBFhHJkGqwpR7PvLiaky77K3f8\ny2F5RxEpNdVgi4iItIhF7Z1sO3p43jFEZIDCDbCj1fM1Q1nqlaqpzeVQxjaXTbQ+uyjH5OKOvuuv\noThZ6xUpb6SsECtvpKyNCDfAFhERaXWaQUQkNtVgi4hkSDXYUo9zfv8wO241ko++Yae8o4iUmmqw\nRUREWoTOYIvEFm6AHa2erxnKUq9UTW0uhzK2uWyi9dlFOSYXt3eqBjtHkbJCrLyRsjYi3ABbRESk\n1S1csZZJYzWLiEhUqsEWEcmQarCllrUvb2Sf827i8fOOYrPNLO84IqWmGmwREZEWsHDFWrYbM0yD\na5HAwg2wo9XzNUNZ6pWqqc3lUMY2l020PrsIx+TCFWuZtGXt8pAiZO2PSHkjZYVYeSNlbUS4AbaI\niEgrW7hc9dci0akGW0QkQ6rBllq+cePjDBm0GZ9769S8o4iUnmqwRUREWsDzKzp1BlskuHAD7Gj1\nfM1QlnqlampzOZSxzWUTrc8uwjG5cLlqsPMWKSvEyhspayPCDbBFRERa2cIVa5k8dkTeMUSkAarB\nFhHJkGqwpS8bNm5ij3Nu4JH/OIqhg3UOTCRvqsEWEREJbsnKdYwbOVSDa5Hgwn2Co9XzNUNZ6pWq\nqc3lUMY2l020PjvvY7I/U/TlnbW/IuWNlBVi5Y2UtRHhBtgiIiKtauGKNUzaUvXXItGpBltEJEOq\nwZa+fO9PT7Jq3UbOPHr3vKOICKrBFhERCa/ey6SLSLGFG2BHq+drhrLUK1VTm8uhjG0um2h9dt7H\n5MIVnUxWDXbuImWFWHkjZW1EuAG2iIhIq1q4fA3b6SqOIuGpBltEJEOqwZbeuDvTzrmRmWe/lS02\nH5x3HBFBNdgiIiKhLVv9MkMHb6bBtUgLCDfAjlbP1wxlqVeqpjaXQxnbXDbR+uw8j8nnV3TWPQc2\nxPv8RMobKSvEyhspayPCDbBFRERaUWUObNVfi7QC1WCLiGRINdjSm0vvmMeC5Ws599g9844iIgnV\nYIuIiAT27Etr2GErXcVRpBWEG2BHq+drhrLUK1VTm8uhjG0uMjM7yszmmNlcMzuzh/s/b2aPmtls\nM7vFzHao9ZzR+uw8j8nnlq1h+3H1D7CjfX4i5Y2UFWLljZS1EeEG2CIi0nxmNgi4GDgamAacYGbT\num12P7C/u+8N/Bq4INuUre25ZTqDLdIqVIMtIpKhotZgm9nBwLnufmSyfBaAu5/fy/bTge+5++v7\nel712fXZuMnZ45wbmP2lIxg2ZFDecUQkoRpsERFpxCRgftXygmRdbz4GXJ9qohJZ1L6WcSOGanAt\n0iLCDbCj1fM1Q1nqlaqpzeVQxjYXmPWwrsc/cZrZh4D9gW/0cv9pZjbLzGbNmzfvFe9zW1tboZcv\nvPDCXF7/uZfWsP1WI/r1+K6fi7T/WiVv98x552mlvBdeeGGh8tRaHqhwJSK77babz5kzJ+8YmWpr\na2PGjBl5x8iU2lwOZWxz9BIRM3srcBHwJndfWut5o/XZeR2TV93zHPc9u5xvvG+fuh8T7fMTKW+k\nrBArb6SsMPA+O9wAW/V8IhJZgQfYg4EngLcAC4GZwInu/kjVNtOpfLnxKHd/sp7nVZ9dn6/f8Dgj\nhw7i9MN2zTuKiFRRDbaIiAyYu28ATgduBB4Dfunuj5jZeWZ2bLLZN4AtgF+Z2QNmdk1OcVvOcy+t\nYUo/pugTkWILN8BWDXY5qM3lUMY2F5m7X+fuU919F3f/arLuHHe/Jvn5re4+wd1fm9yO7fsZ4/XZ\neR2Tzy5bzQ5bjezXY6J9fiLljZQVYuWNlLUR4QbYIiIirea5l9awg85gi7QM1WCLiGSoqDXYaVGf\nXduKNS9z6NdvZfa5R2DW02QuIpKXQtZgp3HZXRERkVby9Iur2WH8CA2uRVpIagPstC67G62erxnK\nUq9UTW0uhzK2uWyi9dl5HJPzXljNzuO36Pfjon1+IuWNlBVi5Y2UtRFpnsE+AJjr7vPc/WXgF8Bx\n1Ru4+63uviZZvBuYnGIeERGRwpn34ip23rp/X3AUkWJLc4DdtMvuVl8VrL29Pfer+uR5FaEi5Mli\nuWsS+qLkyWJ5xowZhcqTxXLXuqLkyWK5bEaNGpV3hH7J4wIY815Yzc5b9/8MdqSLdUCsvJGyQqy8\nkbI2IrUvOZrZ+4Aj3f3UZPkk4AB3P6OHbT9EZf7VN7n7ur6eV1+YEZHI9CVH6e7Ib9/ON9+/D3tN\nGpN3FBHppohfclwATKlangw8332j5LK7/wYcW2twDfHq+ZqhjGe91OZyKGObyyZan531Mblxk/PM\nS6sHVCIS7fMTKW+krBArb6SsjUhzgD0T2NXMdjKzocDxwCuu+pVcdvf7VAbXS1PMIiIiUjjPr1jL\nuJFDGTF0cN5RRKSJUp0H28yOAS4EBgGXuftXzew8YJa7X2NmfwReAyxKHvJcrSuD6c+NIhKZSkSk\nWtucpfzwjnn87NSD8o4iIj0YaJ+d6n+Z3f064Lpu686p+vmtab6+iIhIkQ10ij4RKbZwl0qPVs/X\nDGWpV6qmNpdDGdtcNtH67KyPyUam6Iv2+YmUN1JWiJU3UtZGhBtgi4iItIqBTtEnIsWWag12GlTP\nJyKRqQZbqh30n7fwq08ezJRxI/KOIiI9KOI0fSIiItKL9jXrWbVuA5PGDs87iog0WbgBdrR6vmYo\nS71SNbW5HMrY5rKJ1mdneUw+vriDXSdswWab2YAeH+3zEylvpKwQK2+krI0IN8AWERFpBXOWrGT3\nibEuJS8i9VENtohIhlSDLV3+9bcPMXWbLTj59TvlHUVEeqEabBERkUDmLF7JVJ3BFmlJ4QbY0er5\nmqEs9UrV1OZyKGObyyZan53VMenuPLF4JbtPHD3g54j2+YmUN1JWiJU3UtZGhBtgi4iIRLdwxVqG\nDR3EuJFD844iIilQDbaISIZUgy0Atzy2hMv//AxXfuzAvKOISB9Ugy0iIhLE44tXstsE1V+LtKpw\nA+xo9XzNUJZ6pWpqczmUsc1lE63PzuqYfHzxSnbfduD11xDv8xMpb6SsECtvpKyNCDfAFhERiW72\nghXsPXlM3jFEJCWqwRYRyZBqsKV9zXoO+dotzD73SAYN8CqOIpIN1WCLiIgEMHvhCvacNEaDa5EW\nFm6AHa2erxnKUq9UTW0uhzK2uWyi9dlZHJOzF7Sz96TGy0OifX4i5Y2UFWLljZS1EeEG2CIiIpE9\nOH8Fe08Zm3cMEUmRarBFRDKkGmw56D9v4epPHMQOW43MO4qI1KAabBERkYJb2tFJ54aNbD9uRN5R\nRCRF4QbY0er5mqEs9UrV1OZyKGObyyZan532MfnggnZeM2kMZo1/wTHa5ydS3khZIVbeSFkbEW6A\nLSIiEtXMZ5bxuh3H5R1DRFKmGmwRkQypBrvcjv3enZz9tmkcsJMG2SIRqAZbRESkwDo61/PU0lXs\nM0VXcBRpdeEG2NHq+ZqhLPVK1dTmcihjm8smWp+d5jE565ll7DNlLJsPHtSU54v2+YmUN1JWiJU3\nUtZGhBtgi4iIRHT3vGUctPNWeccQkQyoBltEJEOqwS4v1V+LxKMabBERkYJqX7ueuaq/FimNcAPs\naPV8zVCWeqVqanM5lLHNZROtz07rmLztiRc4cKdxTau/hnifn0h5I2WFWHkjZW1EuAG2iIhINDc9\nspgj9pyYdwwRyYhqsEVEMqQa7PJZt2Ejr/vKH/njF97ENqOG5R1HRPpBNdgiIiIFdPe8ZfzDNlto\ncC1SIuEG2NHq+ZqhLPVK1dTmcihjm8smWp+dxjF586PplIdE+/xEyhspK8TKGylrI8INsEVERKJY\nv3ETNz6yhCOmTcg7iohkSDXYIiIZUg12udzw8GJ+dOc8fvXJQ/KOIiIDoBpsERGRgvnFzOc4/nXb\n5x1DRDIWboAdrZ6vGcpSr1RNbS6HMra5bKL12c08JheuWMsD81dwzGu2bdpzVov2+YmUN1JWiJU3\nUtZGhBtgi4iIRHD1Pc9x3D7bMXxo8y4uIyIxqAZbRCRDqsEuh1XrNvCmC27lV588mJ233iLvOCIy\nQKrBFhERKYgr/vwMB++ylQbXIiUVboAdrZ6vGcpSr1RNbS6HMra5bKL12c04Jl9ctY5L75jHF47Y\nrfFAfYj2+YmUN1JWiJU3UtZGhBtgi4iIFNlXrn2U9+43mZ3Gj8w7iojkRDXYIiIZUg12a7vuoUV8\n7frHueFzhzJi6OC844hIg1SDLSIikqPHF3fw7797mO+dOF2Da5GSCzfAjlbP1wxlqVeqpjaXQxnb\nXDbR+uyBHpNPLlnJyZfN5EvH7snek8c2N1Qvon1+IuWNlBVi5Y2UtRGpDrDN7Cgzm2Nmc83szB7u\n39zMrk7u/6uZ7ZhmHhER6Z367IG55sHnOf4Hd/MvR+3Gsftsl3ccESmA1GqwzWwQ8ARwOLAAmAmc\n4O6PVm3zaWBvd/+kmR0PvMvdP9DX85atnk9EWktRa7DVZ/ffg/NXcNGf5jLvxVV85wPTec3kMXlH\nEpEmG2ifnWaR2AHAXHefB2BmvwCOAx6t2uY44Nzk518D3zMz82jfvBQRiU99dh82bNzEc8vW8MSS\nlTy8sIM/PraEjrXr+fgbd+aiE6brao0i8gppDrAnAfOrlhcAB/a2jbtvMLN2YCvgxd6eNFo9XzO0\ntbUxY8aMvGNkSm0uhzK2ucBS6bPnLungo5fPpGsM3jUS7xqS/335lWP0v9/vr1zubf3fHtjz/T29\nVk9ZVnZ0MGrUKAA2bHLa166nfe161ry8kW3HDGO3CaPYY9vRfOWdezF9+y0ZtJn11vTURfv8RMob\nKSvEyhspayPSHGD31Ot0P8tRzzaY2WnAacniOjN7uMFs0Yynj19gLUptLocytjndq48MXGp99o9P\nOSBSn93rMTkPuCvbLLVE+/xEyhspK8TKGykrDLDPTnOAvQCYUrU8GXi+l20WmNlgYAywrPsTufsP\ngB8AmNmsItYvpkltLge1uRzMrKgFyeqziZU3UlaIlTdSVoiVN1JWGHifneYsIjOBXc1sJzMbChwP\nXNNtm2uAjyQ/vxf4Uxlq+URECkh9tohIk6R2BjupzzsduBEYBFzm7o+Y2XnALHe/BvgRcKWZzaVy\nFuT4tPKIiEjv1GeLiDRPqpeacvfrgOu6rTun6udO4H39fNofNCFaNGpzOajN5VDYNqvPBmLljZQV\nYuWNlBVi5Y2UFQaYN7V5sEVEREREyijcpdJFRERERIos1AC71mV8IzKzKWZ2q5k9ZmaPmNlnk/Xj\nzOxmM3sy+XfLZL2Z2XeTfTDbzPbNtwUDZ2aDzOx+M7s2Wd4pufzyk8nlmIcm61vi8sxmNtbMfm1m\njyfv98Gt/j6b2T8lx/XDZnaVmQ1rtffZzC4zs6XV04cO5H01s48k2z9pZh/p6bWKrFb/XLT3t468\nJ5vZC2b2QHI7NY+cSZZXHWPd7i9Mf1FH1hlm1l61X8/pabss9Pb7t9s2Rdq39eQtxP5N+vp7zOzB\nJOt/9LBNYfqEOvP2r09w9xA3Kl+6eQrYGRgKPAhMyztXE9q1LbBv8vMoKpcqngZcAJyZrD8T+Hry\n8zHA9VTmoz0I+GvebWig7Z8Hfg5cmyz/Ejg++fkS4FPJz58GLkl+Ph64Ou/sA2zvFcCpyc9DgbGt\n/D5TuSjJ08Dwqvf35FZ7n4E3AvsCD1et69f7CoyjMs3yOGDL5Oct825bP/ZBzf65SO9vnXlPBr6X\n977t7Rjrdn9h+os6ss7o6vPzvtHL798C79t68hZi/yb7a4vk5yHAX4GDum1TpD6hnrz96hMincH+\n22V83f1loOsyvqG5+yJ3vy/5eSXwGJWByXFUBmQk/74z+fk44CdecTcw1sy2zTh2w8xsMvA24NJk\n2YDDqFx+GV7d5q598WvgLcn2YZjZaCq/eH4E4O4vu/sKWvx9pvJF6uFWmTN5BLCIFnuf3f12Xj0X\ndH/f1yOBm919mbsvB/pe79gAAAlTSURBVG4Gjko/fdPU0z8X6f0N9fukl2OsWmH6izqyFkYfv3+r\nFWnf1pO3EJL9tSpZHJLcun/przB9Qp15+yXSALuny/gW8sAaqOTPI9Op/M9pgrsvgsqHCtgm2axV\n9sOFwL8Am5LlrYAV7r4hWa5u1ysuzwx0XZ45kp2BF4AfW6Us5lIzG0kLv8/uvhD4L+A5KgPrduBe\nWvt97tLf9zX6+11P/iK9v/Xu7/ckZQG/NrMpPdxfFNGOn4OTP8Vfb2Z75h0GXvX7t1oh920feaEg\n+9cqZaAPAEupnEDodd8WoE+oJy/0o0+INMCu6xK9UZnZFsD/Ap9z946+Nu1hXaj9YGZvB5a6+73V\nq3vY1Ou4L4rBVP5s+j/uPh1YTaV0oDfh22yVuuPjgJ2A7YCRwNE9bNpK73MtvbUxetubdpn1jNST\n5Q/Aju6+N/BH/n6mrYiKtG9ruQ/Ywd33AS4Cfpdznlq/fwu3b2vkLcz+dfeN7v5aKleFPcDM9uq2\nSaH2bR15+9UnRBpg13MZ35DMbAiVD8vP3P03yeolXX+GSv5dmqxvhf3weuBYM3uGyp9mD6NyRnts\nUkoAr2zX39psfVyeueAWAAuq/kf8ayoD7lZ+n98KPO3uL7j7euA3wCG09vvcpb/va/T3uz+XWS/C\n+1szr7u/5O7rksUfAvtllG0gwhw/7t7R9ad4r8y7PsTMxueVp5ffv9UKtW9r5S3a/k1yrADaeHXZ\nW5H6hL/pLW9/+4RIA+x6LuMbTlJv9CPgMXf/VtVd1Zck/gjw+6r1H06+2XwQ0N71p+go3P0sd5/s\n7jtSeR//5O4fBG6lcvlleHWbQ1+e2d0XA/PNbLdk1VuAR2nh95lKachBZjYiOc672tyy73OV/r6v\nNwJHmNmWyZn/I5J1UUS7zHrNvN3qbI+lUu9aVGH6CzOb2FVna2YHUBmHvJRTlt5+/1YrzL6tJ29R\n9q+ZbW1mY5Ofh1M54fJ4t80K0yfUk7fffYLn/E3T/tyofJv3CSrf/v63vPM0qU1voPInkdnAA8nt\nGCp1SLcATyb/jvO/f9P14mQfPATsn3cbGmz/DP4+i8jOwD3AXOBXwObJ+mHJ8tzk/p3zzj3Atr4W\nmJW817+jMltES7/PwH8kndTDwJXA5q32PgNXUakxX0/ljMzHBvK+Ah9N2j4XOCXvdg1gP7yqfwbO\nA44t4vtbR97zgUeozDByK7B7wY6xTwKfrHVcFTDr6VX79W7gkByz9vb7t6j7tp68hdi/wN7A/UnW\nh4FzkvWF7BPqzNuvPkFXchQRERERaaJIJSIiIiIiIoWnAbaIiIiISBNpgC0iIiIi0kQaYIuIiIiI\nNJEG2CIiIiIiTaQBtoRlZp80sw/3sH5HM3u4h/UzzOzaXp7rmYFMxm9ml5vZe2ts02Zmc8zs2B7u\nG25mD5jZy3lfDEBEJBozO9fMvpjx6y00s/N6uf9nZras1u8FaX2Da28iUjxmNtjdL8k7Rz980N1n\ndV/p7muB1yZXtRQRaUlmNsjdN+ado0m+7e7/1dMd7v5BM7s84zxSQDqDLYVjZv9uZo+b2c1mdlXX\n2YnkTPB/mtltwGerz1yY2X5m9qCZ/QX4xz6efrSZ/dbMHjWzS8zsVZ8BM/u8mT2c3D5Xtf7DZjY7\neZ0re3jcl5Mz2vpciUhLSf4y+LiZXZH0g782sxHJfW8xs/vN7CEzu8zMNk/WP2Nm55jZncD7zGwX\nM7vBzO41szvMbPdur7FZ8pixVevmmtkEM9vBzG5JXvsWM9u+h4xtZrZ/8vP4rhMXZnaymf3OzP5g\nZk+b2elJP3+/md1tZuOS7frMJ9IfGghIoSSd43uA6cC7gf27bTLW3d/k7t/stv7HwGfc/eAaL3EA\n8AXgNcAuyWtUv/5+wCnAgcBBwMfNbLqZ7Qn8G3CYu+8DfLbb4y4AtqFyBb5NdTVWRCSW3YAfuPve\nQAfwaTMbBlwOfMDdX0PlL+OfqnpMp7u/wd1/AfwAOMPd9wO+CPx39ZMnfefvgXcBmNmBwDPuvgT4\nHvCT5LV/Bny3n9n3Ak6k8jvgq8Aad58O/AXoKjXsM59If2iALUXzBuD37r7W3VcCf+h2/9XdH2Bm\nY6gMvG9LVr3q7HKVe9x9XvKnyquS1+v++r9199Xuvgr4DXAocBjwa3d/EcDdl1U95t+T1/+E69Ko\nItK65rv7XcnPP6XSX+4GPO3uTyTrrwDeWPWYqwHMbAvgEOBXZvYA8H1g2x5e42rgA8nPx/P3Pv9g\n4OfJz1fy6r67llvdfaW7vwC08/ffLQ8BO/Yjn0hdVIMtRWM17l/dy2PqHdh23677cm+v39drzAT2\nM7Nx3Qber34Ssyn8vWO/JFgduYiUW0/9Z7199mbACnd/bY3t/wL8g5ltDbwT+EqdWQA28PcTh8O6\n3beu6udNVcubqIyF6s33CmZ2IzABmOXup/bnsdLadAZbiuZO4B1mNiw5o/C2Wg9w9xVAu5l1ndH4\nYB+bH2BmOyV10h9IXq/a7cA7zWyEmY2k8qfKO4BbgPeb2VbA/2/v7lmjioIADL9TWAja+APEwlKI\noL2lpZWFaJAQCNhY+QMsLcROAmqRVqIIQcEySRHxg0g0IXYJWGshC0k3FucsbJb9Ui7JXXgfWFh2\nz+UOWwzD7Jxz6c7sVe+BR8C7iDg7JtafmXm5viyuJU2T8xHRHcO7RcmfPygd4Iv181lgrf/CzPwD\n7EXETYAoZgasS+AN8ATYzcxf9asNSkcbSo7vz90A+8CV+v6fTvGYNL4B112v+dziWkdYYKtVMvMz\nsAJsUcYzvlD+zhtnDnhaNzkejFj3gVIMbwN7lETee/9NyjzhJ+Aj8CIzv2bmDmVuby0itijJv/e6\nZeA5sBIRpyeIV5KmzS5wNyK+AeeAxcw8pOTf5Yj4TukID2se3Abmaw7dAW4MWfcSuMPRkcD7wFy9\n9yx9+2Cqx8C9iNgA/ufY00njk8YKR0bVNhFxJjM7dYf6OrBQC9+pFBGrwINBx/T1rNkHrnZnvCWp\nTSLiAvA2My+dcCgnKiIeAp1hx/TVNUuU3+rVccWl9rGDrTZ6VjeZbAKvp7m4rn4DSzHiQTPAKUrn\nR5LUXh1gIUY8aAa4Bhwea1RqHTvYkiRJUoPsYEuSJEkNssCWJEmSGmSBLUmSJDXIAluSJElqkAW2\nJEmS1CALbEmSJKlBfwFH2kKU3OEt+wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot\n",
    "plt.figure(figsize=[12,6])\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(np.linspace(0,nx,nx),phi._value)\n",
    "plt.xlabel('grid block [-]')\n",
    "plt.ylabel('concentration [-]')\n",
    "plt.grid() \n",
    "plt.tight_layout\n",
    "plt.ylim([0,1.2])\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(pv,cout)\n",
    "plt.xlabel('pore volume [-]')\n",
    "plt.ylabel('outlet concentration [-]')\n",
    "plt.ylim([0,1])\n",
    "plt.grid() \n",
    "plt.tight_layout;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "TSVViewer(vars=(phi, phi.grad)).plot(filename=\"myTSV.tsv\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}