{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Lecture 23: PDEs2 - Relaxation Methods\n",
    "#### Phys 503, J.R. Gladden, University of Mississippi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Relaxation Methods\n",
    "Relaxation methods are used for equilibrium boundary value problems.  They require an iterative solution that \"relaxes\" until the solution values over the entire domain stop changing more than some tolerance value. An example of a typical PDE to solve using such a method would be the 2D Laplace equation for the electrostatic potential:\n",
    "$$ \\mu \\left( \\frac{\\partial^2 \\Phi}{\\partial x^2} + \\frac{\\partial^2 \\Phi}{\\partial y^2} \\right)=0 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### The Solver and Imports\n",
    "This solver provides three methods: Jacobi, Gauss-Seidel, and Simulatoneous Overrelaxation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib wx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from pylab import *\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "import time, sys\n",
    "\n",
    "\n",
    "def solver(rhs,u0,eps=1e-3,method='J',omega='opt'):\n",
    "\t#check array dimensions and method options are OK\n",
    "\tif u0.shape != (Nx,Ny): \n",
    "\t\tprint \"X and Y dimensions are not consistent with initial guess for solution\"\n",
    "\t\tprint \"U0 shape is: \", u0.shape\n",
    "\t\tprint \"X and Y sizes are: \", Nx, Ny\n",
    "\t\tsys.exit()\n",
    "\tif method not in ['J', 'GS', 'SOR']:\n",
    "\t\tprint '''\n",
    "\t\tRelxation method must be:\n",
    "\t\tJ: Jacobi, \n",
    "\t\tGS: Gauss-Seidel, or \n",
    "\t\tSOR: Sim. Overrelaxation'''\n",
    "\t\tsys.exit()\n",
    "\n",
    "\tu=u0.copy()\n",
    "\tunew=u0.copy()\n",
    "\titer=0\n",
    "\tchangeSum=0.0\n",
    "\tchange=1.0\n",
    "\t\n",
    "\t#Compute optimal omega and error check\n",
    "\tif method == 'SOR':\n",
    "\t\tif omega == 'opt': \n",
    "\t\t\tr=0.5*(cos(np.pi/float(Nx)) + cos(np.pi/float(Ny)))\n",
    "\t\t\tomega = 2.0/(1.0+sqrt(1.0-r**2))\n",
    "\t\t\n",
    "\t\tprint \"Using SOR method with omega = \", omega\n",
    "\t\tprint\n",
    "\t\t\t\n",
    "\t\tif omega < 1.0 or omega > 2.0:\n",
    "\t\t\tprint \"Relaxation parameter is out of bounds.  Should be between 1 and 2. Quitting...\"\n",
    "\t\t\treturn 1\n",
    "\t\t\n",
    "\tstartTime=time.time()\n",
    "\twhile change > eps:\n",
    "\t\titer +=1\n",
    "\t\tchangeSum=0.0\n",
    "\t\tfor i in range(1,Nx-1):\n",
    "\t\t\tfor j in range(1,Ny-1):\n",
    "\t\t\t\tif method == 'J':\n",
    "\t\t\t\t\tunew[i,j]=0.25*(u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1])\n",
    "\t\t\t\tif method == 'GS':\n",
    "\t\t\t\t\tunew[i,j] = 0.25*(u[i+1,j] + unew[i-1,j] + u[i,j+1] + unew[i,j-1])\n",
    "\t\t\t\tif method == 'SOR':\n",
    "\t\t\t\t\tunew[i,j] = (1.0-omega)*u[i,j]+omega/4.0*(u[i+1,j] + unew[i-1,j] + u[i,j+1] + unew[i,j-1])\n",
    "\t\t\t\tif unew[i,j]!=0.0: changeSum+=abs(1.0-u[i,j]/unew[i,j])\n",
    "\t\tu = unew.copy()\n",
    "\t\t# changeSum has added up deltaU for all interior points, so now need to divide by number\n",
    "\t\t# of interior points to get change per point\n",
    "\t\tchange=changeSum/float(Nx-2)**2\n",
    "\t\tif iter % 100 == 0: \n",
    "\t\t\tprint \"For iteration %i, the sum of change in solution is: %2.3e. \" % (iter,change)\n",
    "\t\n",
    "\tendTime=time.time()\n",
    "\tprint '='*40\n",
    "\tprint \"Solution obtained after %i iterations in %2.3f seconds\" % (iter, endTime-startTime)\n",
    "\tprint '='*40\n",
    "\treturn unew\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Initial Conditions\n",
    "The closer the initial guess is to the final solution, the faster it will converge.  This allows for several different choices for initial conditions based on the boundary values: (1) all zeros, (2) mean value of boundaries, (3) linear (plane) extrapolation between boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def initialGuess(u0,method='mid'):\n",
    "\n",
    "\tif method=='mid':\n",
    "\t\tmidval=(u0.max() - u0.min())/2.0\n",
    "\t\tu0[1:Nx-2,1:Ny-2]=midval\n",
    "\n",
    "\tif method == 'plane':\n",
    "\t\tYslope=(average(u0[Ny-1,:])-average(u0[0,:]))/Ly\n",
    "\t\tXslope=(average(u0[:,Nx-1])-average(u0[:,0]))/Lx\n",
    "\t\tfor i in range(1,Nx-2):\n",
    "\t\t\tfor j in range(1,Ny-2):\n",
    "\t\t\t\tu0[i,j]=Xslope*x[j]+average(u0[:,0])\n",
    "\tif method=='average':\n",
    "\t\tXavg1 = average(u0[:,0])\n",
    "\t\tXavg2 = average(u0[:,Nx-1])\n",
    "\t\tYavg1 = average(u0[0,:])\n",
    "\t\tYavg2 = average(u0[Ny-1,:])\n",
    "\t\tavg1 = (Xavg1+Yavg1)/2.0\n",
    "\t\tavg2 = (Xavg2+Yavg2)/2.0\n",
    "\t\tmidval = (avg1+avg2)/2.0\n",
    "\t\tu0[1:Nx-2,1:Ny-2]=midval\n",
    "\treturn u0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Usage and Visualization\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For iteration 100, the sum of change in solution is: 2.305e-03. \n",
      "For iteration 200, the sum of change in solution is: 1.164e-03. \n",
      "For iteration 300, the sum of change in solution is: 6.595e-04. \n",
      "For iteration 400, the sum of change in solution is: 4.017e-04. \n",
      "For iteration 500, the sum of change in solution is: 2.535e-04. \n",
      "For iteration 600, the sum of change in solution is: 1.632e-04. \n",
      "For iteration 700, the sum of change in solution is: 1.066e-04. \n",
      "========================================\n",
      "Solution obtained after 716 iterations in 3.534 seconds\n",
      "========================================\n"
     ]
    }
   ],
   "source": [
    "Nx=Ny=50\n",
    "Lx=Ly=10.0\n",
    "x=linspace(0,Lx,Nx)\n",
    "y=linspace(0,Ly,Ny)\n",
    "\n",
    "#Apply Boundary Conditions and initial guess\n",
    "u0=zeros((Nx,Ny))\n",
    "u0[0,:]=0.5\n",
    "u0[:,0]=1.0 #+0.25*sin(5*pi*x/Lx)\n",
    "\n",
    "u0=initialGuess(u0,method='mid')\n",
    "\n",
    "ufinal=solver(0,u0,eps=1e-4,method='GS')\n",
    "\n",
    "\n",
    "## Surface Plotting\n",
    "\n",
    "close('all')\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.gca(projection='3d')\n",
    "\n",
    "X,Y=meshgrid(x,y)\n",
    "\n",
    "#surf0=ax.plot_surface(X,Y,u0,rstride=2,cstride=2,linewidth=0.1,alpha=0.2)\n",
    "surf=ax.plot_surface(X,Y,ufinal,rstride=1,cstride=1,cmap=cm.jet,linewidth=0.1)\n",
    "fig.colorbar(surf,shrink=0.5,aspect=5)\n",
    "xlabel('X')\n",
    "ylabel('Y')\n",
    "#zlabel('Z')\n",
    "show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}