BoundaryConditions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "sonic-individual",
   "metadata": {},
   "source": [
    "# Boundary Conditions\n",
    "\n",
    "In this workbook we are going to modify our fast checkerboard *SOR* solver so that it can manage boundary conditions.\n",
    "\n",
    "Remember we are solving\n",
    "$$u_{xx}+u_{yy}=0$$\n",
    "The boundary conditons can be\n",
    "* **Dirichlet** boundary conditions where a value of the unknown, _u_, is specified, e.g. $$u(0,y)=273$$\n",
    "* **Neumann** boundary conditions where a value of the derivative, _u<sub>x</sub>_, is specified. e.g. $$u_x(0,y)=0.$$\n",
    "* **Robin** boundary conditions where the value of a function of the unknown and the derivative is specified, e.g. $$u_x(0,y)-2u(0,y)=1.$$\n",
    "\n",
    "## Neumann boundary conditions\n",
    "We have already implemented a **Dirichlet** boundary condition in our solver, and **Robin** boundary conditions are very rare, so let us consider **Neumann** conditions.\n",
    "\n",
    "Supose we have\n",
    "$$u_x(0)=0$$\n",
    "Since we are already using a 2nd order central finite difference approximation, it seems best to do the same for this partial derivative too,\n",
    "$$\\frac{\\partial u}{\\partial x} \\approx \\frac{u_W-u_E}{2\\Delta x}.$$\n",
    "At the boundary\n",
    "$$\\frac{\\partial u}{\\partial x} =0$$ so\n",
    "$$0 = \\frac{u_W-u_E}{2\\Delta x} \\implies u_W=u_E.$$\n",
    "This is the so-called _symmetry_ boundary condition.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "frank-welsh",
   "metadata": {},
   "source": [
    "### updating the Grid Class\n",
    "We are going to need some flags that tell us if a boundary is a Neumann BC or a Dirichlet BC, so the `Grid` class needs extending.\n",
    "\n",
    "We can define some flags DIRICHLET_BC, NEUMANN_BC and BC_NAME within the class, give every grid left, right, top and bottom flags which are set to *Dirichlet* by default and write a little function to specify boundaries and annother to report them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "tutorial-tucson",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "\n",
    "class Grid:\n",
    "    '''Class defining a 2D computational grid.  The grid object\n",
    "    contains is a regular cartesian grid with a single variable, u.\n",
    "    It stores information about the number of grid points in the i\n",
    "    and j directions, the ordinates of these points and the bottom\n",
    "    left corner of the gird (the origin) and the top right corner\n",
    "    (the extent).\n",
    "    \n",
    "    Written by Prof David Ingram, School of Engineering\n",
    "    (c) 2021 The University of Edinburgh\n",
    "    Licensed under CC-BY-NC.'''\n",
    "    \n",
    "    DIRICHLET_BC = 0\n",
    "    NEUMANN_BC = 1\n",
    "\n",
    "    BC_NAME = ['left', 'right', 'top', 'bottom']\n",
    "    \n",
    "    def __init__(self,ni,nj):\n",
    "        # set up information about the grid\n",
    "        self.origin = (0.0, 0.0)  # bottom left\n",
    "        self.extent = (1.0, 1.0)  # top right\n",
    "        self.Ni = ni # grid points in i direction\n",
    "        self.Nj = nj # grid points in j direction\n",
    "        \n",
    "        # initialse x,y and u arrays\n",
    "        self.u = np.zeros((nj, ni))\n",
    "        self.x = np.zeros((nj, ni))\n",
    "        self.y = np.zeros((nj, ni))\n",
    "        \n",
    "        # boundary conditions (left right top and bottom)\n",
    "        self.BC = [self.DIRICHLET_BC, self.DIRICHLET_BC, \n",
    "                   self.DIRICHLET_BC, self.DIRICHLET_BC]\n",
    "\n",
    "    def set_origin(self,x0,y0):\n",
    "        self.origin = (x0, y0)\n",
    "    \n",
    "    def set_extent(self,x1,y1):\n",
    "        self.extent = (x1, y1)\n",
    "        \n",
    "    def generate(self,Quiet=True):\n",
    "        '''generate a uniformly spaced grid covering the domain from the\n",
    "        origin to the extent.  We are going to do this using linspace from\n",
    "        numpy to create lists of x and y ordinates and then the meshgrid\n",
    "        function to turn these into 2D arrays of grid point ordinates.'''\n",
    "        x_ord = np.linspace(self.origin[0], self.extent[0], self.Ni)\n",
    "        y_ord = np.linspace(self.origin[1], self.extent[1], self.Nj)\n",
    "        self.x, self.y = np.meshgrid(x_ord,y_ord)\n",
    "        if not Quiet:\n",
    "            print(self)\n",
    "\n",
    "    def Delta_x(self):\n",
    "        # calculate delta x\n",
    "        return self.x[0,1]-self.x[0,0]\n",
    "    \n",
    "    def Delta_y(self):\n",
    "        # calculate delta y\n",
    "        return self.y[1,0]-self.y[0,0]\n",
    "    \n",
    "    def find(self,point):\n",
    "        '''find the i and j ordinates of the grid cell which contains \n",
    "        the point (x,y).  To do this we calculate the distance from\n",
    "        the point to the origin in the x and y directions and then\n",
    "        divide this by delta x and delta y.  The resulting real ordinates\n",
    "        are converted to indices using the int() function.'''\n",
    "        grid_x = (point[0] - self.origin[0])/self.Delta_x()\n",
    "        grid_y = (point[1] - self.origin[1])/self.Delta_y()\n",
    "        return int(grid_x), int(grid_y)\n",
    "    \n",
    "    def set_Neumann_bc(self,side):\n",
    "        try:\n",
    "             self.BC[self.BC_NAME.index(side)] = self.NEUMANN_BC\n",
    "        except:\n",
    "             print('error {} must be one of {}'.format(side,self.BC_NAME))\n",
    "     \n",
    "    def set_Dirichlet_bc(self,side):\n",
    "        try:\n",
    "            self.BC[self.BC_NAME.index(side)] = self.DIRICHLET_BC\n",
    "        except:\n",
    "            print('error {} must be one of {}'.format(side,self.BC_NAME))\n",
    "    \n",
    "    def report_BC(self):\n",
    "        '''compile a string listing the boundary conditions on each side.\n",
    "        We build up a string of four {side name}: {BC type} pairs and\n",
    "        return it'''\n",
    "        \n",
    "        # initialise the string\n",
    "        string = ''\n",
    "        \n",
    "        # loop over the sides\n",
    "        for side in range(4):\n",
    "            # add the side name\n",
    "            string = string + self.BC_NAME[side]\n",
    "            # and the bounday condition type\n",
    "            if self.BC[side] == self.DIRICHLET_BC:\n",
    "                string = string + ': Dirichlet, '\n",
    "            elif self.BC[side] == self.NEUMANN_BC:\n",
    "                string = string + ': Neumann, '\n",
    "        return string[:-2] +'.' # lose the last comma and space.\n",
    "    \n",
    "    def __str__(self):\n",
    "        # describe the object when asked to print it\n",
    "        describe = 'Uniform {}x{} grid from {} to {}.'.format(self.Ni, self.Nj, self.origin, self.extent)\n",
    "        boundaries = self.report_BC()\n",
    "        return describe + '\\nBoundaries conditions are - ' + boundaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "royal-reunion",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Uniform 11x11 grid from (0.0, 0.0) to (1.0, 1.0).\n",
      "Boundaries conditions are - left: Dirichlet, right: Neumann, top: Dirichlet, bottom: Dirichlet.\n"
     ]
    }
   ],
   "source": [
    "test = Grid(11,11)\n",
    "test.set_Neumann_bc('right')\n",
    "print(test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "critical-snowboard",
   "metadata": {},
   "source": [
    "## Indirect adressing.\n",
    "The SOR solver loops over all the interior cells in `u_new` using whole array opperations.\n",
    "\n",
    "       u_new.flat[centre[::2]] = (1-omega) * u_new.flat[centre[::2]] + \\\n",
    "                                    omega * C_beta*(u_new.flat[north[::2]]+ \\\n",
    "                                    u_new.flat[south[::2]]+ \\\n",
    "                                    beta_sq*(u_new.flat[east[::2]]+ \\\n",
    "                                             u_new.flat[west[::2]]))\n",
    "\n",
    "We call this **indirect** adressing because the locations in `u_new` are given by the entries in the `centre[]`, `north[]`, `south[]`, `east[]` and `west[]` arrays.  The good news is that this means only need to change the cells refered to in the `east[]` array on the left hand boundary if we want to implement the Neumann boundary condition _u<sub>x</sub>=0_.\n",
    "\n",
    "The bit of the SOR code that needs modifying is the loop over the interior points that builds the lists.\n",
    "\n",
    "    for j in range(1,mesh.Nj-1):\n",
    "        for i in range(1,mesh.Ni-1):\n",
    "            k=i+offset*j  # calculate the k value\n",
    "            # save the south, west, centre, east and north points\n",
    "            centre.append(k)\n",
    "            north.append(k+offset)\n",
    "            east.append(k+1)\n",
    "            south.append(k-offset)\n",
    "            west.append(k-1)\n",
    "\n",
    "If we are next to a boundary **and** it's a Neumann boundary condition we need to adjust the pointer to point to the interior cell. So the west cell would be at k+1 and the east cell at k-1.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ordered-prairie",
   "metadata": {},
   "outputs": [],
   "source": [
    "def SOR(mesh,tol=0.5e-7,maxit=10000):\n",
    "    '''Sucessive over Relaxation method for the solution of\n",
    "    the Laplace equation with checkerboarding.  The solver \n",
    "    assumes a uniform Cartesian mesh and uses the optimal\n",
    "    value of the acceleration parameter, omega.\n",
    "    \n",
    "    The method works on the grid stored in the mesh object.  \n",
    "    It will continue itterating until the relative difference \n",
    "    between u^{n+1} and u^n is less than tol. It will also stop \n",
    "    if we have done more than maxit itterations.\n",
    "    \n",
    "    The solution stored in the mesh.u variable is updated\n",
    "            \n",
    "    Written by Prof David Ingram, School of Engineering\n",
    "    (c) 2021 The University of Edinburgh\n",
    "    Licensed under CC-BY-NC.'''\n",
    "\n",
    "    # calculate the optimal value of omega\n",
    "    lamda = (np.cos(np.pi/mesh.Ni)+np.cos(np.pi/mesh.Nj))**2/4\n",
    "    omega = 2/(1+np.sqrt(1-lamda))\n",
    "\n",
    "    # calculate the coefficients\n",
    "    beta = mesh.Delta_x()/mesh.Delta_y()\n",
    "    beta_sq = beta**2\n",
    "    C_beta = 1/(2*(1+beta_sq))\n",
    "    \n",
    "    # initialise u_new \n",
    "    u_new = mesh.u.copy()\n",
    "    \n",
    "    # build the list of k values    \n",
    "    offset = mesh.u.shape[1] # how far away are j+1 and j-1\n",
    "\n",
    "    # create some empty lists\n",
    "    centre = []; north = []; east = []; south=[]; west = []  \n",
    "    # loop over all the interior points\n",
    "    for j in range(1,mesh.Nj-1):\n",
    "        for i in range(1,mesh.Ni-1):\n",
    "            k=i+offset*j  # calculate the k value\n",
    "            \n",
    "            # save the centre point\n",
    "            centre.append(k)\n",
    "            \n",
    "            # do the north boundary checking for a Neumann bounary condition\n",
    "            if (j == mesh.Nj-2) and (mesh.BC[2]== mesh.NEUMANN_BC):\n",
    "                north.append(k-offset)\n",
    "            else:\n",
    "                north.append(k+offset)\n",
    "                \n",
    "            # do the east boundary checking for a Neumann boundary condition\n",
    "            if (i == mesh.Ni-2) and (mesh.BC[1]== mesh.NEUMANN_BC):\n",
    "                east.append(k-1)\n",
    "            else:\n",
    "                east.append(k+1)\n",
    "                \n",
    "            # do the sound boundary checking for a Neumann boundary condition\n",
    "            if (j==1) and (mesh.BC[3]== mesh.NEUMANN_BC):\n",
    "                south.append(k+offset)\n",
    "            else:\n",
    "                south.append(k-offset)\n",
    "                \n",
    "            # do the west boundary checking for a Neumann boundary condition        \n",
    "            if (i==1) and (mesh.BC[0]== mesh.NEUMANN_BC):\n",
    "                west.append(k+1)\n",
    "            else:\n",
    "                west.append(k-1)\n",
    "                \n",
    "    # itteration\n",
    "    for it in range(maxit):\n",
    "        \n",
    "        # red squares [::2] means from 0 to n in steps of 2, remember\n",
    "        # we are taking alternate values from the lists centre, east,\n",
    "        # north, west and south and using them as indicies into the \n",
    "        # u_new array.\n",
    "        u_new.flat[centre[::2]] = (1-omega) * u_new.flat[centre[::2]] + \\\n",
    "                                    omega * C_beta*(u_new.flat[north[::2]]+ \\\n",
    "                                    u_new.flat[south[::2]]+ \\\n",
    "                                    beta_sq*(u_new.flat[east[::2]]+ \\\n",
    "                                             u_new.flat[west[::2]]))\n",
    "        \n",
    "        # black squares [1::2] means from 1 to n in steps of 2\n",
    "        u_new.flat[centre[1::2]] = (1-omega) * u_new.flat[centre[1::2]] + \\\n",
    "                                    omega * C_beta*(u_new.flat[north[1::2]]+ \\\n",
    "                                    u_new.flat[south[1::2]]+ \\\n",
    "                                    beta_sq*(u_new.flat[east[1::2]]+ \\\n",
    "                                             u_new.flat[west[1::2]]))\n",
    "\n",
    "        # compute the difference between the new and old solutions\n",
    "        err = np.max(np.abs(mesh.u-u_new))/np.max(np.abs(mesh.u))\n",
    "        \n",
    "        # update the solution\n",
    "        mesh.u = np.copy(u_new)\n",
    "        \n",
    "        # converged?\n",
    "        if err < tol:\n",
    "            break\n",
    "    \n",
    "    if it+1 == maxit:\n",
    "        print('Checkerboard Gauss-Seidel itteration failed to converge, error = {}'.format(err))\n",
    "    return it+1,err"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "pointed-palestine",
   "metadata": {},
   "source": [
    "## The original test problem\n",
    "Were going to use the same test problem as before to make sure it still works! \n",
    "\n",
    "$$0=u_{xx}+u_{yy}$$\n",
    "subject to the boundary conditions\n",
    "$$\\begin{align*}\n",
    "u(x,0)&=0&0\\le x\\le 2 \\\\\n",
    "u(x,1)&=0&0\\le x\\le 2 \\\\\n",
    "u(0,y)&=0&0\\le y \\le 1\\\\\n",
    "u(2,y)&=\\sin 2\\pi y&0\\le y \\le 1.\n",
    "\\end{align*}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "white-scanning",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Example929(ni,nj):\n",
    "    # set up a mesh\n",
    "    mesh = Grid(ni,nj)\n",
    "    mesh.set_extent(2.0,1.0)\n",
    "    mesh.generate()\n",
    "    # now the RHS boundary condition\n",
    "    mesh.u[:,-1]=np.sin(2*np.pi*mesh.y[:,-1])\n",
    "    return mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "directed-shipping",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Uniform 161x81 grid from (0.0, 0.0) to (2.0, 1.0).\n",
      "Boundaries conditions are - left: Dirichlet, right: Dirichlet, top: Dirichlet, bottom: Dirichlet.\n",
      "Solver took 2.57 seconds.\n",
      "Converged after 517 itterations, final residual is 4.911869522228329e-15\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAh0UlEQVR4nO3de5gcVZ3/8feHQBK5CEmGhUi4LkEIAgFHcMFFhAABVwLquuEHa7KCrBdcV1gUfvweZFF3QZ5d3FVcGQFFRC7iAlFBQC4qYoBRISHBkBBQEkHIBRSBQOD7+6POQKXTPV3dXT3TnXxez1PP1OVU1Znqnk+fOXVpRQRmZta9NhjuCpiZWWsc5GZmXc5BbmbW5RzkZmZdzkFuZtblHORmZl3OQb6OkjRT0l3DXY9GSRolab6k8cNdl04m6TlJO9VY9tprn47nbyRtObQ1tKHkIG+CpMckvZD+mAaGrwx3vcoiaRtJN0haIWmJpI8MUvbdku6S9IykJyVdLGmzZraVnAT8NCKeSOu/S9Idkp6V9FiNOnxS0qOS/izpIUm7pPnjJc2S9HtJIWmHRo9FLZJ2kHSjpJXp9/6KpA1zy/skLZD0qqSZZe13QERsGhGLC5RbBVwKnF52HaxzOMib9570xzQwnDzcFSrRt4FHga2AdwP/JuldNcpuDnweeBOwG7ANcH6T2wL4CHB5bvrPZEF0WrXCkk4ETkjb3hT4G2BZWvwq8CPgfYPsr1lfBZ4CxgOTgXcCH8stfyBN/6rMneY/LBrwHWCGpFFl1sU6SER4aHAAHgOm1Fj2P8D3ctPnAbcBAsYAPwCeBlam8Qm5sneSheLdwHPA94FxwBXAH4H7gB1y5QP4J2AxWXidD2yQls0E7sqV3RW4FVgBLAA+UKP+m6btbpmb1wdcXvDYvBeY28y2gO2AF4ANqyybAjxWMW8D4HHgkDp12jDVY4dByvxlOjb7pOk3pdfpoBrlHwKOzE2fD1xUpdxdwMw69RuXXuuB1/jzFa9dAB8HFgKP5ubtnFt/Vlr/XuBz+fVTmYXAO4f7b8dDewa3yMt3KrBH6qf8a7LW4ozI/po2AL4BbM/roVXZJTMd+Huylu1fAr9I64wlC4/PVpQ/BugF9gGmAR+qrJCkTchC/DvAX6R9fFXSpCr1V8XPgfG31PvFkwOBeU1uaw9gcUSsLrivCWl4i6THU/fKv0pq+H0dEY8AnwG+LWljsmN+WUTcWWOVLwHTJW0saRvgCLLWfzMuJPvPY2tgRhoqHQ3sB1R7zS4EXiT77+BDVHkPkL139mqyftbphvuTpBsHshb5c8AzueHDueX7kbXufgscO8h2JgMrc9N3Amfmpv8DuCk3/R7g/tx0AFNz0x8DbkvjM0mtMuDvgJ9V7Psi4LM16nUX8GVgNNkHxApgQYHjcijZfxq7NLMt4Dhgdo1l1Vrk+6dj8ENgC2AH4OH8a5HK1W2R58rOAuYCc4BRg5TbDfglsDpt+5uAahzLmYNsZwTwMvDm3LxqLfKDK9YLYOfc+rvmlv0ba7fIrwDOGu6/HQ/tGdwib97REbFFbvj6wIKIuIesu0PANQPzU+vtIkm/lfRH4KfAFpJG5Lb7h9z4C1WmN62ox+O58d+SdQlU2h7YL52QfEbSM2ShuXWN3+04YMe07f8h6+deUqPswO/2drIW//sj4uEmt7US2KzGsmpeSD+/GBHPRMRjZB9QRzawjUpfJ/uP4cuRnShcS2rx/wj4X2AToIes2+y8Jva3JdkHTf51fLxKuWrzaq3/2yrlNiNrcNg6yEHeBpI+DowCfg98OrfoVODNwH4R8UaybghYs+uhUdvmxrdL+6z0OPCTig+eTSPio9U2GBG/jYi/iYgtI2I/sqC6t1YFJO1N1pL9UETc1sK25gA7NnBCbwHwElnr9LVdFlx3LZI2JesyuQQ4W9LYGkXHkh3rr0TEqohYTtYV08wHyNNkrfoJuXnbVilX6/caWL/yfVBpN7ITsLYOcpCXLF369nngeLK+7k9LmpwWb0bWinwmhURlf3czTpM0RtK2wCeBq6uU+QGwi6S/l7RRGt4mabcav8NukjaTNFLS8cBhwH/WKPsWstbpJyLi+61sKyKWAIuAfXPrbyBpNLBRNqnRkkam8s+n3/fTaR8TyC5f/EFu/dFkH6oAo9J0Lf8F9EfEiWTdNV+rUc9lZFfifFTShpK2IOvXnpPb78i0LwEbpXqv9fcWEa+QtezPTv+x7Qp8cJA61lt/EhV97KkPfywwu+h2rbs4yJv3/YrryK9LLclvA+dFxAMRsRD4v8Dl6dKvLwFvILvCZDbNnxzLu4Gsr/Z+svC5pLJARPyJLECnk7XYnyTrBqh1OdrhZF1DK8kuB5waEU8PLEy/71+nyVPJ/r2/JHcs5hXdVhUXkX0ADjiQ7MPvRl4/QXxLbvnJZOcrfk92Yvg7ZJcrDnghLQf4Da93x6xB0jRgKjDwX8opwD6SjqtRz/em8k+Tffi8DHwqt/yWtK/9ya7UeYHX/wOrdDLZZZxPkl16eSVQtVtnkPU3Tet/k+y/g7z/Q3bitpFtWhdRhL9YoltJCmBiRCwa7rqUJX3g/ZrsksInhrs+w0HSecDWEVHt6pVGtzWKrEvlwIh4quXKWUdq5uYCs7ZJrcZql9its1J3ykiyq2XeRnbJ6ollbDsdz13L2JZ1rrpdK5IulfSUpAdrLJek/5a0SNIcSfuUX02zddpmZP3cfybr8/8Psi4z61Ct5KKkGZIWpqHl/7qgQNeKpAPJ+hi/FRFr3cgh6UjgE2Rn7PcD/itdnWBmtk5qNhfTRQ79ZDfxBdn5rbdGxMpW6lO3RR4RPyW7iaOWaWS/TETEbLLrov3kOjNbZ7WQi4cDt0bEihTet5KdNG9JGX3k27DmzQhL0ry1TlRJOons8jA22XiTt7554i4l7N7M1nW/euDXyyKipUfxjt35DfHy868WKvvcEy/NI3vswYC+iOhrYHe1crHW/JYM6cnOdCD6AN46eZ+455afDOXuzaxLbbTVG6vdrdqQl59/lb1PqnUz85p+9q+/ezEielvd51Ap4zrypax5V9mENM/MbH1VKxfbkpdlBPks4IPpLO3bgWfX1+t/zcySWrl4M3BYuht7DNmNeje3urO6XSuSrgQOAnokLSG7rXwjgIj4Gtkdd0eS3d32PPAPrVbKzKyTNZuLEbFC0ufInjsPcE5EDHbStJC6QR4Rx9ZZPvDQezOz9UIruRgRl7LmYyRa5metmJl1OQe5mVmXc5CbmXU5B7mZWZfz0w/NrGO9OHewR9fbAAe5mXUcB3hjHORm1jEc4M1xH7mZdQSHePPcIjezYeUAb51b5GY2bBzi5XCL3MyGnAO8XG6Rm9mQcoiXz0FuZkPGId4eDnIzGxIO8fZxkJtZ27US4k9ef1/9Qus5B7mZtZVDvP0c5GbWNg7xoeEgN7O2cIgPHQe5mXWUbghxSVMlLZC0SNLpVZZfIOn+NDws6Zncsldyy2aVUR/fEGRmpVuXr1CRNAK4EDgUWALcJ2lWRMwfKBMRn8qV/wSwd24TL0TE5DLr5Ba5mZVqPehS2RdYFBGLI+Il4Cpg2iDljwWubGeFHORm1hG6JMQBtgEez00vSfPWIml7YEfg9tzs0ZL6Jc2WdHQZFXLXipmVppO7VDbdZCT777tdobI/43c9kvpzs/oioq+J3U4Hro2IV3Lzto+IpZJ2Am6XNDciHmli269xkJtZKdaxLpVlEdFbY9lSYNvc9IQ0r5rpwMfzMyJiafq5WNKdZP3nLQW5u1bMzBpzHzBR0o6SRpKF9VpXn0jaFRgD/CI3b4ykUWm8BzgAmF+5bqPcIjezlq1jrfFBRcRqSScDNwMjgEsjYp6kc4D+iBgI9enAVRERudV3Ay6S9CpZQ/rc/NUuzXKQm1lHm3/Fz4e7CmuJiBuBGyvmnVUxfXaV9e4G9ii7Pu5aMbOWrE+t8U7lIDezjtWJrfFO5CA3s6a5Nd4ZHORmZl3OQW5mHcndKsU5yM2sKe5W6RwOcjOzLucgN7OO426VxhQK8gIPUd9O0h2Sfi1pjqQjy6+qma0L3K1SvrpBnnuI+hHAJOBYSZMqiv0/4JqI2JvsttSvll1RM+scnfyUw/VRkRZ5kYeoB/DGNL458Pvyqmhm6xN3qzSuSJAXeYj62cDxkpaQPX/gE9U2JOmk9ED1/mXLlzVRXTMzq1TWyc5jgW9GxATgSOBySWttOyL6IqI3Inp7xvWUtGsz6xbuH2+PIkFe5CHqJwDXAETEL4DRgJPabB3k/vHOUyTIizxE/XfAIQCSdiMLcr/aZmZDoG6QR8RqYOAh6g+RXZ0yT9I5ko5KxU4FPizpAbJvi55Z8TB1MzNrk0JfLFHvIerpGy4OKLdqZmZWhO/sNLOO4UsPm+MgNzPrcg5yM7MGFXhsyUxJT0u6Pw0n5pbNkLQwDTPKqI+/fNnMrAG5x5YcSnaD5H2SZqVzhXlXR8TJFeuOBT4L9JLdEf/LtO7KVurkFrmZWWOKPLaklsOBWyNiRQrvW4GprVbILXIzWy+88Q1bcNiexxQqex539Ujqz83qi4i+NF7tsSX7VdnM+yQdCDwMfCoiHq+xbuUjTxrmIDczW9uyiOhtYf3vA1dGxCpJ/whcBhxcTtXW5q4VM7PG1H1sSUQsj4hVafJi4K1F122Gg9zMrDF1H1siaXxu8iiyu+Ihu0P+MEljJI0BDkvzWuKuFTOzBkTEakkDjy0ZAVw68NgSoD8iZgH/lB5hshpYAcxM666Q9DmyDwOAcyJiRat1cpCbWceYdNwBXXF3Z4HHlpwBnFFj3UuBS8usj7tWzMy6nIPczKzLOcjNrCGj99hyuKtgFRzkZjZktj76bcNdhXWSg9zMrMs5yM2so0w6zt9R0ygHuZlZl3OQm1nDfMKzszjIzWxI+YRn+RzkZtZx3E/eGAe5mTXF3Sudw0FuZkPO3SvlcpCbWUdy90pxDnIza5q7VzqDg9zMhkWR7hW3yotxkJtZS9wqH34OcjMbNj7pWQ4HuZm1rJ2tcnev1OcgN7Nh1Y2tcklTJS2QtEjS6VWWnyJpvqQ5km6TtH1u2SuS7k/DrMp1m+EgN7NSrC995ZJGABcCRwCTgGMlTaoo9mugNyL2BK4Fvphb9kJETE7DUWXUyUFuZsOuy1rl+wKLImJxRLwEXAVMyxeIiDsi4vk0ORuY0M4KbdjOjZvZ+mX0Hlvy4tynm1p366PfxpPX31dyjV73hg23YK9xxxQsfWqPpP7cjL6I6Evj2wCP55YtAfYbZGMnADflpkenba8Gzo2I6wtWqiYHuZmVqpUw7yDLIqK31Y1IOh7oBd6Zm719RCyVtBNwu6S5EfFIK/sp1LVSr2M/lflA6tyfJ+k7rVTKzNZPXdLFshTYNjc9Ic1bg6QpwJnAURGxamB+RCxNPxcDdwJ7t1qhukFepGNf0kTgDOCAiNgd+OdWK2Zm3auVE59dEOb3ARMl7ShpJDAdWOPqE0l7AxeRhfhTufljJI1K4z3AAcD8VitUpEVet2Mf+DBwYUSsBMhX3MzWT+tqmEfEauBk4GbgIeCaiJgn6RxJA1ehnA9sCny34jLD3YB+SQ8Ad5D1kbcc5EX6yIt07O8CIOnnwAjg7Ij4UeWGJJ0EnASw3YRtKxeb2TpmHekvX0tE3AjcWDHvrNz4lBrr3Q3sUXZ9yrr8cENgInAQcCzwdUlbVBaKiL6I6I2I3p5xPSXt2sw6WbMt805ulXeaIkFepGN/CTArIl6OiEeBh8mC3czMYd5mRYK8bsc+cD1Za3ygA38XYHF51TSzbucwb5+6QV6wY/9mYLmk+WQd+KdFxPJ2VdrMutP6chv/UCt0Q1CBjv0ATkmDmVlN6+oJ0OHkZ62Y2ZBzy7xcDnIzGxaj99jSgV4SB7mZDSuHeesc5GY27Nw6b42D3Mw6hsO8OQ5yM+sobp03zs8jN7OO5DAvzi1yM7Mu5yA3M+tyDnIzsy7nIDcz63IOcjOzLucgNzPrcg5yM7MGSZoqaYGkRZJOr7J8lKSr0/J7JO2QW3ZGmr9A0uFl1MdBbmbWAEkjgAuBI4BJwLGSJlUUOwFYGRE7AxcA56V1J5F9Oc/uwFTgq2l7LXGQm5k1Zl9gUUQsjoiXgKuAaRVlpgGXpfFrgUMkKc2/KiJWpa/FXJS215JhvbPz2c2WDefuzcxq6ZHUn5vui4i+NL4N8Hhu2RJgv4r1XysTEaslPQuMS/NnV6y7TauVHbYgf2WDl4dr12a2Hhrx6kZs/qeeosWXRURvO+tTJnetmJk1ZimwbW56QppXtYykDYHNgeUF122Yg9zMrDH3ARMl7ShpJNnJy1kVZWYBM9L4+4Hb03cbzwKmp6tadgQmAve2WiE//dDMrAGpz/tk4GZgBHBpRMyTdA7QHxGzgEuAyyUtAlaQhT2p3DXAfGA18PGIeKXVOjnIzcwaFBE3AjdWzDsrN/4i8Lc11v0C8IUy6+OuFTOzLucgNzPrcg5yM+toDyy/brir0PHcR25mHckBXpyD3Mw6igO8cQ5yM+sIDvDmOcjNbFg5wFvnk51mNmwc4uVwkJvZsHCIl8ddK2Y2pBzg5XOL3MyGjEO8PRzkZjYkHOLt4yA3s7ZrJcRvmeMPgHoKBXm9b4zOlXufpJDUNd+sYWbt5RBvv7pBXvAbo5G0GfBJ4J6yK2lm3ckhPjSKtMiLfGM0wOeA84AXS6yfmXUph/jQKRLk1b4xeo1vfZa0D7BtRPxwsA1JOklSv6T+5ctWNFxZM+sODvGh1fLJTkkbAP8JnFqvbET0RURvRPSO6xnb6q7NbB3jEG9OkSCv963PmwFvAe6U9BjwdmCWT3iarZ/W98sMJY2VdKukhennmCplJkv6haR5kuZI+rvcsm9KelTS/WmYXG+fRYJ80G+MjohnI6InInaIiB2A2cBREdFfYNtmtg5xlwoApwO3RcRE4LY0Xel54IMRsTswFfiSpC1yy0+LiMlpuL/eDusGeUSsBga+Mfoh4JqBb4yWdFS99c1s/eAQf8004LI0fhlwdGWBiHg4Iham8d8DTwFbNrvDQs9aqfeN0RXzD2q2MmZmHaJHUr5XoS8i+gquu1VEPJHGnwS2GqywpH2BkcAjudlfkHQWqUUfEasG24YfmmVmLeuG1ni88DIvzn26aPFlEVHzPJ+kHwNbV1l05hr7jAhJMch2xgOXAzMi4tU0+wyyD4CRQB/wGeCcwSrrIDezjnb3vb8b7iqsJSKm1Fom6Q+SxkfEEymon6pR7o3AD4EzI2J2btsDrflVkr4B/Eu9+vhZK2bWkm5ojQ+xWcCMND4DuKGyQLpw5DrgWxFxbcWy8emnyPrXH6y3Qwe5mXWsTmyNF3AucKikhcCUNI2kXkkXpzIfAA4EZla5zPAKSXOBuUAP8Pl6O3TXipk1za3xtUXEcuCQKvP7gRPT+LeBb9dY/+BG9+kWuZl1pC5tjQ8LB7mZNcWt8c7hIDcz63IOcjPrOO5WaYyD3Mwa5m6VzuIgNzPrcg5yM+so7lZpnIPczIaMu1Xaw0FuZg1Z3784ohM5yM3MupyD3Mw6hvvHm+MgN7Mh4f7x9nGQm5l1OQe5mRXmE52dyUFuZtblHORmZl3OQW5mHcFXrDTPQW5m1uUc5GZmJZI0VtKtkhamn2NqlHsl932ds3Lzd5R0j6RFkq5OX9Q8KAe5mVm5Tgdui4iJwG1pupoXImJyGo7KzT8PuCAidgZWAifU26GD3MysXNOAy9L4ZcDRRVeUJOBg4NpG1neQm5mtrUdSf244qYF1t4qIJ9L4k8BWNcqNTtueLenoNG8c8ExErE7TS4Bt6u1wwwYqZ2bWtV5e+TxPXn9f0eLLIqK31kJJPwa2rrLozPxERISkqLGZ7SNiqaSdgNslzQWeLVrBPAe5mVmDImJKrWWS/iBpfEQ8IWk88FSNbSxNPxdLuhPYG/gesIWkDVOrfAKwtF593LViZlauWcCMND4DuKGygKQxkkal8R7gAGB+RARwB/D+wdav5CA3MyvXucChkhYCU9I0knolXZzK7Ab0S3qALLjPjYj5adlngFMkLSLrM7+k3g7dtWJmVqKIWA4cUmV+P3BiGr8b2KPG+ouBfRvZp1vkZtYR9t93u+GuQtdykJuZdblCQS5pqqQF6ZbRte5SknSKpPmS5ki6TdL25VfVzMyqqRvkkkYAFwJHAJOAYyVNqij2a6A3IvYkuyPpi2VX1MyG317jjhnuKlgVRVrk+wKLImJxRLwEXEV2C+prIuKOiHg+Tc4mu/bRzMyGQJEg3wZ4PDdd75bRE4Cbqi2QdNLALa/Ll60oXksz63qH7enWfLuUerJT0vFAL3B+teUR0RcRvRHRO65nbJm7NrN1gK9caU6R68iXAtvmpqveMippCtlzBt4ZEavKqZ6ZmdVTpEV+HzAxPex8JDCd7BbU10jaG7gIOCoiqj5XwMzWDT7h2XnqBnl6cMvJwM3AQ8A1ETFP0jmSBh6Gfj6wKfDdym+7MDMb4H7y9ih0i35E3AjcWDHvrNx4zSeBmZk1Yv99t/MXMTfId3aamXU5B7mZNayVfnJ3r5TPQW5mHceXITbGQW5mTfHVK53DQW5mQ87dK+VykJtZR3L3SnEOcjNrmk96dgYHuZl1rG5slUsaK+lWSQvTzzFVyrwr3Tw5MLwo6ei07JuSHs0tm1xvnw5yM2uJW+VrOR24LSImArel6TWkR39PjojJwMHA88AtuSKnDSyPiPvr7dBBbmYdrQtb5dOAy9L4ZcDRdcq/H7gp950ODXOQm1nL3Cpfw1YR8UQafxLYqk756cCVFfO+kL468wJJo+rtsNCzVszM6tlr3DE8sPy64a5GTS+seI75V/y8aPEeSf256b6I6BuYkPRjYOsq652Zn4iIkBS1diJpPLAH2UMJB5xB9gEwEugDPgOcM1hlHeRmNuwO2/MYbpnTUR8CyyKit9bCwR4UKOkPksZHxBMpqAd7tPcHgOsi4uXctgda86skfQP4l3qVddeKmZXGXSxA9n0NM9L4DOCGQcoeS0W3Sgp/JImsf/3Bejt0kJtZqXzrPucCh0paCExJ00jqlXTxQCFJO5B9+9pPKta/QtJcYC7QA3y+3g7dtWJmHaMDu1gaFhHLgUOqzO8HTsxNP0aVL7KPiIMb3adb5GZWOnexDC0HuZm1hcN86DjIzaxtHOZDw0FuZm3lMG8/B7mZtZ3DvL0c5GY2JHxZYvs4yM1syDjM28NBbmZDaq9xxzjQS+YgN7Nh4TAvj4PczIaNW+flcJCb2bBzmLfGQW5mHcGt8+Y5yM2sozjQG+enH5pZR1ozzE8dtnp0A7fIzcy6nIPczKzLOcjNzLqcg9zMrMs5yM3MulyhIJc0VdICSYsknV5l+ShJV6fl96QvFTUzW+9I+ltJ8yS9Kql3kHJVc1XSjilHF6VcHVlvn3WDXNII4ELgCGAScKykSRXFTgBWRsTOwAXAefW2a2a2jnoQeC/w01oF6uTqecAFKU9XkuXroIq0yPcFFkXE4oh4CbgKmFZRZhpwWRq/FjhEkgbb6IhXNyqwazOz7hIRD0XEgjrFquZqys2DyXIUslw9ut4+i9wQtA3weG56CbBfrTIRsVrSs8A4YFm+kKSTgJPS5KqejXd6sMD+262Hinqup3WAzqiH6/C6TqhHJ9QB4M2tbuCRV566edrKL/cULD5aUn9uui8i+lqtQ06tXB0HPBMRq3Pzt6m3sSG9szMdiD4ASf0RUbP/aKh0Qj06oQ6dUg/XobPq0Ql1GKhHq9uIiKll1AVA0o+BrassOjMibihrP0UVCfKlwLa56QlpXrUySyRtCGwOLC+lhmZmHSYiprS4iVq5uhzYQtKGqVVeLW/XUqSP/D5gYjqTOhKYDsyqKDMLmJHG3w/cHhFRYNtmZuujqrmacvMOshyFLFfrtvDrBnn6VDgZuBl4CLgmIuZJOkfSUanYJcA4SYuAU4C1LlGsosz+plZ0Qj06oQ7QGfVwHV7XCfXohDpA59SjLknHSFoC/BXwQ0k3p/lvknQj1M7VtInPAKekPB1Hlq+D79MNZzOz7uY7O83MupyD3Mysy7UlyFu5pV/SGWn+AkmHt7EOp0iaL2mOpNskbZ9b9oqk+9NQeWK37HrMlPR0bn8n5pbNkLQwDTMq1y2xDhfk9v+wpGdyy0o5FpIulfSUpKr3Dijz36mOcyTtk1tW1nGoV4fj0r7nSrpb0l65ZY+l+fe3eilcgXocJOnZ3HE/K7ds0NeyxDqcltv/g+l9MDYtK/NYbCvpjvS3OE/SJ6uUaft7o+tFRKkDMAJ4BNgJGAk8AEyqKPMx4GtpfDpwdRqflMqPAnZM2xnRpjq8C9g4jX90oA5p+rkhPBYzga9UWXcssDj9HJPGx7SjDhXlPwFc2oZjcSCwD/BgjeVHAjcBAt4O3FPmcShYh/0Htk126/Q9uWWPAT1DdCwOAn7Q6mvZSh0qyr6H7Eq0dhyL8cA+aXwz4OEqfyNtf290+9COFnkrt/RPA66KiFUR8SiwKG2v9DpExB0R8XyanE12vWbZihyLWg4Hbo2IFRGxErgVaOaGhkbrcCxwZRP7GVRE/BRYMUiRacC3IjOb7Fra8ZR3HOrWISLuTvuA9r0nihyLWlp5P7VSh7a8J1I9noiIX6XxP5FdwVF5J2Pb3xvdrh1BXu3W08oXZo1b+oGBW/qLrFtWHfJOIPvEHzBaUr+k2ZKObmL/jdbjfelfxmslDdwkMOTHInUv7Qjcnptd1rGop1Y9yzoOjap8TwRwi6RfKnvURLv9laQHJN0kafc0b8iPhaSNycLxe7nZbTkWyrpY9wbuqVjUae+NjrPef/mypOOBXuCdudnbR8RSSTsBt0uaGxGPtKkK3weujIhVkv6R7D+Vg9u0r3qmA9dGxCu5eUN5LDqCpHeRBfk7crPfkY7DXwC3SvpNatW2w6/Ijvtzko4Ergcmtmlf9bwH+HlE5FvvpR8LSZuSfVj8c0T8sZVtrY/a0SJv5JZ+tOYt/UXWLasOSJoCnAkcFRGrBuZHxNL0czFwJ1kroRl16xERy3P7vhh4ayO/Qxl1yJlOxb/QJR6LemrVs6zjUIikPcleh2kR8dpjJnLH4SngOprr8iskIv4YEc+l8RuBjST1MMTHIhnsPVHKsZC0EVmIXxER/1ulSEe8Nzpa2Z3uZK38xWT/og+ckNm9oszHWfNk5zVpfHfWPNm5mOZOdhapw95kJ44mVswfA4xK4z3AQpo/oVSkHuNz48cAs+P1EzmPpvqMSeNj21GHVG5XspNYasexSNvYgdon+N7Nmie07i3zOBSsw3Zk52X2r5i/CbBZbvxuYGqLfyeD1WPrgdeBLCR/l45LodeyjDqk5ZuT9aNv0q5jkX6vbwFfGqTMkLw3unloz0azs8wPkwXlmWneOWQtX4DRwHfTH829wE65dc9M6y0AjmhjHX4M/AG4Pw2z0vz9gbnpj2QucEKbj8W/A/PS/u4Ads2t+6F0jBYB/9CuOqTps4FzK9Yr7ViQteqeAF4m68s8AfgI8JG0XGQP2n8k7au3DcehXh0uJnuQ/8B7oj/N3ykdgwfSa3Vmi++JevU4OfeemE3ug6Xaa9mOOqQyM8kuPsivV/axeAdZn/uc3HE/cqjfG90++BZ9M7Mu5zs7zcy6nIPczKzLOcjNzLqcg9zMrMs5yM3MupyD3MysyznIzcy63P8HXp+ihzxIcacAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Test problem on a fine grid\n",
    "test = Example929(161,81)\n",
    "print(test)\n",
    "\n",
    "# run the solver\n",
    "start = time.process_time()\n",
    "itt, err = SOR(test,tol=0.5e-14)\n",
    "stop = time.process_time()\n",
    "print(\"Solver took {:.3g} seconds.\".format(stop-start))\n",
    "print('Converged after {} itterations, final residual is {}'.format(itt,err))\n",
    "\n",
    "# plot the solution\n",
    "fig, ax1 = plt.subplots()\n",
    "cmap = plt.get_cmap('PiYG')\n",
    "cf = ax1.contourf(test.x,test.y,test.u,cmap=cmap)\n",
    "fig.colorbar(cf, ax=ax1)\n",
    "ax1.set_title(f'Example 9.29 ({test.Ni} x {test.Nj} grid)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "roman-planner",
   "metadata": {},
   "source": [
    "## A new test problem\n",
    "\n",
    "Now we know we haven't changed the solution to our original test problem we need one that includes Neumann and Dirichlet boundary conditions.\n",
    "\n",
    "### Porous Flow\n",
    "\n",
    "Water seeping through a porous medium can be modelled using the Laplace equation.  If we have a \"box\" _z<sub>0</sub>_ meters tall and _s_ meters wide with a linear pressure gradient _h = z<sub>0</sub> + cx_ applied at the upper boundary, an impermeable base and constant flow boundaries on the left and right.  The pressure in the medium can be modlled as\n",
    "\n",
    "$$h_{zz}+h_{xx}=0$$ with boundary conditions\n",
    "\n",
    "$$h_z(x,0)=0,\\,h_x(0,z)=0\\text{ and }h_x(s,z)=0$$\n",
    "and\n",
    "$$h(x,z_0)=z_0+c x.$$\n",
    "\n",
    "This test problem has an anlytical solution developed by Toth (1962):\n",
    "$$h(x,z) = z_0+\\frac{cs}{2}-\\frac{4cs}{\\pi^2}\\sum_{m=0}^\\infty \n",
    "\\frac{\\cos[(2m+1)\\pi x/s]\\cosh[(2m+1)\\pi z/s]}\n",
    "{(2m+1)^2\\cosh[(2m+1)\\pi z_0/s]}.$$\n",
    "\n",
    "Let us consider a 200 wide by 100 m deep medium with _c_ = 0.176."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "steady-hungary",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Uniform 101x51 grid from (0.0, 0.0) to (200.0, 100.0).\n",
      "Boundaries conditions are - left: Neumann, right: Neumann, top: Dirichlet, bottom: Neumann.\n",
      "Solver took 3.64 seconds.\n",
      "Converged after 1737 itterations, final residual is 4.9461867869058823e-11\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAAEICAYAAACd0wWxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAm0klEQVR4nO3de7xcZX3v8c/XcAsQCRjK4SZBQSoSLBABj6KpeEkRAa8FbQXEQ2m91vYgFCucg5yXSrXWQ9VGwWDlIJSC0IJKpNBYJVgukSRc5C6BSAjhnhQI+Z0/1jO6Msxlz15zWWvW9/16zWvPPLMuz1579m9+81vPs0YRgZmZldOLRt0BMzNrz0HazKzEHKTNzErMQdrMrMQcpM3MSsxB2sysxByk7TckzZf0uXT/IEm3T3I7p0n6bn9713Zf10j6cIfn95R0vSQNoz9VJOkbkv66w/Mhabd0/0uS/nR4vbPaB2lJ90paK+kpSQ+lQLXlqPs1ahHxk4jYY9T96IPTgb+JNCFA0kdT0H5G0vzmhSUdLOk2SWskXS1pl9xz75P0s/TcNf3qoKSZKRA+lbv9de75gey3ISJOiIjTJ7j43wB/JWmTfvfDWqt9kE7eERFbAvsCs4HP9LKypI0G0isrRNL2wO8D3881Pwh8DjinxfIzgIuBvwa2Aa4HLsgtshr4CvD5gXQYpkfElumWD5oD26+kKb0sHxErgNuAw/rdF2vNQTonIh4AfgDsBSDpMEnLJD2WPla/srFsysA/Lelm4GlJG3VZ/jcfGdPjfGlhhqR/TeutlvQTSS3/Nmk7fybpDklPSjpd0stTpvWEpAvzWY6kQyUtTtv+maS9c8/tI+nGtJ0LgM1yz82RtDz3+CRJd6Vlb5H0zi6HcxNJ30nLL5M0O7etHST9s6SHJd0j6eO55/aXdG3q7wpJZzX9Pm9Jme7jks4COpUx3gLcGBH/1WiIiIsj4vvAIy2WfxewLCL+Ka1zGvBqSb+b1v1xRFxIFug7Sq+N6xpv4JL+NB2Hzbqt26yX/aZ9nZiO3YOSPtxUrpgv6euSrpD0NPD7+ddiWuZ/5tb/UItdXAO8vdffwybHQTpH0s7AIcBNkl4BnA98EtgWuAL4l6aPeUeRvVinAy+bwPLt/AWwPK23HfBXQKf5+m8D9gMOBE4E5gF/BOxM9gZzVPp99iHLGP8EeAnwD8BlkjZN/fo+8I9kWeM/Ae/usM+7gIOArYD/BXw3ZartHAZ8j+zYXAaclfr0IuBfgF8AOwIHA5+U9La03vPAnwMzgNem5/8srdvIdD+Tnr8LeF2HPswCeqmrvyr1C4CIeDrt41U9bKPhTOAZ4DOSdgf+D/BH+TeMFu6TtFzSt9Pv2jNJc4FPAW8GdgPmtFjs/cAZwDTgP1qs/5dkb3C7p+00uxV49WT6Z71zkM58X9JjZC/Yfyf7h/pD4PKIWBARz5HV4qYC/z233lcj4v6IWDvB5dt5Dtge2CUinkv14E5B+osR8URELAOWAldGxN0R8TjZJ4F90nLHA/8QEddFxPMRcS5Z4Dgw3TYGvpL2eRHwn+12mLLLByNifURcANwB7N+hj/8REVdExPNkbwSNf+rXANtGxP+OiGcj4m7gm8CRaT83RMSiiFgXEfeSvbG8Ma17CFmme1E6xl8Bft2hD9OBJzs832xL4PGmtsfJgllPImI98EHg42RvUl+MiJvaLL6K7LjsQvbmOw04r9d9Ju8Dvh0RyyJiDdmngWaXRsRP09+y+U2jsf7S9CbVav0nyY6tDYFrqZkjIuLH+QZJOwD3NR5HxHpJ95Nlfw335+5PZPl2ziT7Z7hS2SCEeRHRqf74UO7+2haP/1u6vwtwtKSP5Z7fJPU1gAea3gzuow1JHyTL0Gampi3Jstl28sFzDbBZ+ui/C7BDelNsmAL8JO3nFcCXyc4NbE72Gr0hLbcDuWMeEZGOcTuP0luAfQp4cVPbi+kt0P9GRNwr6WqyN5e/77DcU2T1b4CHJH0UWCFpWkT0uu8dctuCDV+jndry69+Qe9zqNTENeKzHftkkOZNu70GygAKAsui5M/BAbpnoYfk1ZEGnoRFIiYgnI+IvIuJlZGWCT0k6uA+/w/3AGRExPXfbPCLOB1YAO6Z+Nry01UaUjXD4JvBR4CURMZ0sg5/MsLb7gXua+jQtIg5Jz3+d7MTU7hHxYrLST2M/K8iOaaNfyj9u4WbgFT30bRm5j/GStgBentp7JuntZCWbq8jeiCeq8bqazP/nCmCn3ONWx6fTp7QNjjGtXxOvJFcWssFykG7vQuDtyoZkbUxWN34G+Nkkl18MvF/SlFT3a3yEb5zc2y0FncfJ6rLr+/A7fBM4QdIBymwh6e2SpgHXAuuAj0vaWNK7aF++2ILsH/vh1N9jSSdXJ+HnwJPpxNrUdDz2kvSa9Pw04AngqXTCLj8m93LgVZLelbLyj5N7s2thAbBv/mSdshO8m5Fl71MkNTJ8gEuAvSS9Oy3zWeDmiLgtrTsltW8EvCitu3GrHaea8reADwNHA++QdEibZQ+QtIekF0l6CfBV4JpUvuppv2Svw2MlvVLS5mQjVXpxIXCMsvHlmwOntljmjWRlNRsCB+k2IuJ2spNx/5esZvgOsqF6z05y+U+ktseAD7DhsLDdgR+Tfdy+FvhaRFzdh9/heuB/kJ20exS4EzgmPfcs2WiGY8iGeP0h2Um5Vtu5BfhS6ttDZCfkfjrJPj0PHAr8HnAP2bH6FtkJSchOWr2frMTwTXJD4CJiFfBesqFoj5Adt7b9iIiHgH8DDs81f4asJHQS2d9rbWojIh4mO3l6BtnxOoBUK0/+OC3/dbKTqGtTH1uZR1b7vSIiHgGOA76VgnCzlwE/TL/zUrI396Mms9+I+AFZkL+a7O+9KD31TJt+tlr/K2TH7c708zfSyeI92fD1awMkX/TfxpmkPYFzgf27nIwdS8qGgS4FNo2IdX3Y3peAuyLia4U7ZxPiIG02ZpSNYb+C7BzIucD6iDhipJ2ySeta7pB0jqSVkpbm2raRtEDZhIoFkrZO7ZL0VUl3SrpZ0r6D7LyZtfQnwEqyMd7Ps2Fd3/qsTYw8PcXAxZKuTKPFGs/NSe3LJP171+13y6QlvYGsVvqdiGjMxPsisDoiPi/pJGDriPh0OjHyMbIhRwcAfxcRB0zi9zYzq4Q2MfLFEfFEuv9xYM+IOEHSdLLBBHMj4leSficiVnbaftdMOiIWkp1Yyjuc7GMU6ecRufbvRGYRMF2dZ6WZmVVaqxjZCNBJY3QUZCfFL46IX6XlOgZomPxklu0iu9AKZJMWtkv3d2TDgfLLU9sKmkg6nmxGHLDpftmEu4nqeQLYhl5U4AJeWxTbdaGuF9y3tnh+0utOnbK20L43Z02h9acy+f2Pct9T1hQ85/NUgXWL/dqF1n+22Mul4+yhe2FVRGxbZPuzpJjoob03Gyefn5k5LyLmdVtP0hlks04fJ7vQF2Tj9jdWdjXDaWTVhu902k7hGYdp1lfPr8T0S84DkHaN1sMxW5nT6642NG1msfWLFG8OKrZrDiy2+mYHPjrpdfd68ZJC+55FsfX3od2M6sHvu8j6W99QMFot6r5IR9d3X6StG7ov0sm9xQ4717RpP7bDzNiJeoqJR5xj4b8iYnb3JTcUEacAp0g6mWwi2KlkMXc/smvSTAWulbQoIn7ZbjuTHSf9UKOMkX42UvYH2HC20k5sOEOvoDnFVh9lgLZaenS/qaPuwuTtV2z1mbOKrT+n2Oplch6/vXjZcuBHEfF0Gve/kC4Xq5pskL6MbBYV6eelufYPplEeBwKP58oiBcxh5AG64qYWyKJthAp+eqq6OVQzWCu78mHD4WSXOoAsVr4+zXzdnCz1u7XTtrqWOySdT3acZii7vvCpZDO+LpR0HNlHj/elxa8gG9lxJ1lF69gJ/k4dzCm+iX4omkWPuNRhNTWbYiWP/ShU9pg5q3jZA0oTBVpqEyMPkbQH2eUd7gNOAIiIWyX9kOy6MuuBb0XE0pYbTroG6Yg4qs1TL7gAUJrR9ZFu25y4Of3ZTM2zaLNR6legLqs2MfLsDsufSQ8X3CrhtTvm0NcPOf0I0K5FV9oSChZHq67nU15NCtamoXh9us5KEqSnMZDqU1kyaJc6Kq1okC988rAMf38H6pEpSZAusTHIon3S0MrCgbp34xukxyWLLoFZBcdIWwkULXlAX7JpcKDu1XgG6X4F6DHIom1MlKHkAQ7UIzB+QXrcAnRZ/jmtkFJMaulHNt1HM2c5WE/EeAXpspQ4Gsag1NEPRadlW8n0KZtucKDubHy+LbyfAbosWXQf+KRhfyxh1ujfbA6k+LU8ik5uaSg4yaVZI1APazx1YzxZFYxHJl22DBr6k0W71GE14xLIC1U7kx5EcB6jLNqspZJm03n5QD3OsxUnorpBuswB2rVoa+HR/aYWv3RpP0oe/TTAQN3wgsy6ZkG7euWOaTPLWd4YUx4jPab6OdKjzycSbUPVCtKDDM5lK3P0oR7tk4a/VfvrdwyaA/XAVCNIDzp77meAdqnDBq1fJ5RLNm7aWitvkG4E5kGXNsqWQVtp9SMbL8WklkFxNj0Q5QjSL9pkw6Bc1Zpzv7JoD72zbsqaTdcwUEs6R9JKSUtzbadLulnSYklXStohtX8gtS+R9DNJHb86C8oSpEfFWbRZ/9UvUM8H5ja1nRkRe0fE7wH/Cnw2td8DvDEiZgGnk76Mu5P6Bul+B+iS1aJ90tAmbBC16RoF6ohYCKxuansi93ALIFL7zyKi8c+5iOzLujuqZ5AucwbtUsfY61tdup+vFQfqvpN0hqT7gQ/w20w67zjgB922U93JLJM1iABdsiy6XzxG2gprBOoBT3jp1SZTYeZuE1x4CTMk5edozouIrmWKiDgFOEXSycBHyb6gFgBJv08WpF/fbTv1zKTNJmmsx1sPckhetbPqVRExO3frGqCbnAe8u/FA0t7At4DDI+KRbivXK0iXPYvu08dX16PLr5Qlj0GrdqDuiaTdcw8PB25L7S8FLgb+OCJ+OZFt1afcUeY6tFlZ9OviS+2UtPxRhKTzya58OkPScrKyxiGS9gDWA/cBJ6TFPwu8BPiaJIB1EdHxM0w9gvSgAvSY1qKtYvp90aVBB2oYq2AdEUe1aD67zbIfBj7cy/bHv9xRlQy6Sh9bezDyC+XnjHU9uapqVAKZrPEN0gcw2ADtLNoK6usU8X6/yQ/zuh775W72AuMZpAedPZc4QPfrpKGH37VXm4x8FBdg2g8H7SbjV5OuSnkjb0xLHWaFtQrUNcsfxieTHnR5o6HEWbTV3CDe7H0505EbjyBdxey5wVl0rVXi0qUO1CNV7SA9rOy5oeRZtCex2MDe9B2oR6aaNelRZM4lD9A2XEuYVarhhUMxjPHTw7I5Ez8xOeI/c7Uy6WFnzoPmUocxgJKHX1djpfyZdBmCcs2yaA+/s7bGKZuuiPIE6TIE41YGFaD7nO24Hm0b6PdU8TwH6qEqR7lji1F3wOqirBNRKjHKI88nEoemUJCW9OeSlklaKul8SZtJ2lXSdZLulHSBpE361dmhq1mZw8bMoGvTDtRDMekgLWlH4OPA7IjYC5gCHAl8AfjbiNgNeJTs2weqZ5AB2id2xkJZs/KhcqAeuKLljo2AqZI2IhvUsgJ4E3BRev5c4IiC+7AuXI+2toaREDhQD9Skg3REPAD8DfArsuD8ONnVYR+LiHVpseXAjq3Wl3S8pOslXc9zD0+2G4PhLNpGoHJ16TwH6oEpUu7YmuxrYXYFdiA7/Td3outHxLzGd4ax8baT7Ub/1bwO3c/hd7Wb7FFWw0oMahqoJZ0jaaWkpbm2MyXdJulmSZdImp7aN5Z0rqQlkm5NX1LbUZFyx5uBeyLi4Yh4jux7u14HTE/lD4CdgAcK7GO4ah6gzQqbTR2D9XxemKAuAPaKiL2BXwKNYPxeYNOImEU25/FPJM3stPEiQfpXwIGSNlf2ZV0HA7cAVwPvScscDVxaYB/DM4wAPYCMxvXo0er3ycOBlTyGXWarUaCOiIXA6qa2K3Nl30VkCStAAFukRHYq8CzwRKftT3oyS0RcJ+ki4EZgHXATMA+4HPiepM+ltpbf9VUqzqDN+q/Mk142Z+JvJPOZISn/m8yLiHk97O1DwAXp/kVkZeIVqRd/HhGr260IBWccRsSpZN+Mm3c3sH+R7Y4lnzC0URvkLMR2GoGwrMF6YlZ1+0bvdiSdQpbEnpea9geeJzuPtzXwE0k/joi7222jHDMOR6nCWbRLHZNT9vHNlR7l0U6Nyh8Nko4BDgU+EBGRmt8P/DAinouIlcBP6XJ06h2khxWgK5JF+8JKNTDK12KNTipKmgucCBwWEWtyT/2KbC4JkrYg+4vc1mlb9Q3SFc6grVzKnpm/wKiThjEL1pLOB64F9pC0XNJxwFnANGCBpMWSvpEW/3tgS0nLgP8Evh0RN3fafnmugjdMwwzQo/6HsEp6dL+pbH3D2lF3Y7DygbrCNeuIOKpFc8sBExHxFNkwvAmrXyY9Jhm069FWSNmSh9mMXYbdL/XJpEcRnMv2j2BWBd0C9fxhdKI86pFJj0n2bPUy8FEeTiIqYfyD9KgC9AD/AQZR6vDIjmIqd/KwwYG69MY7SDuDtoobyzHT1pPxDNIHMdoA7eyk9FfAq2zmOwh+vZba+J04dPZs1rtRTBkfpS2pzJvT+GTSo86eGwb8h6/r0Lub2GfUXRgZlzzqrfpBuizB2Wqv8iWUimSWdVPdckcZA3NFX+Qe2WG/UbeyRwVUK5M+iFpnznUtddiQSx4VTTbGVfkz6aoEZL+wzWwAyhGkp1GdYNyKA7QlS5g1sOGHQ73okssepVGtckeNudRhQ+fkoxQcpIvyC9nGmV/fI+cgXXN1HtlRxSFzIxkz7UDdkaRzJK2UtDTXdqak2yTdLOkSSdOb1nmppKck/WW37TtIFzGkF69LHTZyDtSdzAfmNrUtAPaKiL2BXwInNz3/ZeAHE9m4g/Rk+UVrbQw6Q/cMxHKJiIXA6qa2KyNiXXq4CNip8ZykI4B7gGUT2X45RneYWfmN0YiP5zcXj+632QSXXjtDUv4LvuZFxLwedvch4AIASVsCnwbeAnQtdYCD9OQMMYuuYqmj7FfAswLGKFD3YFVETOqLvSSdAqwDzktNpwF/GxFPSZrQNhyke+Uyh5XASL+otp6BumeSjgEOBQ6OiEjNBwDvkfRFYDqwXtJ/RcRZ7bbjIF1jdR7ZMWiDnNRSCg7UHUmaC5wIvDEi1jTaI+Kg3DKnAU91CtDgE4e9GXIWXcVShw3PyE8g+lMlAJLOB64F9pC0XNJxwFlkc6kXSFos6RuT3b4z6YnyC3IsjX3GO2iN/4saZ9URcVSL5rMnsN5pE9m+M2mzCht5Nt3gJGZgHKQnYgQvQJc6qq+KMxoLcaAeCAfpbsb0heeThjYQY/r/MkoO0p34BWcVUJqSR8OB+H+njxykS8ilDhsLDtR94SDdjl9g1gfDqkuXLptucFZdmIfgteIXlVl/lWyo3lqmsoS9Jrj0zwfal24KZdKSpku6KF039VZJr5W0jaQFku5IP7fuV2eHYsQBehilDp80HE+lzabznFn3rGi54++AH0bE7wKvBm4FTgKuiojdgavS42rwi6eWajdUrgwcrCds0kFa0lbAG0gzayLi2Yh4DDgcODctdi5wRLEuWpV49t4LDfNNoBLZdN6BOGB3UaQmvSvwMPBtSa8GbgA+AWwXESvSMr8Gtmu1sqTjgeMB2OqlBbrRJyV4kXhUh9Vaq//BktSwR6lIuWMjYF/g6xGxD/A0TaWNdHm+aLEuETEvImZHxGy22LZAN/qgBAHarF8ql013cmCLW80UCdLLgeURcV16fBFZ0H5I0vYA6efKYl0csJL80Z1Fm1krkw7SEfFr4H5Je6Smg4FbgMuAo1Pb0cClhXo4SCUJ0MPkkR2jMeyTk2OVTddc0XHSHwPOk7QJcDdwLFngvzBdU/U+4H0F9zEYNQzQZlY9hYbgRcTiVFfeOyKOiIhHI+KRiDg4InaPiDdHxOruWxqykgVolzpsEJxND4ekcyStlLQ013Zmmj9ys6RLJE3PPXeypDsl3S7pbd22X69p4TU98WDdeay0FTAfmNvUtgDYKyL2Bn4JnAwgaU/gSOBVaZ2vSZrSaeP1CdIOzjZio3gjcDY9eBGxEFjd1HZlRKxLDxcBO6X7hwPfi4hnIuIe4E5g/07br8e1O0ocoIdZ6vBJQ7PMGjbv4U3z5zMkXZ9rmBcR83rY3YeAC9L9Hdlw9Pfy1NbW+AfpEgdos2F4dL+pbH3D2lF3o8pWRcTsyawo6RRgHXDeZHc+vkG6AsHZJwzNxpekY4BDgYPTxD6AB4Cdc4vtlNraGs+adAUCtNXTqE5QujY9XJLmAicCh0XEmtxTlwFHStpU0q7A7nS5Fup4ZdIOziPliytZHUk6H5gDzJC0HDiVbDTHpsACSQCLIuKEiFgm6UKyiX/rgI9ExPOdtj8+QbpiAXrYpQ6fNDTXpgcjIo5q0Xx2h+XPAM6Y6ParX+7w2GezCXPZo3qqGaQrfkUsnzAsp2HViz1xxnpRrSBd4cBsVhbOpqul3DVpB2Qzq7lyBOktqE1AHkWpwycNrZlPIlZHtcodZmOiDHVplz2qwUF6iHzC0Mx6VY5yh5mNRF3LHmuZyk3sM+puTIgz6SFxFm1mk+EgPeZ80rA3w6wVl6EuDa5Nl52DtPWFr9tRbQ7U5eUgPQQudZjZZDlIW2VU5URPL8pS8gBn02XlID1gzqKtShyoy8dBeoz5pKFZ9TlID5CzaKsiZ9O9kXSOpJWSluba3itpmaT1kmY3Lb+3pGvT80skbdZp+w7SZiNWprp0gwN1T+YDc5valgLvAhbmGyVtBHwXOCEiXkX2jS7Pddq4g/SAOIuurjIGzVFwoJ6YiFgIrG5quzUibm+x+FuBmyPiF2m5R7p9fZaDtJlZZzMkXZ+7HV9gW68AQtKPJN0o6cRuK/jaHWPKJw2tH8b12h5r2LyXT0yrImJ298UmZCPg9cBrgDXAVZJuiIir2q3gTHoAXOqwXpW5xOKyR18tBxZGxKqIWANcAezbaQUHaTPryoG6b34EzJK0eTqJ+Ebglk4rOEj3mbNos3qRdD5wLbCHpOWSjpP0TknLgdcCl0v6EUBEPAp8GfhPYDFwY0Rc3mn7rklbYb64Un8sYVapj+W41qeLioij2jx1SZvlv0s2DG9CnEn3UVmyaJ80tEFx2WP4HKTNrCcO1MPlIN0nZcmirT/KPNqiDByoh6dwkJY0RdJNkv41Pd5V0nWS7pR0gaRNinfTrB6q9ObgQD0c/cikPwHcmnv8BeBvI2I34FHguD7so9ScRVtdOVAPXqEgLWkn4O3At9JjAW8CLkqLnAscUWQf1hufNLRhc6AerKKZ9FeAE4H16fFLgMciYl16vBzYsdWKko5vzIXn8YcLdmN0nEWbOVAP0qTHSUs6FFgZETdImtPr+hExD5gHoFfMjsn2w2zclH28dDtVGke99vmpLHmiGvX/IpNZXgccJukQYDPgxcDfAdMlbZSy6Z2AB4p3s5ycRZttqEqBuiomXe6IiJMjYqeImAkcCfxbRHwAuBp4T1rsaODSwr00s8pw6aO/BjFO+tPApyTdSVajPnsA+xi5MmbRPmnYX6McDleloXitOFD3T1+u3RER1wDXpPt3A/v3Y7tmVl0uffSHZxxOQhmzaLMyckZdnIO0FVLFUQg2XI/uN9XBugAH6R45i7ZhqXpdupkD9eQ4SJvZ0IxjVi3pHEkrJS3Ntb1X0jJJ6yXNzrW/RdINkpakn2/qtn0H6R6UOYv2yA6rkjEL1vOBuU1tS4F3AQub2lcB74iIWWRDlP+x28YdpM1sZMYhUEfEQmB1U9utEXF7i2VviogH08NlwFRJm3bavoP0BJU5i7bxNW516VYaWXWJA/aMxnWG0u34Pm333WTfcfhMp4X8HYdmHVT1OhpV1RyoBzXOOp6ewtpFW0908VURMbv7YhMn6VVkl3V+a7dlHaQnwFm02Wi0zq6rPUEmXeL5EuCDEXFXt+Vd7hgDPmloVg2SpgOXAydFxE8nso6DdBfOom3U6lCXrjJJ5wPXAntIWi7pOEnvlLQceC1wuaQfpcU/CuwGfFbS4nT7nU7bd7nDzKyAiDiqzVOXtFj2c8Dnetm+M+kOnEWb2ag5SJuZlZiDdBvOoq1MXJeuLwfpivPIDrPx5iDdgrPoifEkD7PBc5A268KlBhslB+kmzqKtrPxmUU8eJ53jAG1WE08Di0bdiYlxJl1hPmloNv4cpBNn0WZWRg7SZhXiunT9OEjjLNrMystB2sysxGofpJ1Fm1mZ1T5Im1WN69L1UusgXeUs2sPvzMpB0jmSVkpammvbRtICSXekn1un9q0k/YukX0haJunYbtuvbZCucoC24XP2ah3MB+Y2tZ0EXBURuwNXpccAHwFuiYhXA3OAL0napNPGaxukzcz6ISIWAqubmg8Hzk33zwWOaCwOTJMkYMu03rpO269lkHYWbVXnzH6oZki6Pnc7fgLrbBcRK9L9XwPbpftnAa8EHgSWAJ+IiPWdNuRrd5hZ/TwJ/GTCS6+KiNmT3VVEhKRID98GLAbeBLwcWCDpJxHxRLv1a5dJO4s2syF4SNL2AOnnytR+LHBxZO4E7gF+t9OGahekx4FHdpiV3mXA0en+0cCl6f6vgIMBJG0H7AHc3WlDtSp3OIvuH38ri1lG0vlkIzVmSFoOnAp8HrhQ0nHAfcD70uKnA/MlLQEEfDoiVnXafm2CtAO0jZslzPKbZQlExFFtnjq4xbIPAm/tZfuTLndI2lnS1ZJuSYOyP5HaWw7iNjOz3hWpSa8D/iIi9gQOBD4iaU/aD+IeGWfRZlZVkw7SEbEiIm5M958EbgV2pP0gbjMz61FfRndImgnsA1xH+0Hczesc3xgczuMP96MbLTmLNrMqKxykJW0J/DPwyeYB2RERZNMgXyAi5kXE7IiYzVbbFu2G1cRN7DOyfZdxll8Z+2T9VShIS9qYLECfFxEXp+Z2g7iHbhyzaI+RNquXIqM7BJwN3BoRX8491W4Q91CNY4A2s/opMk76dcAfA0skLU5tf0X7QdxmZuXwNNkZtAqYdJCOiP8gmzHTygsGcQ+Ts2gzGxe+doeZWYmNXZB2Fm114xEe422sgrQDtJmNm7EK0mZm42ZsgrSzaDMbR2MTpOvAE1nM6mcsgrSzaDMbFUnnSFopaWmureMlmyW9RtI6Se/ptv2xCNJmZiM0H5jb1Nb2ks2SpgBfAK6cyMYrH6SdRZt5GN4oRcRCYHVTc6dLNn+M7JpHE7quUaWDtAO0mZVUy0s2S9oReCfw9YluqDbfcWhm9hvrn4Un753o0jMkXZ97PC8i5k105YgISY1LNn+F7Mtn12fXqOuuskHaWbSZDcmqiJjd4zoPSdo+IlY0XbJ5NvC9FKBnAIdIWhcR32+3oUqXO8zMSqrlJZsjYteImBkRM4GLgD/rFKChokHaWbSNik/QWTNJ5wPXAntIWp4u0/x54C2S7gDenB5PSuXKHQ7QZlYmEXFUm6c6XrI5Io6ZyPYrmUmbmdVFpYK0s2iz9lyKGU+VCtJWDrPwNUTMhqUyQdpZtJnVUSWCtAO0r4BnVleVCNJmZnVV+iDtLNrM6qxy46TNzIp7Erhm1J2YkFJn0s6izazuShukHaDNzEocpM3MrKRB2lm0mVmmlEHazMwypQvSzqLNzH6rdEHazMx+q1RB2lm0mdmGShWkzcyqRtI5klZKWppr20bSAkl3pJ9bp3ZJ+qqkOyXdLGnfbtsvTZB2Fm1mFTUfmNvUdhJwVUTsDlyVHgP8AbB7uh0PfL3bxksTpM3MqigiFgKrm5oPB85N988Fjsi1fycyi4Dp6dvE2yrHtTvuuGHV2q22uW+EPZgBrBrh/hva9uPnJekHDLUvpf+7DJn7kdml+Cbu/REcO2OCC28m6frc43kRMa/LOttFxIp0/9fAdun+jsD9ueWWp7YVtFGKIB0R245y/5Kuj4jZo+yD++F+uB/DExHN5YlB7iskxWTXd7nDzKz/HmqUMdLPlan9AWDn3HI7pba2HKTNzPrvMuDodP9o4NJc+wfTKI8DgcdzZZGWSlHuKIFu9aVhcT825H5syP0oIUnnA3OAGZKWA6cCnwculHQccB/wvrT4FcAhwJ3AGuDYrtuPmHSpxMzMBszlDjOzEnOQNjMrsVoFaUk7S7pa0i2Slkn6RGo/TdIDkhan2yFD6Mu9kpak/V2f2lpOJR1gH/bI/c6LJT0h6ZPDOh6Dnk5bsB9nSrot7esSSdNT+0xJa3PH5hsD7kfbv4Wkk9PxuF3S2wbcjwtyfbhX0uLUPrDjYUlE1OYGbA/sm+5PA34J7AmcBvzlkPtyLzCjqe2LwEnp/knAF4bYnylkg+53GdbxAN4A7Ass7XYMyE62/AAQcCBw3YD78VZgo3T/C7l+zMwvN4Tj0fJvkV63vwA2BXYF7gKmDKofTc9/CfjsoI+Hb9mtVpl0RKyIiBvT/SeBW8lm+5RFu6mkw3AwcFdEDG3mZwx4Om2RfkTElRGxLj1cRDaedaDaHI92Dge+FxHPRMQ9ZKMF9h90PySJbKTC+f3Yl3VXqyCdJ2kmsA9wXWr6aPpoe86gywxJAFdKukHS8amt3VTSYTiSDf/xhn08GnqdTjsMHyLL4ht2lXSTpH+XdNAQ9t/qbzGq43EQ8FBE3JFrG/bxqJVaBmlJWwL/DHwyIp4guxLVy4HfI5tD/6UhdOP1EbEv2VWxPiLpDfknI/ssOZTxkZI2AQ4D/ik1jeJ4vMAwj0E7kk4B1gHnpaYVwEsjYh/gU8D/k/TiAXahFH+LnKPY8M182MejdmoXpCVtTBagz4uIiwEi4qGIeD4i1gPfpE8fGzuJiAfSz5XAJWmf7aaSDtofADdGxEOpT0M/Hjl9m05blKRjgEOBD6Q3DFJ54ZF0/wayWvArBtWHDn+LURyPjYB3ARfk+jfU41FHtQrSqZ52NnBrRHw5156vbb4TWNq8bp/7sYWkaY37ZCepltJ+KumgbZAdDft4NOnbdNoiJM0FTgQOi4g1ufZtJU1J919Gdl3guwfYj3Z/i8uAIyVtKmnX1I9BX6DwzcBtEbE817+hHo9aGvWZy2HegNeTfXy+GVicbocA/wgsSe2XAdsPuB8vIzsz/wtgGXBKan8J2QXC7wB+DGwzhGOyBfAIsFWubSjHg+yNYQXwHFlN9bh2x4BsVMffk2VqS4DZA+7HnWQ138br5Btp2Xenv9li4EbgHQPuR9u/BXBKOh63A38wyH6k9vnACU3LDux4+JbdPC3czKzEalXuMDOrGgdpM7MSc5A2MysxB2kzsxJzkDYzKzEHaTOzEnOQNjMrsf8PyCEIBqrKEvMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def porous(ni,nj):\n",
    "    # set up a mesh\n",
    "    mesh = Grid(ni,nj)\n",
    "    mesh.set_extent(200.0,100.0)\n",
    "    mesh.generate()\n",
    "    # set the boundary conditions\n",
    "    mesh.set_Neumann_bc('left')\n",
    "    mesh.set_Neumann_bc('right')\n",
    "    mesh.set_Neumann_bc('bottom')\n",
    "\n",
    "    # now the top boundary condition\n",
    "    mesh.u[-1,:]=100.0 + 0.176 * mesh.x[-1,:]\n",
    "    return mesh    \n",
    "\n",
    "ptest = porous(101,51)\n",
    "print(ptest)\n",
    "\n",
    "# run the solver\n",
    "start = time.process_time()\n",
    "itt, err = SOR(ptest,tol=0.5e-10,maxit=100000)\n",
    "stop = time.process_time()\n",
    "print(\"Solver took {:.3g} seconds.\".format(stop-start))\n",
    "print('Converged after {} itterations, final residual is {}'.format(itt,err))\n",
    "\n",
    "# plot the solution Theleft, right and bottom cells contain nothing useful.\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "cmap = plt.get_cmap('jet')\n",
    "cf = ax1.contourf(ptest.x[1:,1:-1],ptest.y[1:,1:-1],ptest.u[1:,1:-1],levels=21, cmap=cmap)\n",
    "fig.colorbar(cf, ax=ax1)\n",
    "ax1.set_title(f'Porous media head ({ptest.Ni} x {ptest.Nj} grid)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "metric-physics",
   "metadata": {},
   "source": [
    "# Further tests\n",
    "We can compare this with the analytical solution, and conduct a mesh refinement study.  We're going to do the 2nd and need the integrate function again.  I'm interested in the integral of head at z=20.0 m from 25.0 to 175.0 m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "artificial-camping",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate as integrate\n",
    "from refinement_analysis import refinement_analysis\n",
    "\n",
    "def integrate_u_dx(mesh,x0,x1,y):\n",
    "    '''Calculate U=\\int_{1.0}^{1.95)u(x,0.25) dx  using the\n",
    "    u value stored on the grid and simpsons rule'''\n",
    "    \n",
    "    # find the left and right grid points\n",
    "    i0,j = mesh.find((x0,y))\n",
    "    i1,j = mesh.find((x1,y))\n",
    "\n",
    "    # add 1 to i1 as we are going to use it as an upper \n",
    "    # bound forarray slicing\n",
    "    i1 = i1 +1\n",
    "    \n",
    "    # integrate\n",
    "    return integrate.simps(mesh.u[j,i0:i1],mesh.x[j,i0:i1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "collectible-venture",
   "metadata": {},
   "source": [
    "We then need to set up a sequence of grids based on grid doubling. The number of grid points in each direction should be given by $$N =2^n m + 1$$ where n is the grid index and m is an integer.  For our problem we are going to use\n",
    "$$N_i=2^n 20 + 1\\text{ and }N_j=2^n 10 + 1$$\n",
    "\n",
    "The following script creates a sequence of grids and solves them all.  The upper limit on the `range(5,-1,-1)` is the grid index of the finest grid.  Grid 5 takes about 20 minutes to run on my macbook pro. Grid 6, takes a much more reasonable three and a half minutes. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "opponent-brazilian",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Uniform 641x321 grid from (0.0, 0.0) to (200.0, 100.0).\n",
      "Boundaries conditions are - left: Neumann, right: Neumann, top: Dirichlet, bottom: Neumann.\n"
     ]
    }
   ],
   "source": [
    "import datetime # just seconds may not be enough\n",
    "\n",
    "# we need some lists u and dx values\n",
    "U_val = []\n",
    "dx_val = []\n",
    "run_time = []\n",
    "n_pts =[]\n",
    "for grid_index in range(5,1,-1):\n",
    "    ni = 20*2**grid_index + 1\n",
    "    nj = 10*2**grid_index + 1\n",
    "    n_pts.append(ni*nj)\n",
    "    \n",
    "    # set up the problem\n",
    "    test = porous(ni,nj)\n",
    "    print(test)\n",
    "    \n",
    "    # solve it with 12 d.p. precision and a lot of itterations\n",
    "    start = time.process_time()\n",
    "    itt, err = SOR(test,tol=0.5e-10,maxit=1000000)\n",
    "    stop = time.process_time()\n",
    "    print(\"The solver converged after {} itterations, it took {}, final residual is {}\" \\\n",
    "          .format(itt,datetime.timedelta(seconds=int(stop-start)),err))\n",
    "\n",
    "    # save dx and the integral\n",
    "    dx_val.append(test.Delta_x())\n",
    "    U_val.append(integrate_u_dx(test,75.0,125.0,25.0)/50)\n",
    "    run_time.append(stop-start)\n",
    "          \n",
    "    print('Integrated value is ',U_val[-1],'\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "norman-stylus",
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets to the refinement analysis\n",
    "analysis = refinement_analysis(dx_val,U_val)\n",
    "analysis.report(r'\\int_{75}^{125}h(x,25) dx')\n",
    "analysis.plot(True,r'$\\int_{75}^{125}h(x,25) dx$')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "internal-geometry",
   "metadata": {},
   "source": [
    "### Analysis\n",
    "\n",
    "The refinement analysis shows\n",
    "* _p_ = 1.0, so our 2nd order method has an apparent order of accuracy of 1.  Which is very disapointing.\n",
    "* The GCI Ratios show that the solution on none of the girds is safely in the asymptotic region where the error is proprtional to Δx<sup>2</sup>.  This means need either\n",
    "    * Our choice of monitor quanity is poor,\n",
    "    * The SOR solution is not sufficiently converged (i.e. `tol` is too high), or\n",
    "    * The grids are too coarse.\n",
    "\n",
    "The We can try running the simulation on a finer grid, but this will need significantly longer to run.  Doubling the mesh size will increase the runtime by at least eight (2<sup>4</sup>) times. So we can expect it to run for almost three hours! \n",
    "\n",
    "We may need a faster solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "nuclear-inventory",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}