{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples on Using MOIRA\n",
    "---\n",
    "__Mokbel Karam, James C. Sutherland and Tony Saad__\n",
    "\n",
    "__Department of Chemical Engineering, University of Utah__\n",
    "\n",
    "Below are the necessary libraries used in the examples below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "# %config InlineBackend.figure_format = 'svg'\n",
    "import matplotlib.pyplot as plt\n",
    "# plt.rcParams['figure.figsize'] = [3.5, 4]\n",
    "plt.rc(\"font\", family='serif')\n",
    "# plt.rc('text', usetex=True)\n",
    "import numpy as np\n",
    "from sympy import symbols, pretty_print\n",
    "from IPython.display import display, Math, clear_output\n",
    "from matplotlib import cm, animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Example Using FTUS on Advection PDE (1D)\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following PDE\n",
    "\\begin{equation}\n",
    "u_t=-a u_x ; a>0\n",
    "\\end{equation}\n",
    "to be solved using forward Euler discretization in time and upwind discretization in space(FTUS). The difference equation is\n",
    "\\begin{equation}\n",
    "\\frac{u_{i}^{n+1}-u_{i}^{n}}{\\Delta t}=-a\\left(\\frac{u_{i}^{n}-u_{i-1}^{n}}{\\Delta x}\\right).\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "init_cell": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial}{\\partial t} u{\\left (x,t \\right )} =  - a \\frac{\\partial}{\\partial x} u{\\left (x,t \\right )} + \\frac{a}{2} \\left(- \\Delta x + \\Delta{t} a\\right) \\frac{\\partial^{2}}{\\partial x^{2}}  u{\\left (x,t \\right )} + \\frac{a}{6} \\left(- \\Delta x^{2} + 3 \\Delta x \\Delta{t} a - 2 \\Delta{t}^{2} a^{2}\\right) \\frac{\\partial^{3}}{\\partial x^{3}}  u{\\left (x,t \\right )}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from src.MOIRA import DifferentialEquation, i, n # import the differential equation and the indicies we need to use\n",
    "\n",
    "a= symbols('a') # define positive a constant for the advection velocity\n",
    "\n",
    "DE1 = DifferentialEquation(dependentVar='u',independentVars=['x']) # construct a one D differential equation with 'x' and 'u' as independent and dependent variables, respectively\n",
    "\n",
    "# method I of constructing the rhs:\n",
    "advectionTerm1 = DE1.expr(points=[1, 0], direction='x', order=1, time=n+1) # construct the first upwind derivative of the dependent variabe 'u' \n",
    "                                                                        # with respect to the independent variable 'x' using the stencil -1, and 0 evaluated at time tn = n\n",
    "\n",
    "# set the rhs of the differential equation\n",
    "DE1.set_rhs(- a * advectionTerm1 )\n",
    "\n",
    "# compute the modified equation up to two terms\n",
    "DE1.modified_equation(nterms=3)\n",
    "\n",
    "# display the modified equation in latex form\n",
    "display(Math(DE1.latex()))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "init_cell": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial}{\\partial t} f{\\left (y,t \\right )} =  - b \\frac{\\partial}{\\partial y} f{\\left (y,t \\right )} + \\frac{b}{2} \\left(\\Delta y - \\Delta{t} b\\right) \\frac{\\partial^{2}}{\\partial y^{2}}  f{\\left (y,t \\right )}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from src.MOIRA import DifferentialEquation, j, n # import the differential equation and the indicies we need to use\n",
    "\n",
    "b = symbols('b') # define positive a constant for the advection velocity\n",
    "\n",
    "DE2 = DifferentialEquation(dependentVar='f',independentVars=['y'], indices=[j])# construct the first upwind derivative of the dependent variabe 'f' \n",
    "                                                                        # with respect to the independent variable 'y' using the stencil -1, and 0 evaluated at time tn = n\n",
    "    \n",
    "\n",
    "# method II of constructing the rhs:\n",
    "advectionTerm2 = (DE2.f(time=n, y=j) - DE2.f(time=n, y = j-1)) / DE2.dy \n",
    "    \n",
    "DE2.set_rhs(- b * advectionTerm2 ) # set the rhs of the Differential equation DE\n",
    "\n",
    "DE2.modified_equation(nterms=2)  # compute the modified equation up to 2 terms\n",
    "\n",
    "display(Math(DE2.latex())) # display the latex form of the modified equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example Using BTCS on Advection-Diffusion PDEs (1D)\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following PDE\n",
    "\\begin{equation}\n",
    "u_t=-cu_x+\\gamma u_{xx}\n",
    "\\end{equation} \n",
    "to be solved using backward Euler discretization in time and central discretization in space(BTCS). The difference equation is\n",
    "\\begin{equation}\n",
    "\\frac{u_{i}^{n+1}-u_{i}^{n}}{\\Delta t}=-c\\left(\\frac{u_{i+1}^{n+1}-u_{i-1}^{n+1}}{2\\Delta x}\\right)+ \\gamma \\left(\\frac{u_{i+1}^{n+1}-2u_i^{n+1}+u_{i-1}^{n+1}}{\\Delta x^2}\\right).\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "init_cell": true,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial}{\\partial t} u{\\left (x,t \\right )} =  - c \\frac{\\partial}{\\partial x} u{\\left (x,t \\right )} + \\left(\\frac{\\Delta{t} c^{2}}{2} + \\gamma\\right) \\frac{\\partial^{2}}{\\partial x^{2}}  u{\\left (x,t \\right )} - c \\left(\\frac{\\Delta x^{2}}{6} + \\frac{\\Delta{t}^{2} c^{2}}{3} + \\Delta{t} \\gamma\\right) \\frac{\\partial^{3}}{\\partial x^{3}}  u{\\left (x,t \\right )}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from src.MOIRA import DifferentialEquation, i, n # import the differential equation and the indicies we need to use\n",
    "\n",
    "c, g = symbols('c gamma') # define positive a constant for the advection velocity\n",
    "\n",
    "DE3 = DifferentialEquation(dependentVar='u',independentVars=['x'], indices=[i])# construct the first upwind derivative of the dependent variabe 'f' \n",
    "                                                                        # with respect to the independent variable 'y' using the stencil -1, and 0 evaluated at time tn = n\n",
    "    \n",
    "\n",
    "# method II of constructing the rhs:\n",
    "advectionTerm3 = (DE3.u(time=n+1, x=i+1) - DE3.u(time=n+1, x=i-1)) / DE3.dx /2\n",
    "\n",
    "diffusionTerm3 = (DE3.u(time=n+1, x=i+1) - 2* DE3.u(time=n+1, x=i) + DE3.u(time=n+1, x=i-1)) / DE3.dx / DE3.dx\n",
    "    \n",
    "DE3.set_rhs(- c * advectionTerm3 + g*diffusionTerm3 ) # set the rhs of the Differential equation DE\n",
    "\n",
    "DE3.modified_equation(nterms=3)  # compute the modified equation up to 2 terms\n",
    "\n",
    "display(Math(DE3.latex())) # display the latex form of the modified equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example Using FTUS on Advection PDEs (2D)\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the advection equation in a two dimensional space\n",
    "\\begin{equation}\n",
    "u_{t}=-bu_{x}-cu_{y};\\quad(b,c)\\in\\mathbb{R}^{+}\n",
    "\\end{equation}\n",
    "to be discretized using a forward Euler method in time and UPWIND in space (FTUS)\n",
    "\\begin{equation}\n",
    "\\frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{\\Delta t}=-b\\frac{u_{i,j}^{n}-u_{i-1,j}^{n}}{\\Delta x}-c\\frac{u_{i,j}^{n}-u_{i,j-1}^{n}}{\\Delta y},\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with a modified equation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "init_cell": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\partial}{\\partial t} u{\\left (x,y,t \\right )} =  - c \\frac{\\partial}{\\partial y} u{\\left (x,y,t \\right )} - b \\frac{\\partial}{\\partial x} u{\\left (x,y,t \\right )} + \\frac{c}{2} \\left(\\Delta y - \\Delta{t} c\\right) \\frac{\\partial^{2}}{\\partial y^{2}}  u{\\left (x,y,t \\right )} - \\Delta{t} b c \\frac{\\partial^{2}}{\\partial x\\partial y}  u{\\left (x,y,t \\right )} + \\frac{b}{2} \\left(\\Delta x - \\Delta{t} b\\right) \\frac{\\partial^{2}}{\\partial x^{2}}  u{\\left (x,y,t \\right )}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from src.MOIRA import DifferentialEquation, i, n # import the differential equation and the indicies we need to use\n",
    "\n",
    "b, c= symbols('b c') # define positive a constant for the advection velocity\n",
    "\n",
    "DE1 = DifferentialEquation(dependentVar='u',independentVars=['x','y']) # construct a one D differential equation with 'x' and 'u' as independent and dependent variables, respectively\n",
    "\n",
    "# method I of constructing the rhs:\n",
    "advectionTerm1 = - b * DE1.expr(points=[-1, 0], direction='x', order=1, time=n)  - c * DE1.expr(points=[-1, 0], direction='y', order=1, time=n)\n",
    "\n",
    "# set the rhs of the differential equation\n",
    "DE1.set_rhs( advectionTerm1 )\n",
    "\n",
    "# compute the modified equation up to two terms\n",
    "DE1.modified_equation(nterms=2)\n",
    "\n",
    "# display the modified equation in latex form\n",
    "display(Math(DE1.latex()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "nx = 64\n",
    "ny = 64\n",
    "Lx = 1.0\n",
    "Ly = 1.0\n",
    "x = np.linspace(0,Lx,nx)\n",
    "y = np.linspace(0,Ly,ny)\n",
    "dx = Lx/(nx-1)\n",
    "dy = Ly/(ny-1)\n",
    "# create a grid of coordinates\n",
    "xx,yy = np.meshgrid(x,y)\n",
    "\n",
    "u0 = lambda x,y : np.exp(-(x-0.2)**2/0.01)*np.exp(-(y-0.2)**2/0.01)\n",
    "\n",
    "cx = 1.0\n",
    "cy = 1.0\n",
    "\n",
    "dt = 1/128\n",
    "tend = 0.5 #s\n",
    "t = 0\n",
    "\n",
    "cflx = cx * dt/dx\n",
    "cfly = cy * dt/dy\n",
    "\n",
    "# setup the initial condition\n",
    "sol = []\n",
    "u = np.zeros([ny+2,nx+2]) # we will ghost cells to simplify the implementation of periodic BCs\n",
    "u[1:-1, 1:-1] = u0(xx,yy)\n",
    "# set periodic boundaries\n",
    "u[:,0] = u[:,-3] #x-minus face\n",
    "u[:,-1] = u[:,2] #x-plus face\n",
    "u[0,:] = u[-3,:] #y-minus face\n",
    "u[-1,:] = u[2,:] #y-plus face\n",
    "sol.append(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "init_cell": true
   },
   "outputs": [],
   "source": [
    "while t < tend:\n",
    "    un = sol[-1]\n",
    "    unew = un.copy()\n",
    "    unew[1:-1,1:-1] = un[1:-1,1:-1] - cflx * (un[1:-1,1:-1] - un[1:-1,:-2]) - cfly * (un[1:-1,1:-1] - un[:-2,1:-1])\n",
    "    # set periodic boundaries\n",
    "    unew[:,0] = unew[:,-3] #x-minus face\n",
    "    unew[:,-1] = unew[:,2] #x-plus face\n",
    "    unew[0,:] = unew[-3,:] #y-minus face\n",
    "    unew[-1,:] = unew[2,:] #y-plus face\n",
    "\n",
    "    sol.append(unew)\n",
    "    t += dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "init_cell": true,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAHmCAYAAACWKUEfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+0XFV99/HPNwn5SYIEEhDCL0l4bDQ/0FgMlAc0IGj8\nUdAuKFpKsYBo24XaUrDWSrGIlFLlqUqzRJGCKCKiFkRjrKyCEQk1ITVCAii2GJKQgAnc3JCQ7/NH\nzg3DMPfec+acs88+c96vtbKamTszZ9+Bet7sveeMubsAAAAgjah6AAAAALEgjAAAABKEEQAAQIIw\nAgAASBBGAAAACcIIAAAgESyMzGx/M/uCmd03yM/Hmtm/mNnFZvZFMzsi1NgAAEBcquqGkDNGvyfp\nW5JskJ9fIOnX7v5JSf8s6dpQAwMAANGppBuChZG73yJpyxAPWShpafLYlZLmmNmkEGMDAABxqaob\nRuV9gQJN1YvfgM3JfZvbH2hm50o6V5L2mLDHa/d95T5BBggAQNXW3v/Ek+4+JdTxXm3mz5Twuo9J\nP5fU33LXIndflOElUndDFjGF0XpJE1tuT0rue4nkjVskSQfMe7n/6bKzyx8dAAARuNQueyzk8Z6R\n9NESXvccqd/d5+V4idTdkEWln0ozs8kt0163S5qf3D9L0gp3z1V9AACgd4TohpCfSjtO0h9JermZ\nfdTMxkm6SNL7k4d8RtIhZvZRSR+W9N5QYwMAAHGpqhuCLaW5+12S7mq7+8KWn2+V9IFQ4wEAAPGq\nqhu4wCMAAECCMAIAAEgQRgAAAAnCCAAAIEEYAQAAJAgjAACABGEEAACQIIwAAAAShBEAAECCMAIA\nAEgQRgAAAAnCCAAAIEEYAQAAJAgjAACABGEEAACQIIwAAAAShBEAAECCMAIAAEiMqnoAAAAgXhMl\nHVf1IAJixggAACDBjBGAxtrRv0OfnfF5vf+h92mP8XtUNobb3/ddbXxwo3bu2Kk3XHa8Dn/TKzo+\n9lc/ekx3vO+72nP/PXffN+/PXquZ7/qdUMMFeh5hBKBRrj/+Bs05a7bmnDVbo8aO0nkrz6ksiiTp\nro//p+TS2T85SxtXb9QXX/9lnf+Lc7Xnfnt2fPwxFx2tOWfNDjxKoDlYSgPQaGNfNrayY/tO18++\nsFxz3ztHkrTPEfto/yP308obfl7ZmICmY8YIQGMsufg/9MTydXrm8h9rxXUPaN2K9dq+dbvOuPN0\nHXT0NN34ppv02F2/1sn/8iat+c7DeuqRp/TWaxfqN/et1YO3PiTf6Trt2+/ShCkTJElr71+r73/w\nB5JJI0aN0Js/e5L2feW+qcfz1KNPaevGrdr3lfvsvm/KzH21dtnaQZ/z4G2rteLLD2jnjp06/KRX\n6JiLjtaIUSPkO13f/cCdWr9yg2yEafIRk3XSZ07U6Amju3/DgAYijAA0xoJPvkGPL31891KaJF19\n6GclSSNHj9SZP3qPLrXL9Py253XGnafr3k//VLe9+1t659dP1fwPH6Wvvv3rWn7tCh1z0dHq/22/\nvnLyV3Xq107RYW88VGtuf1hfe8ctev8vzpONMH3vgsVat3xdx3HsN3c/nfTpE/XsumclSWP2GrP7\nZ2NeNlYbVj3Z8Xlj9hqjafMP1PwPH6Ud/Tv01bferP6n+nXiP52gh+98RE//6rc66+4zJUk3n3KL\n+jb0EUZARoQRALQ57MTDJElTXj1F2/t2aNrrD5Qk7Td7ip569GlJ0pp/f1ij9xytw954qCRpxsLp\n+uYZ39Lj9z6uafOn6aRPn5j6eGb24ju88+NefuT+evmR+0uSRu85Wkf/9Xzd8q5bdcKVCzRu77Fa\nv3KDHl38Sx224FCdetPva+TokRl+awASYQQALzFm4q5ZlhGjRmj0xBdmXEaMGqHnn3tekrT5f7do\n66Z+XX/8Dbt/Pn7KePVt3Jr6OOOn7lqS63+6XxOSv297ul/jp45P9fxJB0/S9r7t6tvQp2nzp2nh\nojfrx59aqu+c/e96zXlH6piLj5bJhn8hALsRRgDQhUkHTdLEaRN15o/es/u+bZu3adTYXf+zmmYp\nbfLhe2vc5HHa+NDG3WG0YdWTmv6W6R2f99Or79Nrzj1y9zGeXfesRo4eqXH7jFP/b/t16PGHaMZb\npmvTI0/pppO/qokHTtTcP5lT5K8N9DzCCECjjJ44Wtv7tmvjmk362aKfdf06R7x1uhZ/cLF+c99v\ndMDrDtBzzz6n699wo8648zRNmDIh1VKajTAdec5cLf/iAzr42IO1cc0mrVu+Xqfc+A5JUt+Tfbr5\nlFt02rf+QOMmj9Pa/3pCP7/5F5pz5iztfH6n7vt/y/TqM16lESNH6KFvrlb/0/066oLf1eTD99bE\naRPlzw+yJgdgUIQRgEaZe/YcLbnoP7Tiuge05fEt6tu4Vd+/YLHedu1C/fDiH0mSbj39Nr39y2/T\n9y9YrGeeeEa3n3eHpi+crhXXrdSO/h36yVX36vUfOkp/eMdpWvzhJXJ3yaXjLjl29yfW0jru48fq\n9vd9V198/XXauWOnTrnpHbsv4Lh96w5tfHCjtvdt17jJ4/Sac+bq7n/4sVZ8aYWee2a79j9yP514\n5QJJ0oHzD9QP/uqHWv2dNXpuy3OaOnuqZp85q9D3DmgCc6/3f1EcMO/l/qfLzq56GAAABHGpXXa/\nu88LdbxZZn5rCa97hBT090iLCzwCAAAkCCMAAIAEe4wA9KSZWlXI66zSzEJeB0A9EEYAaq+oCMr6\n2kQT0HsIIwC1U2YIZdE+DkIJqD/CCEAtxBJDQyGUgPojjABEqQ4hNBxCCagfwghAdHohijohlID4\nEUYAotGrQTSYgd+XQALiQRgBqFzTgqhd6+9PJAHVIowAVKbpQdQJs0iIzZgx0oyDS3jhNSW8ZgEI\nIwDBEUTDYxYJqAZfCQIgKKIou5laxfsGBEIYAQiGk3s+BBJQPsIIQBCc0ItDIAHlYY8RgFJxAi8P\nG7WB4jFjBKA0RFEYzCABxSGMAJSCE3V4BBKQH2EEoHCcnKtFIAHdI4wAoEcRR0B2hBGAQnEyjguz\nR0A2hBGAwnACjheBBKRDGAEoBCfdeiCQgKERRgDQQMQR0BlhBCA3TrL1xOwR8FKEEYBcOLHWH/8M\ngRfwlSAAGq3oKKjr13Pw9SLALoQRgK7VbaYhxHg7HaNOsTFTq2o1XqBohBGAnld1wLUfP/bwII7Q\nZIQRgK5UHRvDiXl8rWOLNUBYWkNTEUYAekrMQdRJ7JHE7BGahjACkFmM8RHjmLKKNZKIIzQJH9cH\nUHu9EEXtYrvGUGzjAcrCjBGA2mrCiTq2vT7MHjXQBEnzSnjdNSW8ZgGYMQKQSSwxEss4QhmYsYnh\n945hDEBZCCMAtdP0E3MMgVT18YGyEEYAaoUT8guqfi9iCDSgaIQRgNSqPglWffwYxRAnVR8fKBJh\nBKAWOPkOrepA4p8PegVhBCB6nHTTqzKQ+OeEXsDH9QFELYaT7axNqzM/Z+XkI0oYSXpVfayej/Oj\n7ggjAKnEECihdBNCaV8jZDBVdQ0k4gh1RhgBiFbIGCsihro5TohQqiKQiCPUFWEEIEqhoihUEKU5\nftmRFDpWiCPUEWEEoJGqDqJOQkRS6Nmj2L7SBBgOn0oDEJ0yZ4tmbVodZRS1K3ucofeMNWmPGuqN\nMALQGHUIonZlBhJxBLwUS2kAolLWybOOUdSqrGW2KpbWWFZDzIKGkZmdIOlUSeslubtf0vbzwyRd\nKek+SXMlfcXdvx1yjABequ7/pV/3KGo38PsUHUgECxBwKc3Mxku6RtIH3f3jkmab2YK2h10o6W53\nv1zSpyT9U6jxAehNvRZFrYr+3UIFcN1DG70t5B6j+ZIec/dtye17JC1se8w6SVOSv0+RdH+gsQGI\nQNEnzF6OogFF70EijtB0IcNoqqQtLbc3J/e1ukrSUWZ2laSPSfpSpxcys3PNbJmZLevb0FfKYAHU\nWxOiqFWRgUQcoclC7jFaL2liy+1JyX2trpP0BXe/ycymSFpjZq9w902tD3L3RZIWSdIB817u5Q0Z\nQB2FjCJbnu5xPrfccQyYtWl1IXuPQm3KZm8TYhMyjJZKOsTMxiTLacdI+pyZTZa0w903SzpI0trk\n8U9J2ikuKQA0QlGzB2VHUdoQSvu8MoKpyM3ZIcKFOEJMgoWRu/eZ2fmSrjazDZIecPclZnaFpE2S\nLpf0QUkXmNnRkg6T9BF3fzLUGAGgk25jKOtrFx1JRc4eEUdoiqAf13f3xZIWt913Ycvf75Z0d8gx\nAegdZcwWlRlFQx2rqEgqavaIOEJTcIFHAOggZBANdfwiA4k4QlfGadeVBYt2UwmvWQD27wDoCUXN\nFtny6qOoVZHjKeKTayE+Scan1VAlwggAEjEFUbuBQCpijHWII6AqhBGAysVwoo05itoVEUixx1EM\n/06gmQgjALWX9yRfpyhqRRwBxSOMAKDG8s4eEUfAixFGAIbVy58SqutsUTviCCgGYQQAPSLP7FHe\nT6wRL+gVhBGAWstzMg8yW7Ss7U8AVc0elRlHhBdC4QKPABqplChKEz5DPWZeUQPJd4HIPBeDLPMC\njVz8ESEwYwQAeRQ5G1TCzFKepbVuMXOEOiOMAKBbZS6NFRhJVcQRUFeEEYDGKWQZLdB+od3Hynm8\n0HHErBHqijACkAp7OyJQURx1izhCHRFGAJBVyNmiTsfOcfxu4oj9RqiKmZ1gZp8zs4+b2d91+PlZ\nZvYTM/tR8ueP8h6TMAJQuVrNRlUZRa1yBFIvxRF6l5mNl3SNpA+6+8clzTazBR0eerq7H5/8+be8\nxyWMACCtWKKoVY44yhpIMW7GJrp62nxJj7n7tuT2PZIWdnjcn5nZX5rZx8xsct6DEkYAGqXrfTYx\nRtGAgLNHbMZGgfY1s2Utf85t+/lUSVtabm9O7mt1l6RPufuV2vX/BV/POyjCCEBqtVryaqJA8RZj\nHKGWnnT3eS1/FrX9fL2kiS23JyX37ebuv3T3DcnNH0o6zsxG5hkUYQQAvaSLOAq556isOCK6etJS\nSYeY2Zjk9jGSbjezyWY2SZLM7JNmNvAtHjMk/crdn89zUL4SBACGE/MyWifLlPnrRWx5d18fEhO+\nMqS3uHufmZ0v6Woz2yDpAXdfYmZXSNok6XJJT0j6vJn9UtIsSe/Je1zCCAB6UYA46vY71QgYpOXu\niyUtbrvvwpa/f6boY7KUBiCTsk5onChLEGBZjSU19BpmjADU2srJR0T5MfJoRDxzhJqYoMz/DtUZ\nM0YAGiWafTTLB/lThkAbsrNi1ggxIowAILShoqOsUCp5A3lss3bEEbpFGAHIjH1GOWSNnSIDKWMc\nhdhvRMAgNoQRANRBTeKoGyypISaEEYDaq83G37yRUdTsUYnLarEtqQFZEUYAusJyWkZFzrwEmMVp\nxawRmoQwAtA40XwyLY+8sVLikhqzRqgzwghA12Ka3Yl6Oa2sGZfAcZRFTBuxmTVCFoQRANRZwGW1\nOi+pAWkRRgCiE9NMlKT4r/qbJ1gimzUqC8GFtAgjALnEFDFZltN6Yp9RUTLEUYhZI6BKhBGAKMUU\nXJLyzRqFiLBIg4W9RqgbwghAbtFFTEo9N2sUaEmNWSP0MsIIQLS6Ca6oP53WUMwaoU4IIwBIK/bl\nNIlZIyAnwghAIWJaTmMTdnyYNUJdEEYAohZTcEnq/VmjDJg1Qi8ijAAUJqaIYdaoSyVe16gbzO4g\nNMIIQPRiCi5J9Zg1ihAXfEQdEEYAChVTxEQ7a9RDccRyGnoNYQSgFmIKLknxf01IHg1ZTmPWCJ0Q\nRgAKF1PEMGsUl5iW04BORlU9AABIa5Vmlv5f+T43w/LQPOWbXZmraL/KIwtbXn5UztSqUoK7rNft\nKeOa9QEFZowAlCKmk02pV8POu6TWoBMOUAeEEYBaCRFcmf/ruIg4KjKQIo8tltMQM8IIQGnKipgQ\n36FWydJB3mMWHVhAAxFGABoj6iW1Ad3GTQOCiE+nIQTCCECpYpo1yir4klqruUoXSZHMEnE9I/QK\nwghA6eq8EbvSOBowV51DKYIg6hb7jBArwghAbXUbXKUuqUnlX/wxklmiXsJyGgYQRgCCaNSSmtTb\nV8auEAGDshFGABopyKfUiCOgdggjAMHENmtEHAFoRxgBCIo4qoG6jbcgLNNBIowA9JDo46ihwQHU\nCWEEILiYPr7fra6vjE0c5cbMDspEGAGoRN2X1KQejaOYxwYEQBgB6Dm1iaMeipBKvluuBMxGgTAC\nUJkyl9RqEUdSXHEU01iAihBGACoV436jSuKoyiip+vhARAgjAD0rT3QFjyOpmkAhiIAXIYwAVC7G\nJTWpojiSwgUSUQS8BGEEIAq9FkeFBlIZAVPga/bKxusBbMBuNsIIQCOEjiOp4GAoKpLYTwQMaVTV\nAwCAAas0s9T/Ws/z+isnH6FZm1Znfp7PlWx5V4ccXERh02uzRQBhBCAqZcdRHnniSCohkCqWJ4q6\nnYVDeFtHjdHKyYeU8MrZ/38pBJbSAEQn1v1GUr4TOrMrQPwIIwBR6uU46oVA6oXfAeiEMALQSFXG\nkVTvsKjz2IHhEEYAolX2VbGLiCNmj4DeQhgBiFrscSQVM3tUl0AqYpxsvEbMCCMA0WtCHEnxB1Is\nY4vx+/XQOwgjALXQlDiS4gukIsfDbBFiRxgBqI26xFHRgVRVJMUWaEAIhBEAtFilmVHNHg0IGUll\nHYfZItRB0Ctfm9kJkk6VtF6Su/slbT83SX+e3DxU0svc/eyQYwQQt1BXxi7iOAMh0M3VsofSHi1F\nXFGbmSFgl2BhZGbjJV0j6VXuvs3MvmFmC9x9ScvD3iPpaXe/PnnO7FDjA1AfdYojqfuvEkkrTdS0\nxlMVEcRsEeoi5FLafEmPufu25PY9kha2Pebdkiab2V+Y2WWSngk4PgA1EuqTSUUdp8i9R92ocr9S\nkb83n0hD2UKG0VRJW1pub07ua3WIpEnufrWk6yTdaWYj21/IzM41s2VmtqxvQ19Z4wUQubrFkdS8\nmZOm/b6ov5BhtF7SxJbbk5L7Wm2WdK8kufvq5DEHtb+Quy9y93nuPm/8lPElDRdAHdQ1jpoQDEX/\njswWIYSQYbRU0iFmNia5fYyk281ssplNSu5bIukVkpTcN1LSEwHHCKCGQsYRgZROr/5e6H3Bwsjd\n+ySdL+lqM/uEpAeSjdcXSXp/8rBPSZprZh+R9M+S/tjd+0ONEUB9hZxNKPpYvRZIvfS7oHmCflzf\n3RdLWtx234Utf/+tpPNCjglA7wj1abWyjlXWx/tDKiuKWEZDKEHDCADKFjqOJJUWSFK9IomZIvQC\nrnwNoOeEnl0o83h1WGYre4zMFiEkZowA9KSQM0cDx5OKnz0aEOMsUuzB1i1CrNkIIwA9K3QchTpm\nlZEUOoaIFIRGGAHoaVXFkVTe7FGr9lApI5SqmhkiilAFwghAz6sijqo67lARkyaaenV5DN3bqnEl\nRWocS8LtCCMAjVBlHElhZo+GU6foqWq2iFkq8Kk0AI1R9JWr63LsuuF9QpUIIwCNU+WJl5P+0Hh/\nUDXCCEAjVR1HBMBL8Z4gBoQRgMaq+kRMIL0ghvchhjGgeoQRgEaL4WTY9EBq8u+O+BBGABovljCJ\nZRwhNe33Rfz4uD4AJKr6SH+ncQyIYTxliC2IYhsPqsOMEQC0iO0E2WuzSL32+6D3MGMEAG1iuijj\ngF6YRYo1iGIdF6pBGAHAIGJZWmtXt0giPFAnhBEADCHWOBoQayTVJYbqMk6EQxgBwDBiXFrrpNNJ\nPuSYiQz0AsIIAFKKffaok8FiJc/v0SsB1Cu/B4pFGAFABnWZPRoOUQB0xsf1AaALhEW98c8PgyGM\nAKBLXJOnnvhnhqEQRgCQEydaoHcQRgBQAGaP6oF/RhgOYQQABSKQ4sU/F6RBGAFACTgJx4V/HkiL\nj+sDQEl65aP9dUcU5dOvsSW9h7eV8Jr5MWMEACVjeQ2oD2aMACAQZpDCIkbRDWaMACAwZpDKx/uL\nbhFGAFARAqkcvKfIgzACgIoRSMXhfURehBEARIJAyof3DkVg8zUARKb1BM9G7eERRCgSYQQAESOS\nhkYUoWipwsjMvu3uby97MACAwfFx/xcQRM1gZidIOlXSeknu7pe0/XyspCslPS5phqTL3X11nmOm\nnTH6v2Z2t6SbJd3g7pvyHBQA0L0mzyIRRM1hZuMlXSPpVe6+zcy+YWYL3H1Jy8MukPRrd7/CzGZJ\nulbSsXmOm3bz9dclnSzpWUnfMLOvmdnJZmZ5Dg4AyGdgw3avB0MTfscG2tfMlrX8Obft5/MlPebu\n25Lb90ha2PaYhZKWSpK7r5Q0x8wm5RlUqhkjdz8n+eu1kq41s1dqV6UtMrMvS/qSuz+aZyAAgHx6\ncSaJGOppT7r7vCF+PlXSlpbbm5P70jxmc7eDSrvHaPfUlZn9nqSzJf2BpO2S9pV0lZmNlnSxu6/o\ndjAAgGK0B0VdQokQQov1kia23J6U3Jf1MZmk3WN0hZndrF1BdLikJZLOkfTNgSkuMztU0tckHZVn\nQACA4nUKjlhiiRjCIJZKOsTMxiStcYykz5nZZEk73H2zpNu1a8ntP5M9RiuS+7uWNoyOlLSXpC9L\nus7d/6fDY/aQtF+ewQAAwhksSMoKJgIIWbh7n5mdL+lqM9sg6QF3X2JmV0jaJOlySZ+RdKWZfVTS\ndEnvzXvctGF0j7sPt8t7rqR/zDkeAEDFCBjEwt0XS1rcdt+FLX/fKukDRR4z7ebrYT/65u5fzz8c\nAACA6vBdaQAAAAnCCAAAIEEYAQAAJAgjAACABGEEAACQSPtxfQAA0EBbNU4rNavqYQTDjBEAAECC\nMAIAAEgQRgAAAAnCCAAAIEEYAQAAJAgjAACABGEEAACQIIwAAAAShBEAAECCMAIAAEgQRgAAAAnC\nCAAAIEEYAQAAJAgjAACABGEEAACQIIwAAAAShBEAAECCMAIAAEgQRgAAAAnCCAAAIEEYAQAAJAgj\nAACAxKiqBwBpR/8OfXbG5/X+h96nPcbvUdkYbn/fd7XxwY3auWOn3nDZ8Tr8Ta8Y9PErb/xv/eSq\nn8pMOuS4g3XClQtkZgFHDAAIYavGaZVmVj2MYJgxqsj1x9+gFdc9IEkaNXaUzlt5TmVRJEl3ffw/\nJZfO/slZOuUr79Ctp9+mZ9Y90/Gx6/97vRZ/eIne/b3T9d6f/onW/tcTWva5+wOPGACA4hFGkRj7\nsrGVHdt3un72heWa+945kqR9jthH+x+5n1be8POOj19+7QpNf8vhGr/veNkI09yz5+j+a34WcsgA\nAJSCpbQKLLn4P/TE8nV65vIfa8V1D2jdivXavnW7zrjzdB109DTd+Kab9Nhdv9bJ//ImrfnOw3rq\nkaf01msX6jf3rdWDtz4k3+k67dvv0oQpEyRJa+9fq+9/8AeSSSNGjdCbP3uS9n3lvqnH89SjT2nr\nxq3a95X77L5vysx9tXbZ2o6P/819a3XE22e86LEbfr5B27du16gxo/TdD9yp9Ss3yEaYJh8xWSd9\n5kSNnjC6y3cLAIBwCKMKLPjkG/T40sc156zZmnPWbEnS1Yd+VpI0cvRInfmj9+hSu0zPb3teZ9x5\nuu799E9127u/pXd+/VTN//BR+urbv67l167QMRcdrf7f9usrJ39Vp37tFB32xkO15vaH9bV33KL3\n/+I82QjT9y5YrHXL13Ucx35z99NJnz5Rz657VpI0Zq8xu3825mVjtWHVkx2f9+y6Z1/yWLnU9+RW\nrV+5Xk//6rc66+4zJUk3n3KL+jb0EUYAgFogjCJ22ImHSZKmvHqKtvft0LTXHyhJ2m/2FD316NOS\npDX//rBG7zlah73xUEnSjIXT9c0zvqXH731c0+ZP00mfPjH18V6yedozPFaS3DVu77Fav3KDHl38\nSx224FCdetPva+TokanHAABAlYKGkZmdIOlUSeslubtfMsjj3i3pBkkT3b3zDuAGGDNx1yzLiFEj\nNHriCzMuI0aN0PPPPS9J2vy/W7R1U7+uP/6G3T8fP2W8+jZuTX2c8VN3Lcn1P92vCcnftz3dr/FT\nxw/6+P6n+3ff3vZ0v2S7jrvXwXtp4aI368efWqrvnP3ves15R+qYi4+WiU+sAQDiFyyMzGy8pGsk\nvcrdt5nZN8xsgbsvaXvc70gN+lxgTpMOmqSJ0ybqzB+9Z/d92zZv06ixu/7RpllKm3z43ho3eZw2\nPrRxdxhtWPWkpr9lesfnHfC6l2vjQ5t2396w6klNedUU7TFuD/X/tl+HHn+IZrxlujY98pRuOvmr\nmnjgRM39kzlF/coAAJQm5IzRfEmPufu25PY9khZK2h1GSTxdKOk8SR8JOLbgRk8cre1927VxzSb9\nbFH3n+g64q3TtfiDi/Wb+36jA153gJ579jld/4Ybdcadp2nClAmpltJshOnIc+Zq+Rcf0MHHHqyN\nazZp3fL1OuXGd0iS+p7s082n3KLTvvUHGjd5nI780zm68cSb1LexT+P2HqcV1z2g177vSEnSQ99c\nrf6n+3XUBb+ryYfvrYnTJsqfH2JNDgCAiIQMo6mStrTc3pzc1+ofJP29uz831MUCzexcSedK0l4H\nTyp4mGHMPXuOllz0H1px3QPa8vgW9W3cqu9fsFhvu3ahfnjxjyRJt55+m97+5bfp+xcs1jNPPKPb\nz7tD0xdO14rrVmpH/w795Kp79foPHaU/vOM0Lf7wErm75NJxlxy7+xNraR338WN1+/u+qy++/jrt\n3LFTp9z0Du25/56SpO1bd2jjgxu1vW+7xk0ep6mvnqoTrlygG9/0VdkI0yH/9yDNe/9rJUkHzj9Q\nP/irH2r1d9bouS3PaersqZp95qxC3zsAAMpi7mH+a97MFkj6iLsvSG5/SNI0d/9QcvsgSZdKejB5\nyicl/Z2kO9x92WCve8C8l/ufLju71LEDABCLS+2y+919XqjjjZ33Kj9k2U2Fv+5qmxP090gr5IzR\nUkmHmNmYZDntGEmfM7PJkna4+/9IOmvgwWb2SUlXNXnzNQAACCvYla/dvU/S+ZKuNrNPSHog2Xh9\nkaT3Dzxd6uB0AAAV9klEQVTOzKaY2UeTmxea2YGhxggAAJot6Mf13X2xpMVt913YdnuDpE8kfwAA\nAILhu9IAAAAShBEAAECCMAIAAEgQRgAAAAnCCAAAIEEYAQAAJIJ+XB/xm6lVqR63iu/5BQD0IMKo\nYdKGT97XIZwAAHVEGPW4okKoiOMSSwCA2BFGPaaqEEqjdWxEEgDUw7Zt47T60VlVDyMYwqjmYg6h\nobSPm1ACAMSAMKqpugbRYJhNAgDEgDCqkV6LocEQSQCAqhBGkWtKDA1m4PcnkAAAIXCBx4g1PYpa\nzdQq3g8AQOmYMYoQATA4ZpAAAGUijCJCEKVHIAEAysBSWgRYJuoe7xsAoEiEUcU4sedHWAIAisJS\nWkU4kReP5TUAQF7MGFWAKCoX7y8AoFvMGAXECTscZo8AAN0gjAKpQxRlHWMdomOmVtVinACAOBBG\nJYsxiIoa01CvE1OMEEcAgLQIoxLFEEVVjaHTcauME5bWAABpEEYlqTKKYgiyTmL4clhmjwAAQyGM\nSlBFmMQaQ4OpMpKIIwDAYAijgoUMlLrF0GCqWOYijgAAnXAdowIRRfmEvoJ1L76HAIB8CKOChDrJ\nNuHrL0L+jr3+XgIAsiGMChDi5NqEIGoX6ndu2vsKABgce4xyKvukGuqkPWvT6q6et3LyEQWP5KVC\n7AdizxEAQCKMcqlrFHUbQVleq+hgCrFBmzgCgA76JN1nVY8iGMKoS2VGUdGvXWQIdXvMokKp7Hgh\njgCg2QijLtQhiqqIoaEMjKeIQCp79og4AoDmIowyij2KYguidq3jyxtJBAwAoGiEUQSaEESdFDGL\nVFYcEV0A0Ex8XD+DMmaL8r7mrE2raxlFrfL+DmXN4vExfgBoHmaMUootioqOIVve/XN9bjFjyDOD\nVNa+I2aOAKBZCKMUei2K8kRQ2tfLE0uzNq3uenmNkAEA5EEYVaDbKMoTREXHUJbjdRNJeWePiowj\nYgsAmoM9RsMoerYoZBTZ8hf+VCnPOLqNwVj+uQEA6oUZoyHEcHLtNohiNTC2LLNI3c4eMdMDAMiK\nGaNAQkRRDLNDaXUz1qo/fcesEQD0PmaMBlH1STBLBBQSQ8tyPHde90+15dlnj7LMHDFrBADIgjAK\nIGtklR5FeSIozetlDKWsy2tVxhGhBQC9jaW0DoqcLSorirpaNlum4qNoqONkPFaW3yfrslqV/0wB\nAPXBjFGJyoyi1EKEUNrjp5hJyrK0lud6RwAAdMKMUZuqZgMKj6JQs0NZpBxPltmwLDNHzBoBAIbD\njFGLqk6cpURRt9Ieo9srWw+MrcDZoywzR+wRAgAMhTAqQWVR1E0QdfuJtk7PyxJLKQMp5jgisgCg\n9xBGiaJmi2oRRWVd66j1ddNG0jIVFkcAAOTFHqPIDRtFWfYSLVd5UZTnWCnGnyYOq9hvxF4jAOgt\nhFGBip4tShVFaYQMom6PXUEcAQDQjjBS+P/qDxZFVQZRuzRjKSiO0mK2BwDQjj1GBSnyJFtYFGXV\nzebtrF8HslxD7z9KsedoOKGvb8QmbAA97RlJ91Y9iHAaP2MU22xR7ijKureny6tU53qNnL9jkUtq\nzBoBAFoxY1SAtCfX3Ptf0kRREa+TR9rrFA2MdbDZo2FeJ7ZPqjFrBADhmNlkSZdLelTSDEkfcfd1\nHR73K0m/Sm4+7u7vHu61Gz9jFJMhZ0KKiKKQV8NOe6wce4aGmzli1ggAetZlkn7g7pdLuk3SlYM8\n7jp3Pz75M2wUSQ0PoyJOiEXNFuXaVFzWxR+LkDeOYvtaEwBADBZKWpr8/Z7kdifHmtmFZnapmR2d\n5oVZSquDoeKgqI/0lynN8tpQm7KH2JA93JJa2o3YLIUBQHD7mlnrWWqRuy8auGFm35O0X4fnfUzS\nVElbktubJe1tZqPcfUfbYy9295+a2XhJ/2Vmb3X3h4caFGEUga6X0AqOojVrsj1ekmbMyPDgAj5x\nFjPiCgAyedLdBz0ruPtJg/3MzNZLmijpaUmTJD3VIYrk7j9N/m+fmS2XdIykIcOo0UtpeQXbdN1J\nQVG0Zs0Lf7qR+fndht4QzytqrxEAoDZulzQ/+fsxyW2Z2QgzOzj5+wIzO7nlOdMlPTLcCzc2jGLZ\ncJtrw3WO5+WJodyvWUIcFSGWfycAAMP6iKQTzeyjkk6V9JfJ/bOVRJKk9ZLOMbOPmNm/SLrV3e8e\n7oVZSitZ17MVJcVD0TE02DGGXWIballtuItAdlDUXqO8WE4DgPK5+yZJ53S4f7mkWcnfV0p6Z9bX\nbuyMUV5FzC509Um0yKOo9VjDHq+bGaAYNpMDAHoWYdQrIoqiQo7bRTQW8T1qLKcBQLMRRiUqfBmt\nyxN/EVF0V8ufQo/PrBEAICKN3GMUw6xAkd8SP1QoZImitNEz8Ljj0r/00PuOBttv1MVeo6GE/nJZ\nAED9MGPUhVLDqsDZojKiqP053TyvTIUGZ5diCG8AQHcIo5IEu3bOICGVNoqKiJu0r9HVktpgocNy\nGgCgBIRRBQqb1cgRB2XM9uSOowDSBCszPgDQXIRRTApaRhsuPspc/sr12gXNAsWwnAYAqCfCCIUb\nLo4yzxoROgCAQBoXRk1fJolts3Qu7DMCABSscWHUM3Juui5bTwUYAKAxCKOM0sw4FfqJtB5dRho0\n4HpkFqjpM5MAUFeEUWBVbgwOPYvDrBEAoG4aeeXrKPXITElpCr4K9nBmapVWaWa4AwJArPrUqHNU\n0DAysxMknSppvSR390vafv7XkvaXtFa7viTiY+7+YMgx1lks+4sAAKirYGFkZuMlXSPpVe6+zcy+\nYWYL3H1Jy8P2lPQhd3czO03SP0p6W6gx1t2MGcQRAAB5hNxjNF/SY+6+Lbl9j6SFrQ9w9791d28Z\n2zMBxwcAABou5FLaVElbWm5vTu57CTMbLemPJX1gkJ+fK+lcSdrr4EnFjhIAADRWyBmj9ZImttye\nlNz3IkkUfV7S37j7I51eyN0Xufs8d583fsr4UgYLAACaJ2QYLZV0iJmNSW4fI+l2M5tsZpOk3fuQ\n/lXSVe5+v5m9M+D4AABAwwVbSnP3PjM7X9LVZrZB0gPuvsTMrpC0SdLlkm6Q9GpJh5mZJE2Q9I1Q\nYwQG8FF9AGimoB/Xd/fFkha33Xdhy99PDTmeqMxTo64TAQBAjLjydWAe8CKF7Y6r8/EqfN8AAM1B\nGGWUZoll5eQjijtgxiCYMaO4Q+cxXBQNOs55RY+kGizFAUA9NS6MeuaE1WVAhJ41AgCgThoXRihf\n17NFWfXI7BIAIB6EUUwKOtEPFx7RzhoN9vtnXE6sch8XAKDeCKM6GOxEnyOkyoqjaKMrUej+LwBA\nzyGMKhBiRiPNclXREZPm9WLZHD6UntmHBgDIjDAqSeEzE13MGqWNo7yBlPY1hhxP1mU09hcBAErQ\nyDDKOyNQ6oxCwSf8tDM0WQPpOGV7TldR1IUY9hcx4wQA9RX0ytd4gc+VbHnGJ82V1Ok5w1w1e8YM\nac2adIcoY49Q18tnBc8Wsb8IADCcRs4YhdL1ibibE/8wz6lqb8+wx41sSYzZHgBoNsKoQl0t++RY\nKgodR7miqIvZohiW0QAA9UYYdan0mYUSZo2kcHEU26fPQi2jMeMEAPVGGJWslBPyUDMjEcRRqtfv\nZrZoCMwWAQCK0NjN16s0UzO1quphDL0Je6hN1YNtxB7ueYnWeEm7MTvt6w2r2ygqeT8Ssz0A0MEz\nkpZWPYhwmDHKIe2JdLhZo1JmOzJExIwZL/zJoqvnlRBFw71/fBoNAJBWY2eMaqPbWaOB52qI53dQ\n2jLbcKFW8VJYEbNFzDgBQP01esYo5Mkw16xR3v04VX8kPm8UMVsEAAik0WFUK0XEUehASnNMoggA\nEBHCKKBS9xqlfW6IOEobYTmiqEgsowEABhBGBSjypNj1kpq0KzSyzB4VHR5ZXjNnFDFbBAAoQ+PD\nqKioKWqv0bCKmIlpf71uQqn9eVmeW3IUZcFMDwCgFZ9Kq8DKyUdo1qbVg/582C+YTXGdot3xkfWL\nastcvipgo3iaKAo9W0RcAUDvaPyMkRR+1iiNYQMgy5JV1VeFzrLEN4Sio4igAQC0I4wqkuYEniqO\nYg6kLMes+pICXSKuAKC3EEaJKmaNCokjKfv+nrIDKesxUoyf2SIAQAjsMarYcPuNpBR7jqTsV7lu\nD42se5GGeq20UgYdUQQACIUwalHUF8uW8QW1qeJISrcxu5OQy2wZZrhi3Gw9gMACgN7DUlpJil5S\nkzJ8TL2Kq1ynlWGWqIwv1yVmAABDIYzaVHXizBJHtQykDGPJEkRVLaERWADQm1hKK1HWJbU0+40G\npF5ak14cJN0ss3WriyhLG0VZl8+IIgBAGswYdVDlSTTLCb+rpaayvg4k5+tnmQnj6z4AAGVhxiiA\nsmeOpAyzR606xUuWGaWC4qqspbMBzBYBANIijAZR9CfLyowjKWcgtQq4JynrjFfVUQQA6H0spQ2h\n6JNqN8tqWWOgrE9zFambMcYQRUQWAPQ+wiiwbk6u3URBjIHU7ZiIIgBAKCylDaOMizV285oDcZBl\neU16cYjkXmbrQp4463aTNREDAOgWYZRCLHEkZd971CpUJBUxUxVTFBFaAJptm7RjTdWDCIYwqlCe\nOJKyzx61GixesgRTGUt1eT6KTxQBAPIijFIqY9Yo7+vmmT0aTFX7kvJem4iAAQAUgTDKINY4GlB0\nJIVQxMUay4oiYgsAmodPpWUU80m4m4/3V6Woscb8zwMAUD/MGHWhzJkjSblfu4g9SGUpKtzKDBei\nCACaizDqUllxVORrx7DMVsYMFlEEACgLYRSposOrPVDKDKWylvPKjhaiCABAGOVQ5qzRwOtL+ZfW\nOhksXtIGU+i9TEQLACAEwiinsuNo4BhSOYHULrbN26GCiPACAEh8Kq0QIU/eTTmBh/xdm/KeAgCG\nRxgVJOTJtZdP5KHjr5ffSwBAdiylFSjEslrrsQaEOmaZqggUoggA0I4wKljIOGo9plTPQCKIAAAx\nIYxKUFWo1GUWqcowIYoAAEMhjEpUxexR67EHxBBJMQRJDGMAAMSNMCpZlXHUOoZ2Ia6/FJMYxwQA\niA9hFEAMcdSuKaHQlN8TAFAMPq4fSJOuQRQL3m8AQFaEUWCcrMPgfQYAdIOltArU+eP1sSOIAAB5\nMGNUIU7ixeL9BADkxYxRxWLcmF03BBEAlGmLpLuqHkQwhFEEWFrrDkEEACgaYRSR2C7KGCuCCABQ\nFvYYRYqT/0txyQMAQNmYMYoYS2y7EEMAgFAIoxpo6hIbQQQACI0wqplen0UihgAAVSKMaqrXZpEI\nIgBADAijHtAeFbGHEhEEAIgVYdSDOoVHVbFEBAEA6oQwaojBAqXIYCKCAAB1Rxg1HDEDAMALuMAj\nAABAgjACAABIEEYAAAAJwggAACBBGAEAACQIIwAAgARhBAAAkCCMAAAAEkEv8GhmJ0g6VdJ6Se7u\nl7T9fKykKyU9LmmGpMvdfXXIMQIAgLiZ2QhJ50i6VNIb3f2/B3nceyQdKel5SY+4+78O99rBwsjM\nxku6RtKr3H2bmX3DzBa4+5KWh10g6dfufoWZzZJ0raRjQ40RAADUwhxJ90rqG+wBZjZN0l9KOtLd\n3czuM7MfuvuaoV445FLafEmPufu25PY9kha2PWahpKWS5O4rJc0xs0nhhggAAGLn7j9z9+XDPOwk\nSfe7uye3l0p683CvHXIpbaqkLS23Nyf3pXnM5tYHmdm5ks5Nbm671C7rOIWGwuwr6cmqB9EAvM/l\n4z0uH+9x+f5P2MM99j3pnH1LeOGxZras5fYid180cMPMvidpvw7P+5i7fzvF66fpjpcIGUbrJU1s\nuT0puS/rY5S8cYskycyWufu8YoeKVrzHYfA+l4/3uHy8x+Vri4nSufvJIY/XctyTcr7EeknTW25P\nkvTwcE8KuZS2VNIhZjYmuX2MpNvNbHLLctnt2rXkpmSP0Qp33/zSlwIAAHgxMxthZgcnN78n6bVm\nZsnt+ZK+O9xrBJsxcvc+Mztf0tVmtkHSA+6+xMyukLRJ0uWSPiPpSjP7qHZV3ntDjQ8AANSDme0t\n6QOS9pJ0rpl9xd1/Imm2pH+TNMvd/9fMrpT0z2b2vKQvDLfxWpLshT1J9WRm57auSaJ4vMdh8D6X\nj/e4fLzH5eM9LlftwwgAAKAoXPkaAAAgQRgBAAAkgn4lSB58nUj5UrzHfy1pf0lrJc3TrmtJPBh8\noDU23Hvc8rh3S7pB0kR3fybgEHtCin+XTdKfJzcPlfQydz876CBrLsV7fJh2/W/yfZLmSvpKymvP\nIGFm+0v6hKQ57v66Dj/nvFeCWoQRXydSvpTv8Z6SPpRcWv00Sf8o6W1VjLeOUr7HMrPfkTSzkkH2\ngJTv83skPe3u1yfPmV3FWOsq5Xt8oaS73f2fzexISTdLIoyy+T1J39KusOyE814J6rKUxteJlG/Y\n99jd/7bl0uojJDGTkc2w73FywrlQUseZJKSS5n8v3i1pspn9hZldJv5dzirNe7xO0pTk71Mk3R9o\nbD3D3W/Ri6/c3I7zXgnqEkZ5vk4E6aR+/8xstKQ/lvTRAOPqJWne43+Q9Pfu/lywUfWeNO/zIZIm\nufvVkq6TdKeZjQwzvJ6Q5j2+StJRZnaVpI9J+lKgsTUJ570S1GIpTQV+nQgGler9S6Lo85L+xt0f\nCTS2XjHke2xmB0naW9JpL1yoVR8yszvcPehXANRcmn+XN2vXN3PL3Vcn/5V9kKRfhRhgD0jzHl+n\nXRfUu8nMpkhaY2avcPdNgcbYBJz3SlCXGSO+TqR8w77HyTLPv0q6yt3vN7N3VjTWuhryPXb3/3H3\ns9z9cne/PHnMVURRZmn+92KJpFdIUnLfSElPBB9pfaV5jw/Srg9qSNJTknaqPuecaHHeK19tLvBo\nZidKepekDZK2u/slA18n4u6Xm9k47dqdv1a7vk7kMnbnZ5PiPb5V0qsl/SZ5yoROn5TA4IZ7j5PH\nTJF0nqRLkz//6u6PVzXmOkrx7/Jekq6Q9JikwyV9w93vqG7E9ZPiPf497doc/F+SDpN0v7tfU92I\n68fMjpN0pqSTtWum/p+0a/8h570S1SaMAAAAysa0JgAAQIIwAgAASBBGAAAACcIIAAAgQRgBAAAk\nCCMAAIAEYQQAAJAgjAAAABKEEYDUzOxLZtZnZk+Y2SlmdrKZ/dLMfmFmp1Q9PgDIqy5fIgsgAu7+\nJ2b235L+VtIySZskPS7pre7+dKWDA4AC8JUgADIxs5GSfqxd31D/sKTvuftt1Y4KAIpBGAHIzMxe\nLel+SYvd/a1VjwcAisIeIwDd+LWkJyS9zsz2qXowAFAUwghANz4l6XztiqNPVzwWACgMm68BZGJm\nx0sa5e53mNl6ST8xs5vc/Y6KhwYAuTFjBCA1M/uEpK9Jeo2Z7SXpBEl9kq43sysqHRwAFIDN1wAA\nAAlmjAAAABKEEQAAQIIwAgAASBBGAAAACcIIAAAgQRgBAAAkCCMAAIAEYQQAAJD4/+IA2yRH7c9Q\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c5c98d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "levs = np.linspace(-1,1,20)\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.contourf(xx,yy,sol[-1][1:-1,1:-1]+sol[0][1:-1,1:-1],cmap=cm.jet,levels=levs,vmax=1.0,vmin=-1.0)\n",
    "plt.text(0.1, 0.4, 'time=0.0s', fontsize=12)\n",
    "plt.text(0.6, 0.92,'time={}s'.format(tend), fontsize=12)\n",
    "plt.xlabel('x', fontsize=14)\n",
    "plt.ylabel('y', fontsize=14)\n",
    "cbar = plt.colorbar()\n",
    "plt.clim(-1,1)\n",
    "cbar.set_ticks(np.linspace(-1,1,5))\n",
    "# plt.savefig('./advection_2d.pdf',transparent=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Solver\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "init_cell": true,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "24e95a88228e47f0920d38f3e93f6404",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Tab(children=(VBox(children=(HBox(children=(FloatText(value=0.01, description='dt:'), FloatText(value=10.0, de…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib notebook\n",
    "from numerical.layout import *\n",
    "plty=PlotStyling()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "63px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}