Lecture-20-Introduction-to-FiPy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 20: Introduction to FiPy - Getting to Know the Diffusion Equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Objectives:\n",
    "\n",
    " * Understand how to create the diffusion equation in FiPy.\n",
    " * Be able to change variables in the equation and observe the effects in the diffusion equation solution.\n",
    " * Understand how to save the results to a data file.\n",
    " \n",
    "First thing we'll do is use the qt backend for interacting with matplotlib.  There is documentation in matplotlib about backends.  One or more of the available backends may be installed with your python distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Setting up Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib osx\n",
    "from fipy import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib\n",
    "from fipy import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    " * L is the LENGTH of our box.  You can set this to any value you choose however, appropriate scaling of the problem would admit 1 as the length of choice.\n",
    " * nx is the number of grid points we wish to have in our solution.  If we have fewer we sacrifice accuracy, if we have more, the computational time increases.  You should always check that your solution does not depend on the number of grid points and the grid spacing!\n",
    " * dx is the spacing between grid points.  Similar comments as above in \\#2\n",
    " * timeStepDuration is the amount of 'time' at each step of the solution.  Accuracy and stability of the solution depend on choices for the timeStepDuration and the grid point spacing and the diffusion coefficient.  We will not discuss stability any further, just know that it is something that needs to be considered.\n",
    " * steps is the number of timesteps in dt you wish to run.  You can change this value to observe the solution in different stages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Initializing the Simulation Domain and Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "L = 1.\n",
    "nx = 200\n",
    "dx = L / nx\n",
    "timeStepDuration = 0.001\n",
    "steps = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Here we set the diffusion coefficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "D11 = 0.5   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note:  The 'mesh' command creates the mesh (gridpoints) on which we will solve the equation.  This is specific to FiPy.\n",
    "\n",
    "At this point, if you are in the IPython notebook I would suggest you try the following in the cell below:\n",
    "\n",
    " * Put your cursor to the right of the \"(\" and hit TAB.  You will get the docstring for the Grid1D function.\n",
    " * Do this again with your cursor to the right of the \"d\" in Grid1D().\n",
    " * And again after the \"G\".\n",
    "\n",
    "This is a powerful way to explore the available functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "mesh = Grid1D(dx = dx, nx = nx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We assign to 'c' objects that are 'CellVariables'.  This is a special type of variable used by FiPy to hold the values for the concentration in our solution.  We also create a viewer here so that we can inspect the values of 'c'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "c = CellVariable(name = \"c\", mesh = mesh)\n",
    "viewer = MatplotlibViewer(vars=(c,),datamin=-0.1, datamax=1.1, legend=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note:  This command sets 'x' to contain a list of numbers that define the x position of the grid-points.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CellVariable(value=array([[ 0.0025,  0.0075,  0.0125,  0.0175,  0.0225,  0.0275,  0.0325,\n",
       "         0.0375,  0.0425,  0.0475,  0.0525,  0.0575,  0.0625,  0.0675,\n",
       "         0.0725,  0.0775,  0.0825,  0.0875,  0.0925,  0.0975,  0.1025,\n",
       "         0.1075,  0.1125,  0.1175,  0.1225,  0.1275,  0.1325,  0.1375,\n",
       "         0.1425,  0.1475,  0.1525,  0.1575,  0.1625,  0.1675,  0.1725,\n",
       "         0.1775,  0.1825,  0.1875,  0.1925,  0.1975,  0.2025,  0.2075,\n",
       "         0.2125,  0.2175,  0.2225,  0.2275,  0.2325,  0.2375,  0.2425,\n",
       "         0.2475,  0.2525,  0.2575,  0.2625,  0.2675,  0.2725,  0.2775,\n",
       "         0.2825,  0.2875,  0.2925,  0.2975,  0.3025,  0.3075,  0.3125,\n",
       "         0.3175,  0.3225,  0.3275,  0.3325,  0.3375,  0.3425,  0.3475,\n",
       "         0.3525,  0.3575,  0.3625,  0.3675,  0.3725,  0.3775,  0.3825,\n",
       "         0.3875,  0.3925,  0.3975,  0.4025,  0.4075,  0.4125,  0.4175,\n",
       "         0.4225,  0.4275,  0.4325,  0.4375,  0.4425,  0.4475,  0.4525,\n",
       "         0.4575,  0.4625,  0.4675,  0.4725,  0.4775,  0.4825,  0.4875,\n",
       "         0.4925,  0.4975,  0.5025,  0.5075,  0.5125,  0.5175,  0.5225,\n",
       "         0.5275,  0.5325,  0.5375,  0.5425,  0.5475,  0.5525,  0.5575,\n",
       "         0.5625,  0.5675,  0.5725,  0.5775,  0.5825,  0.5875,  0.5925,\n",
       "         0.5975,  0.6025,  0.6075,  0.6125,  0.6175,  0.6225,  0.6275,\n",
       "         0.6325,  0.6375,  0.6425,  0.6475,  0.6525,  0.6575,  0.6625,\n",
       "         0.6675,  0.6725,  0.6775,  0.6825,  0.6875,  0.6925,  0.6975,\n",
       "         0.7025,  0.7075,  0.7125,  0.7175,  0.7225,  0.7275,  0.7325,\n",
       "         0.7375,  0.7425,  0.7475,  0.7525,  0.7575,  0.7625,  0.7675,\n",
       "         0.7725,  0.7775,  0.7825,  0.7875,  0.7925,  0.7975,  0.8025,\n",
       "         0.8075,  0.8125,  0.8175,  0.8225,  0.8275,  0.8325,  0.8375,\n",
       "         0.8425,  0.8475,  0.8525,  0.8575,  0.8625,  0.8675,  0.8725,\n",
       "         0.8775,  0.8825,  0.8875,  0.8925,  0.8975,  0.9025,  0.9075,\n",
       "         0.9125,  0.9175,  0.9225,  0.9275,  0.9325,  0.9375,  0.9425,\n",
       "         0.9475,  0.9525,  0.9575,  0.9625,  0.9675,  0.9725,  0.9775,\n",
       "         0.9825,  0.9875,  0.9925,  0.9975]]), mesh=UniformGrid1D(dx=0.005, nx=200))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = mesh.cellCenters\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Setting Initial Conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "c.setValue(0.0) \n",
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "c.setValue(0.2, where=x < L/2.) \n",
    "c.setValue(0.8, where=x > L/2.)\n",
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note: Boundary conditions come in two types.  Fixed flux and fixed value.  The syntax is: \n",
    "\n",
    " * FixedValue(mesh.getFacesLeft(), VALUE)\n",
    " * FixedFlux(mesh.getFacesLeft(), FLUX)\n",
    "\n",
    "Fixed value boundaries, can set VALUE or FLUX as a float "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "boundaryConditions=(FixedValue(mesh.facesLeft,0.2),  \n",
    "                   FixedValue(mesh.facesRight,0.8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This line defines the diffusion equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "eqn1 = TransientTerm() == ImplicitDiffusionTerm(D11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "for step in range(1000): \n",
    "    eqn1.solve(c, boundaryConditions = boundaryConditions, dt = timeStepDuration)\n",
    "    viewer.plot()  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "These lines print out a text file with the final data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from fipy.viewers.tsvViewer import TSVViewer  \n",
    "TSVViewer(vars=(c)).plot(filename=\"output.txt\")  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "c\r\n",
      "x\tc\r\n",
      "0.0025\t0.201499999758458\r\n",
      "0.0075\t0.204499999278958\r\n",
      "0.0125\t0.207499998803152\r\n",
      "0.0175\t0.210499998333946\r\n",
      "0.0225\t0.213499997860399\r\n",
      "0.0275\t0.216499997387362\r\n",
      "0.0325\t0.219499996922307\r\n",
      "0.0375\t0.222499996462214\r\n"
     ]
    }
   ],
   "source": [
    "!head output.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Concentration \n",
    "# Dependent Diffusion\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Standard Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib osx\n",
    "from fipy import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### The parameters of our system.\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "* L is the LENGTH of our box.  You can set this to any value you choose\n",
    "however, 1 is the easiest.\n",
    "* `nx` is the number of grid points we wish to have in our solution.  If we have fewer we sacrifice accuracy, if we have more, the computational time increases and we are subject to roundoff error.  You should always check that your solution does not depend on the number of grid points and the grid spacing!\n",
    "* `dx` is the spacing between grid points.  Similar comments as above.\n",
    "* timeStepDuration is the amount of 'time' at each step of the solution.\n",
    "\n",
    "Accuracy and stability of the solution depend on choices for the `timeStepDuration` and the grid point spacing and the diffusion coefficient.  We will not discuss stability any further, just know that it is something that needs to be considered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "L = 1.\n",
    "nx = 50\n",
    "dx = L / nx\n",
    "timeStepDuration = 0.001\n",
    "steps = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note: In the first equation, the diffusion coefficient is constant, concentration independent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "D1 = 3.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note: You have seen all of the following code before.  This time we are solving two simultaneous equations, eqn1 and eqn2.\n",
    "\n",
    "Note:  The 'mesh' command creates the mesh (gridpoints) on which we will solve the equation.  This is specific to FiPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "mesh = Grid1D(dx = dx, nx = nx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We assign to 'c' objects that are 'CellVariables'.  This is a special type of variable used by FiPy to hold the values for the concentration in our solution.  c1 for eqn1 and c2 for eqn2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "c1 = CellVariable(\n",
    "    name = \"c1\",\n",
    "    mesh = mesh,\n",
    "    hasOld = True)\n",
    "\n",
    "c2 = CellVariable(\n",
    "    name = \"c2\",\n",
    "    mesh = mesh,\n",
    "    hasOld = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note:  This command sets 'x' to contain a list of numbers that define the x position of the grid-points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x = mesh.cellCenters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "These lines set diffusant in the initial condition.  They are set to the same values for easy comparison.  Feel free to change these values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "c1.setValue(0.8)\n",
    "c1.setValue(0.2, where=x > L/3.)\n",
    "c2.setValue(0.8)\n",
    "c2.setValue(0.2, where=x > L/3.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "viewer = Matplotlib1DViewer(vars = (c1,c2), limits = {'datamin': 0., 'datamax': 1.})\n",
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Boundary conditions can be either fixed flux or fixed value.  Here, fixed value is used for simple comparison between the diffusion profiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "boundaryConditions=(FixedValue(mesh.facesLeft,0.8), FixedValue(mesh.facesRight,0.2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "boundaryConditions=(FixedFlux(mesh.facesLeft,0.0), FixedFlux(mesh.facesRight,0.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note: In the second equation, the diffusion coefficient is non-constant and is a function of concentration in the system.  So we use D22_0 and D22_1 as the end points of our function. The function, given by \"D\" is simply a linear interpolation of the two D values.  You are, of course, free to try other functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "D22_0 = 3.0\n",
    "D22_1 = 0.5\n",
    "D2 = (D22_1 - D22_0)*c2 + D22_0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note: These are the two diffusion equations.  The first equation is as in previous code, using a concentration independent D1.  The second equation uses a non-constant D described above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "eqn1 = TransientTerm() == ImplicitDiffusionTerm(D1) \n",
    "eqn2 = TransientTerm() == ImplicitDiffusionTerm(D2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The following is an if loop that waits for user input before executing.  Iterates for number of steps stated earlier for each equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "for step in range(10):\n",
    "    c1.updateOld()\n",
    "    c2.updateOld()\n",
    "    eqn1.solve(c1, boundaryConditions = boundaryConditions, dt = timeStepDuration)\n",
    "    eqn2.solve(c2, boundaryConditions = boundaryConditions, dt = timeStepDuration)\n",
    "    viewer.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}