ion-transport-diagnosis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equations that we are trying solve are those of the paper titled \"Diffuse-charge dynamics in electrochemical systems\" by Bazant, Thornton and Ajdari published in https://journals.aps.org/pre/abstract/10.1103/PhysRevE.70.021506"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from fipy import *\n",
    "m = Grid1D(nx=100,Lx=2.0)\n",
    "v = 1\n",
    "epsilon= 0.05\n",
    "delta = 0.1\n",
    "x = m.cellCenters[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The set of equations is "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\frac{\\partial c}{\\partial t}\n",
    "&= \\varepsilon\\frac{\\partial}{\\partial x} \\left(\\frac{\\partial c}{\\partial x} + \\rho\\frac{\\partial \\phi}{\\partial x}\\right)\\\\\n",
    "\\frac{\\partial \\rho}{\\partial t} \n",
    "&=  \\varepsilon\\frac{\\partial}{\\partial x} \\left(\\frac{\\partial\\rho}{\\partial x}  + c \\frac{\\partial \\phi}{\\partial x}\\right)\\\\\n",
    "-\\varepsilon^2 \\frac{\\partial^2\\phi}{\\partial x^2} &= \\rho\n",
    "\\end{split}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the initial conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "c(x,0) &= 1\\\\\n",
    "\\rho(x,0) &= 0\\\\\n",
    "\\phi(x,0) &= v(x-1)\\\\\n",
    "\\end{split}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the boundary conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\frac{\\partial c}{\\partial x} + \\rho\\frac{\\partial \\phi}{\\partial x} &= 0\\\\\n",
    "\\frac{\\partial \\rho}{\\partial x} + c\\frac{\\partial \\phi}{\\partial x} &= 0\\\\\n",
    "\\text{where both fluxes are zero at} &\\quad\\text{$x = 0$ and $x=2$}\\\\\n",
    "\\end{split}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\phi = +v - \\varepsilon\\delta \\frac{\\partial \\phi}{\\partial x} \\quad\\text{at $x = 2$}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\phi = -v + \\varepsilon\\delta \\frac{\\partial \\phi}{\\partial x} \\quad\\text{at $x = 0$}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "##Initial conditions\n",
    "c = CellVariable(name=\"c\", mesh=m , value=1.0)\n",
    "\n",
    "rho = CellVariable(name=\"rho\", mesh=m , value=0.0)\n",
    "\n",
    "phi = CellVariable(name=\"phi\",mesh=m)\n",
    "phi.setValue(v*(x-1), where=m.cellCenters)\n",
    "\n",
    "#vi =Viewer((phi, c, rho),datamin=-1.1, datamax=1.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##Boundary conditions\n",
    "c.faceGrad.constrain(-rho.faceValue*phi.faceGrad, where=m.facesLeft)\n",
    "c.faceGrad.constrain(-rho.faceValue*phi.faceGrad, where=m.facesRight)\n",
    "\n",
    "rho.faceGrad.constrain(-c.faceValue*phi.faceGrad, where=m.facesLeft)\n",
    "rho.faceGrad.constrain(-c.faceValue*phi.faceGrad, where=m.facesRight)\n",
    "\n",
    "phi.faceGrad.constrain((phi.faceValue+v)/(epsilon*delta), where=m.facesLeft)\n",
    "phi.faceGrad.constrain((phi.faceValue-v)/(-epsilon*delta), where=m.facesRight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##Equation\n",
    "\n",
    "Drho = epsilon*rho\n",
    "Dc   = epsilon*c\n",
    "\n",
    "eq1 = TransientTerm(var=c) == DiffusionTerm(coeff=epsilon, var=c) +DiffusionTerm(coeff=Drho, var=phi)\n",
    "eq2 = TransientTerm(var=rho) == DiffusionTerm(coeff=epsilon, var=rho) +DiffusionTerm(coeff=Dc, var=phi)\n",
    "eq3 = DiffusionTerm(coeff=epsilon*epsilon, var=phi)+ImplicitSourceTerm(coeff=1.0, var=rho) == 0\n",
    "eqns = eq1 & eq2 & eq3 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##Solve\n",
    "#vi =Viewer((phi, c, rho),datamin=-1.1, datamax=1.1)\n",
    "#from builtins import range\n",
    "\n",
    "res = 1e+10\n",
    "restol= 1e-13\n",
    "\n",
    "for t in range(100):\n",
    "    c.updateOld()\n",
    "    rho.updateOld()\n",
    "    phi.updateOld()\n",
    "    while res > restol:\n",
    "        res = eqns.sweep(dt=1e-15)\n",
    "        print(res)\n",
    "    \n",
    "#if __name__=='__main__':\n",
    "#        vi.plot()\n",
    "#        raw_input(\"No-flux -- stationary. Press <return> to proceed..\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However the solution method above gave some errors, originated with the warning \"self.factorizeFnc(*args, **kwargs)\" "
   ]
  },
  {
   "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": 2
}