FiPy bug

README.md

Bug with cylindrical gradients in FiPy

Untitled.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  2.66666667e+10   8.88888889e+09   5.33333333e+09 ...,   1.45799162e+07\n",
      "    1.45639905e+07   1.45480997e+07]\n",
      " [  0.00000000e+00   0.00000000e+00   0.00000000e+00 ...,   0.00000000e+00\n",
      "    0.00000000e+00   0.00000000e+00]]\n",
      "[[ 0.  0.  0. ...,  0.  0.  0.]\n",
      " [ 0.  0.  0. ...,  0.  0.  0.]]\n",
      "[ 100000.  100000.  100000. ...,  100000.  100000.  100000.]\n"
     ]
    }
   ],
   "source": [
    "import fipy\n",
    "\n",
    "L = 89e-3  # m\n",
    "\n",
    "diameterChamber = 13.75e-3\n",
    "\n",
    "radius = diameterChamber/2\n",
    "\n",
    "dz = 1e-3\n",
    "\n",
    "dr = 0.75e-5\n",
    "\n",
    "nz = round(L/dz)\n",
    "\n",
    "nr = round(radius/dr)\n",
    "\n",
    "#mesh = Grid2D(dx=dx, nx=nx, dy=dy, ny=ny)\n",
    "\n",
    "mesh = fipy.CylindricalGrid2D(dz=dz, nz=nz, dr=dr, nr=nr)\n",
    "\n",
    " \n",
    "\n",
    "pressure_0 = fipy.CellVariable(mesh=mesh, name='Pressure Gas', hasOld=1)\n",
    "\n",
    "pressure_0.setValue(1e5)\n",
    "\n",
    "print pressure_0.grad.value\n",
    "print pressure_0.leastSquaresGrad\n",
    "print pressure_0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating the gradient for constant $\\phi$\n",
    "\n",
    "$$ \\int \\nabla \\phi dV = \\int \\nabla \\vec{n} \\phi dA $$\n",
    "\n",
    "so this integral is\n",
    "\n",
    "$$ \\phi \\theta \\left(r + \\frac{\\Delta r}{2} \\right) \\Delta z - \\phi \\theta (r - \\frac{\\Delta r}{2}) \\Delta z - \\sin \\frac{\\theta}{2} \\phi \\Delta r \\Delta z - \\sin \\frac{\\theta}{2} \\phi \\Delta r \\Delta z$$ \n",
    "\n",
    "$\\theta$ is the small angle approximation, which drops out when we divide by the volume (which we need to do to calculate $\\nabla \\phi$ eventually). The FiPy formulation does not have the $\\sin \\frac{\\theta}{2}$ terms from the cylindicrial directed faces! This works out to be zero if we approximate for $\\sin$ terms. Currently in FiPy we have an error of $\\phi / r$ as we don't have the $\\sin$ terms, which is bad (after dividing through by the volume, $\\Delta z \\Delta r r \\theta$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:fipy]",
   "language": "python",
   "name": "conda-env-fipy-py"
  },
  "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": 2
}