2013-12-03-Alternating_Direction_Implicit_Method.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---\n",
      "layout: default\n",
      "title: \"The Alternating Direction Implicit Method\"\n",
      "date: 2013-12-03\n",
      "author: Georg R Walther\n",
      "tags: wip\n",
      "---"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**This post is work in progress**"
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "The Alternating Direction Implicit Method"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Just as the [Crank-Nicolson (CN) method](http://georg.io/2013/12/03/Crank_Nicolson.html) for reaction-diffusion\n",
      "systems with one space dimension, the \n",
      "[Alternating Direction Implicit (ADI) method](http://en.wikipedia.org/wiki/Alternating_direction_implicit_method)\n",
      "is used commonly for reaction-diffusion systems with two space dimensions.\n",
      "\n",
      "The ADI method has been described thoroughly many times, for instance by\n",
      "[Dehghan](http://www.sciencedirect.com/science/article/pii/S0377042700004520).\n",
      "\n",
      "In the following discussion we will consider a reaction-diffusion system similar to the one we\n",
      "studied [previously](http://georg.io/2013/12/03/Crank_Nicolson.html) and we will use analogous\n",
      "[Neumann boundary conditions](http://en.wikipedia.org/wiki/Neumann_boundary_condition).\n",
      "\n",
      "$$\\frac{\\partial u}{\\partial t} = D \\left( \\frac{\\partial^2 u}{\\partial x^2} + \\frac{\\partial^2 u}{\\partial y^2} \\right) + f(u),$$\n",
      "\n",
      "$$\\frac{\\partial u}{\\partial x}\\Bigg|_{x = 0, L_x} = 0,$$\n",
      "\n",
      "$$\\frac{\\partial u}{\\partial y}\\Bigg|_{y = 0, L_y} = 0,$$\n",
      "\n",
      "where $u(x,y,t)$ is our concentration variable, $D$ is the diffusion coefficient of $u$, $f$ is the reaction term, and $L_x$ and $L_y$ are the extent of our domain in the $x$ and $y$ direction respectively.\n",
      "\n",
      "Since ADI and CN are somewhat related, let us try to derive the most basic properties of the ADI in relation to CN."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Grid Construction"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Analogous to the $(x,t)$-grid\n",
      "[used in CN](http://georg.io/2013/12/03/Crank_Nicolson.html#finite_difference_methods),\n",
      "we need to construct an $(x,y,t)$-grid for this problem with two space dimensions.\n",
      "\n",
      "In the simplest case, we construct a [regular grid](http://en.wikipedia.org/wiki/Regular_grid) as follows\n",
      "\n",
      "$$t_n = n \\Delta t,~ n = 0, \\ldots, N-1,$$\n",
      "\n",
      "$$x_j = j \\Delta x,~ j = 0, \\ldots, J-1,$$\n",
      "\n",
      "$$y_i = i \\Delta y,~ i = 0, \\ldots, I-1,$$\n",
      "\n",
      "where $N$, $J$, and $I$ are the number of grid points in the $t$-, $x$-, and $y$-direction.\n",
      "\n",
      "$\\Delta t$, $\\Delta x$, and $\\Delta y$ are defined as follows\n",
      "\n",
      "$$\\Delta t = \\frac{T}{N-1},~ \\Delta x = \\frac{L_x}{J-1},~ \\Delta y = \\frac{L_y}{I-1},$$\n",
      "\n",
      "where $T$ is the total amount of time we are interested in.\n",
      "\n",
      "[As before]((http://georg.io/2013/12/03/Crank_Nicolson.html#finite_difference_methods) we will refer to the\n",
      "numerical approximations of the unknown analytic solution $u(x,y,t)$ in our grid points as\n",
      "\n",
      "$$U(j \\Delta x, i \\Delta y, n \\Delta t) \\approx u(j \\Delta x, i \\Delta y, n \\Delta t),$$\n",
      "\n",
      "and we use the shorthand $U(j \\Delta x, i \\Delta y, n \\Delta t) = U_{j,i}^n.$\n",
      "We also refer to the grid point $(j \\Delta x, i \\Delta y, n \\Delta t)$ as $(j,i,n)$."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Motivation for ADI"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "When integrating the above reaction-diffusion equation numerically, we can still make use of the\n",
      "[CN stencil](http://georg.io/2013/12/03/Crank_Nicolson.html#the_cranknicolson_stencil).\n",
      "\n",
      "Applying the CN stencil to our reaction-diffusion equation on our $(x,y,t)$-grid, we obtain:\n",
      "\n",
      "$$\\frac{U_{j,i}^{n+1} - U_{j,i}^n}{\\Delta t} = \n",
      "\\frac{D}{2 \\Delta x^2} \\left( U_{j+1,i}^n -\n",
      "2 U_{j,i}^n + U_{j-1,i}^n + U_{j+1,i}^{n+1} - 2 U_{j,i}^{n+1} + U_{j-1,i}^{n+1}\\right) +\n",
      "\\frac{D}{2 \\Delta y^2} \\left( U_{j,i+1}^n -\n",
      "2 U_{j,i}^n + U_{j,i-1}^n + U_{j,i+1}^{n+1} - 2 U_{j,i}^{n+1} + U_{j,i-1}^{n+1}\\right) +\n",
      "f(U_{j,i}^n).$$\n",
      "\n",
      "Let us define $\\sigma_x = \\frac{D \\Delta t}{2 \\Delta x^2}$ and $\\sigma_y = \\frac{D \\Delta t}{2 \\Delta y^2}$.\n",
      "\n",
      "To [reorder our stencil](http://georg.io/2013/12/03/Crank_Nicolson.html#reordering_stencil_into_linear_system) into a linear system\n",
      "we need to define a new index that combines indices $i$ and $j$:\n",
      "\n",
      "$$k = j + i J$$\n",
      "\n",
      "This new index $k$ flattens the two-dimensional spatial part of our $(x,y,t)$-grid into a (one-dimensional) vector\n",
      "and our \"flattening methodology\" is analogous to the [row-major order](http://en.wikipedia.org/wiki/Row-major_order)\n",
      "representation of matrices.\n",
      "\n",
      "To illustrate, when we use the $k$ index grid points $(j,i,n)$ and $(j,i+1,n)$ become $(k,n)$ and $(k+J,n)$ respectively.\n",
      "\n",
      "Our stencil therefore becomes\n",
      "\n",
      "$$U_{k}^{n+1} - U_{k}^n = \n",
      "\\sigma_x \\left( U_{k+1}^n -\n",
      "2 U_{k}^n + U_{k-1}^n + U_{k+1}^{n+1} - 2 U_{k}^{n+1} + U_{k-1}^{n+1}\\right) +\n",
      "\\sigma_y \\left( U_{k+J}^n -\n",
      "2 U_{k}^n + U_{k-J}^n + U_{k+J}^{n+1} - 2 U_{k}^{n+1} + U_{k-J}^{n+1}\\right) +\n",
      "\\Delta t f(U_{k}^n),$$\n",
      "\n",
      "and reordering this expression, we obtain\n",
      "\n",
      "$$\n",
      "-\\sigma_y U_{k-J}^{n+1}\n",
      "-\\sigma_x U_{k-1}^{n+1}\n",
      "+ (1 + 2 \\sigma_x + 2 \\sigma_y) U_{k}^{n+1}\n",
      "-\\sigma_x U_{k+1}^{n+1}\n",
      "-\\sigma_y U_{k+J}^{n+1} =\n",
      "\\sigma_y U_{k-J}^n\n",
      "+\\sigma_x U_{k-1}^n\n",
      "+(1 - 2 \\sigma_x - 2 \\sigma_y) U_{k}^n\n",
      "+\\sigma_x U_{k+1}^n\n",
      "+\\sigma_y U_{k+J}^n\n",
      "+\\Delta t f(U_{k}^n).$$\n",
      "\n",
      "If we [wrote this out in matrix notation](http://georg.io/2013/12/03/Crank_Nicolson.html#reordering_stencil_into_linear_system), we would\n",
      "see \n",
      "[banded matrices](http://publib.boulder.ibm.com/infocenter/clresctr/vxrx/index.jsp?topic=%2Fcom.ibm.cluster.essl.v5r2.essl100.doc%2Fam5gr_bandma.htm) \n",
      "both on the left- and right-hand side (analogous to matrices $A$ and $B$\n",
      "[in the one-dimensional case](http://georg.io/2013/12/03/Crank_Nicolson.html#reordering_stencil_into_linear_system)):\n",
      "\n",
      "On the left-hand side of this liner system, matrix $A$ contains a\n",
      "tridiagonal core [as before](http://georg.io/2013/12/03/Crank_Nicolson.html#reordering_stencil_into_linear_system) with\n",
      "non-zero elements, $-\\sigma_y$, in the $J$-th [sub- and superdiagonals](http://en.wikipedia.org/wiki/Diagonal#Matrices).\n",
      "Matrix $B$ on the right-hand side looks similar but with $\\sigma_y$ in the $J$-th sub- and superdiagonals.\n",
      "\n",
      "Numerical inversion of banded matrices, such as $A$, is [expensive](http://en.wikipedia.org/wiki/Alternating_direction_implicit_method)\n",
      "with standard methods - at the very least we would not be able to use the standard\n",
      "[Thomas algorithm](http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm) for tridiagonal matrices.\n",
      "\n",
      "The [ADI stencil](http://en.wikipedia.org/wiki/Alternating_direction_implicit_method) makes use of a trick that allows us to invert\n",
      "tridiagonal matrices, just as with the \n",
      "[CN stencil](http://georg.io/2013/12/03/Crank_Nicolson.html#the_cranknicolson_stencil),\n",
      "instead of banded matrices.\n",
      "\n",
      "The ADI stencil also suggests a straightforward way to implement a parallelized variant of the ADI method and we will discuss and\n",
      "demonstrate this below."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The ADI Stencils"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The trick used in constructing the ADI method is to split our time step $\\Delta t$ into two and apply two different stencils in each\n",
      "half time step:\n",
      "therefore to increment time by one time step $n \\rightarrow n+1$ in grid point $(j,i,n)$, we first compute $U_{j,i}^{n+1/2}$\n",
      "($n \\rightarrow n+1/2$) and then compute $U_{j,i}^{n+1}$ ($n+1/2 \\rightarrow n+1$).\n",
      "Both of these stencils are chosen such that the resulting linear system is tridiagonal.\n",
      "\n",
      "The two [ADI stencils](http://en.wikipedia.org/wiki/Alternating_direction_implicit_method) \n",
      "bring in the $x$- and $y$-direction at the next time point\n",
      "([implicit](http://en.wikipedia.org/wiki/Finite_difference_method#Implicit_method) v\n",
      "[explicit](http://en.wikipedia.org/wiki/Finite_difference_method#Explicit_method) finite difference stencils) alternatingely -\n",
      "inspiring the name of the ADI method:\n",
      "\n",
      "$$\\frac{U_{j,i}^{n+1/2} - U_{j,i}^n}{\\Delta t / 2} = \n",
      "\\frac{D}{\\Delta x^2} \\left( U_{j+1,i}^{n+1/2} - 2 U_{j,i}^{n+1/2} + U_{j-1,i}^{n+1/2} \\right) +\n",
      "\\frac{D}{\\Delta y^2} \\left( U_{j,i+1}^n - 2 U_{j,i}^n + U_{j,i-1}^n \\right) +\n",
      "f(U_{j,i}^n),$$\n",
      "\n",
      "$$\\frac{U_{j,i}^{n+1} - U_{j,i}^{n+1/2}}{\\Delta t / 2} = \n",
      "\\frac{D}{\\Delta x^2} \\left( U_{j+1,i}^{n+1/2} - 2 U_{j,i}^{n+1/2} + U_{j-1,i}^{n+1/2} \\right) +\n",
      "\\frac{D}{\\Delta y^2} \\left( U_{j,i+1}^{n+1} - 2 U_{j,i}^{n+1} + U_{j,i-1}^{n+1} \\right) +\n",
      "f(U_{j,i}^{n+1/2}).$$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The ADI Linear Systems"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "[As before](http://georg.io/2013/12/03/Crank_Nicolson.html) we define $\\sigma_x = \\frac{D \\Delta t}{2 \\Delta x^2}$ and \n",
      "$\\sigma_y = \\frac{D \\Delta t}{2 \\Delta y^2}$.\n",
      "\n",
      "Reordering both stencils, we obtain the following systems of linear equations\n",
      "\n",
      "$$\n",
      "-\\sigma_x U_{j-1,i}^{n+1/2} + (1+2\\sigma_x) U_{j,i}^{n+1/2} - \\sigma_x U_{j+1,i}^{n+1/2} =\n",
      "\\sigma_y U_{j,i+1}^n + (1-2\\sigma_y) U_{j,i}^n + \\sigma_y U_{j,i-1}^n + \\frac{\\Delta t}{2} f(U_{j,i}^n),\n",
      "$$\n",
      "\n",
      "$$\n",
      "-\\sigma_y U_{j,i+1}^{n+1} + (1+2\\sigma_y) U_{j,i}^{n+1} - \\sigma_y U_{j,i-1}^{n+1} =\n",
      "\\sigma_x U_{j-1,i}^{n+1/2} + (1-2\\sigma_x) U_{j,i}^{n+1/2} + \\sigma_x U_{j+1,i}^{n+1/2} + \\frac{\\Delta t}{2} f(U_{j,i}^{n+1/2}).\n",
      "$$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Family of Linear Systems in the $x$-Direction"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Consider the stencil pictured in \n",
      "[Dehghan Figure 1](http://www.sciencedirect.com/science/article/pii/S0377042700004520#FIG1)\n",
      "and let us think of the matrix defined by $\\left(U_{j,i}^n\\right)$ \n",
      "($j=0,\\ldots,J-1$, $i=0\\ldots,I-1$, $n$ held constant) as a concentration plane.\n",
      "\n",
      "We realize that for fixed index $i$ the first of our two equations above define a\n",
      "system of linear equations similar to the linear system we obtained for\n",
      "[CN with one space dimension](http://georg.io/2013/12/03/Crank_Nicolson.html#reordering_stencil_into_linear_system).\n",
      "As for CN with one space dimension we need to amend certain entries of our linear system to\n",
      "accommodate our Neumann boundary conditions on the edges of our grid.\n",
      "In fact we need to accommodate boundary conditions in both the $x$- and $y$-direction.\n",
      "\n",
      "Our first equation does not need to be modified for grid points that lie \"within\" the\n",
      "spatial part of our grid, i.e. $(j,i,n)$ for $j=1,\\ldots,J-2$, $i=1,\\ldots,I-2$, and $n=0,1,2,\\ldots$:\n",
      "\n",
      "$$\n",
      "-\\sigma_x U_{j-1,i}^{n+1/2} + (1+2\\sigma_x) U_{j,i}^{n+1/2} - \\sigma_x U_{j+1,i}^{n+1/2} =\n",
      "\\sigma_y U_{j,i+1}^n + (1-2\\sigma_y) U_{j,i}^n + \\sigma_y U_{j,i-1}^n + \\frac{\\Delta t}{2} f(U_{j,i}^n).\n",
      "$$\n",
      "\n",
      "Let us take a look at what happens for grid points on the far left and right with $j=0,J-1$ and $i=1,\\ldots,I-2$:\n",
      "\n",
      "$$\n",
      "j=0:~ (1+\\alpha_x) U_{0,i}^{n+1/2} - \\alpha_x U_{1,i}^{n+1/2} =\n",
      "\\alpha_y U_{0,i-1}^n + (1-2\\alpha_y) U_{0,i}^n + \\alpha_y U_{0,i+1}^n + \\Delta t f(U_{0,i}^n),\n",
      "$$\n",
      "\n",
      "$$\n",
      "j=J-1:~ -\\alpha_x U_{J-2,i}^{n+1/2} + (1+\\alpha_x) U_{J-1,i}^{n+1/2} =\n",
      "\\alpha_y U_{J-1,i-1}^n + (1-2\\alpha_y) U_{J-1,i}^n + \\alpha_y U_{J-1,i+1}^n + \\Delta t f(U_{J-1,i}^n).\n",
      "$$\n",
      "\n",
      "On the left-hand side of the linear system defined by these three expressions (imagine varying $j=0,1,\\ldots,J-1$)\n",
      "we can already make out essentially the same matrix as we constructed for\n",
      "[CN with one space dimension](http://georg.io/2013/12/03/Crank_Nicolson.html#reordering_stencil_into_linear_system).\n",
      "Note that this matrix is the same for all indeces $i=0,1,\\ldots,I-1$.\n",
      "\n",
      "This completes incorporating our boundary conditions in the $x$-direction.\n",
      "Before we move on to doing the same with our boundary conditions in the $y$-direction let us\n",
      "state clearly that our $(x,y)$-concentration plane has\n",
      "[right-handed orientation](http://en.wikipedia.org/wiki/Cartesian_coordinate_system#In_two_dimensions)\n",
      "meaning that $x$ increases *left to right* and $y$ increases *bottom to top*.\n",
      "Our grid has the same orientation so that index $j$ incrases *left to right* and\n",
      "index $i$ increases *bottom to top*.\n",
      "\n",
      "Let us now take a look at this equation for $j=1,\\ldots,J-2$ and $i=0$ (bottom) and $i=I-1$ (top):\n",
      "\n",
      "$$\n",
      "i=0:~ -\\alpha_x U_{j-1,0}^{n+1/2} + (1+2\\alpha_x) U_{j,0}^{n+1/2} - \\alpha_x U_{j+1,0}^{n+1/2} =\n",
      "(1-\\alpha_y) U_{j,0}^n + \\alpha_y U_{j,1}^n + \\Delta t f(U_{j,0}^n),\n",
      "$$\n",
      "\n",
      "$$\n",
      "i=I-1:~ -\\alpha_x U_{j-1,I-1}^{n+1/2} + (1+2\\alpha_x) U_{j,I-1}^{n+1/2} - \\alpha_x U_{j+1,I-1}^{n+1/2} =\n",
      "\\alpha_y U_{j,I-2}^n + (1-\\alpha_y) U_{j,I-1}^n + \\Delta t f(U_{j,I-1}^n).\n",
      "$$\n",
      "\n",
      "To rewrite these equations in compact matrix notation, let us define a horizontal slice (*left to right*)\n",
      "through our concentration plane as vector\n",
      "\n",
      "$$\\mathbf{U}_{x,i}^n = \\begin{bmatrix}U_{0,i}^n, & \\ldots, & U_{J-1,i}^n \\end{bmatrix}$$\n",
      "\n",
      "We can now combine all of the above and write the first of our two families of ADI linear systems compactly:\n",
      "\n",
      "$$A \\mathbf{U}_{x,i}^{n+1/2} = \\mathbf{b}_i + \\mathbf{f}\\left( \\Delta t \\mathbf{U}_{x,i}^n \\right),~i=0,\\ldots,I-1,$$\n",
      "\n",
      "where \n",
      "\n",
      "$$A = \\begin{bmatrix}\n",
      "1+\\alpha_x & -\\alpha_x & 0 & 0 & 0 & \\cdots & 0 & 0 & 0 & 0\\\\\n",
      "-\\alpha_x & 1+2\\alpha_x & -\\alpha_x & 0 & 0 & \\cdots & 0 & 0 & 0 & 0 \\\\\n",
      "0 & -\\alpha_x & 1+2\\alpha_x & -\\alpha_x & \\cdots & 0 & 0 & 0 & 0 & 0 \\\\\n",
      "0 & 0 & \\ddots & \\ddots & \\ddots & \\ddots & 0 & 0 & 0 & 0 \\\\\n",
      "0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\alpha_x & 1+2\\alpha_x & -\\alpha_x \\\\\n",
      "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\alpha_x & 1+\\alpha_x\n",
      "\\end{bmatrix}.$$\n",
      "\n",
      "The form of vector $\\mathbf{b}_i$ depends on $i$:\n",
      "\n",
      "$$i=I-1:~ \\mathbf{b}_{I-1} = \n",
      "\\begin{bmatrix}\n",
      "(1-\\alpha_y) U_{0,I-1}^n + \\alpha_y U_{0,I-2}^n \\\\\n",
      "\\vdots \\\\\n",
      "(1-\\alpha_y) U_{J-1,I-1}^n + \\alpha_y U_{J-1,I-2}^n\n",
      "\\end{bmatrix}$$\n",
      "\n",
      "$$i=I-2,\\ldots,1:~ \\mathbf{b}_i = \n",
      "\\begin{bmatrix}\n",
      "\\alpha_y U_{0,i+1}^n + (1-2\\alpha_y) U_{0,i}^n + \\alpha_y U_{0,i-1}^n \\\\\n",
      "\\vdots \\\\\n",
      "\\alpha_y U_{J-1,i+1}^n + (1-2\\alpha_y) U_{J-1,i}^n + \\alpha_y U_{J-1,i-1}^n\n",
      "\\end{bmatrix}$$\n",
      "\n",
      "$$i=0:~ \\mathbf{b}_0 = \n",
      "\\begin{bmatrix}\n",
      "\\alpha_y U_{0,1}^n + (1-\\alpha_y) U_{0,0}^n \\\\\n",
      "\\vdots \\\\\n",
      "\\alpha_y U_{J-1,1}^n + (1-\\alpha_y) U_{J-1,0}^n\n",
      "\\end{bmatrix}.$$\n",
      "\n",
      "The reaction term vector is\n",
      "\n",
      "$$\\mathbf{f}\\left( \\Delta t \\mathbf{U}_{x,i}^n \\right) = \n",
      "\\begin{bmatrix}\n",
      "\\Delta t f(U_{0,i}^n), & \\ldots, & \\Delta t f(U_{J-1,i}^n)\n",
      "\\end{bmatrix}$$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Family of Linear Systems in the $y$-Direction"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us first define a vertical (*from top to bottom*) slice through our concentration plane\n",
      "(note our comment on \n",
      "[right-handedness](http://en.wikipedia.org/wiki/Cartesian_coordinate_system#In_two_dimensions) above) as:\n",
      "\n",
      "$$\\mathbf{U}_{y,j}^n = \\begin{bmatrix}U_{j,I-1}^n, & U_{j,I-2}^n, & \\ldots, & U_{j,0}^n \\end{bmatrix}$$\n",
      "\n",
      "Following an equivalent procedure to above for the second family of ADI linear systems we obtain:\n",
      "\n",
      "$$C \\mathbf{U}_{y,j}^{n+1} = \\mathbf{d}_j + \\mathbf{f}\\left( \\Delta t \\mathbf{U}_{y,j}^{n+1/2} \\right),~j=0,\\ldots,J-1,$$\n",
      "\n",
      "where \n",
      "\n",
      "$$C = \\begin{bmatrix}\n",
      "1+\\alpha_y & -\\alpha_y & 0 & 0 & 0 & \\cdots & 0 & 0 & 0 & 0\\\\\n",
      "-\\alpha_y & 1+2\\alpha_y & -\\alpha_y & 0 & 0 & \\cdots & 0 & 0 & 0 & 0 \\\\\n",
      "0 & -\\alpha_y & 1+2\\alpha_y & -\\alpha_y & \\cdots & 0 & 0 & 0 & 0 & 0 \\\\\n",
      "0 & 0 & \\ddots & \\ddots & \\ddots & \\ddots & 0 & 0 & 0 & 0 \\\\\n",
      "0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\alpha_y & 1+2\\alpha_y & -\\alpha_y \\\\\n",
      "0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\\alpha_y & 1+\\alpha_y\n",
      "\\end{bmatrix}.$$\n",
      "\n",
      "The form of vector $\\mathbf{d}_j$ depends on $j$:\n",
      "\n",
      "$$j=0:~ \\mathbf{d}_0 = \n",
      "\\begin{bmatrix}\n",
      "(1-\\alpha_x) U_{0,I-1}^{n+1/2} + \\alpha_x U_{1,I-1}^{n+1/2} \\\\\n",
      "(1-\\alpha_x) U_{0,I-2}^{n+1/2} + \\alpha_x U_{1,I-2}^{n+1/2} \\\\\n",
      "\\vdots \\\\\n",
      "(1-\\alpha_x) U_{0,0}^{n+1/2} + \\alpha_x U_{1,0}^{n+1/2}\n",
      "\\end{bmatrix}$$\n",
      "\n",
      "$$j=1,\\ldots,J-2:~ \\mathbf{d}_j = \n",
      "\\begin{bmatrix}\n",
      "\\alpha_x U_{j-1,I-1}^{n+1/2} + (1-2\\alpha_x) U_{j,I-1}^{n+1/2} + \\alpha_x U_{j+1,I-1}^{n+1/2} \\\\\n",
      "\\alpha_x U_{j-1,I-2}^{n+1/2} + (1-2\\alpha_x) U_{j,I-2}^{n+1/2} + \\alpha_x U_{j+1,I-2}^{n+1/2} \\\\\n",
      "\\vdots \\\\\n",
      "\\alpha_x U_{j-1,0}^{n+1/2} + (1-2\\alpha_x) U_{j,0}^{n+1/2} + \\alpha_x U_{j+1,0}^{n+1/2}\n",
      "\\end{bmatrix}$$\n",
      "\n",
      "$$j=J-1:~ \\mathbf{d}_{J-1} = \n",
      "\\begin{bmatrix}\n",
      "\\alpha_x U_{J-2,I-1}^{n+1/2} + (1-\\alpha_x) U_{J-1,I-1}^{n+1/2} \\\\\n",
      "\\alpha_x U_{J-2,I-2}^{n+1/2} + (1-\\alpha_x) U_{J-1,I-2}^{n+1/2} \\\\\n",
      "\\vdots \\\\\n",
      "\\alpha_x U_{J-2,0}^{n+1/2} + (1-\\alpha_x) U_{J-1,0}^{n+1/2}\n",
      "\\end{bmatrix}$$\n",
      "\n",
      "The reaction term vector is\n",
      "\n",
      "$$\\mathbf{f}\\left( \\Delta t \\mathbf{U}_{y,j}^{n+1/2} \\right) = \n",
      "\\begin{bmatrix}\n",
      "\\Delta t f(U_{j,I-1}^{n+1/2}), & \\Delta t f(U_{j,I-2}^{n+1/2}), & \\ldots, & \\Delta t f(U_{j,0}^{n+1/2})\n",
      "\\end{bmatrix}$$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Parallelism in the ADI"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To summarize, the ADI stencil generates two families of linear systems that we need to solve iteratively:\n",
      "\n",
      "$$A \\mathbf{U}_{x,i}^{n+1/2} = \\mathbf{b}_i + \\mathbf{f}\\left( \\Delta t \\mathbf{U}_{x,i}^n \\right),~i=0,\\ldots,I-1,$$\n",
      "\n",
      "$$C \\mathbf{U}_{y,j}^{n+1} = \\mathbf{d}_j + \\mathbf{f}\\left( \\Delta t \\mathbf{U}_{y,j}^{n+1/2} \\right),~j=0,\\ldots,J-1.$$\n",
      "\n",
      "Upon closer inspection, we realize that we can solve the $I$ linear systems of the first family in parallel.\n",
      "To see this let us take a look at two arbitrary linear systems out of this family:\n",
      "\n",
      "$$A \\mathbf{U}_{x,i_1}^{n+1/2} = \\mathbf{b}_{i_1} + \\mathbf{f}\\left( \\Delta t \\mathbf{U}_{x,i_1}^n \\right),$$\n",
      "\n",
      "$$A \\mathbf{U}_{x,i_2}^{n+1/2} = \\mathbf{b}_{i_2} + \\mathbf{f}\\left( \\Delta t \\mathbf{U}_{x,i_2}^n \\right).$$\n",
      "\n",
      "As we can see, there is no interdependence between the linear systems for indeces $i_1$ and $i_2$ that would prevent\n",
      "us from solving these in parallel.\n",
      "\n",
      "The same reasoning can be applied to our second family of linear systems.\n",
      "We can solve the linear systems of all indeces $j$ of that family in parallel."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "An ADI Example in Python"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy\n",
      "from matplotlib import pyplot"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "J = 5\n",
      "I = 10\n",
      "\n",
      "print numpy.array([[str((j,i)) for j in range(J)] for i in range(I-1,-1,-1)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[['(0, 9)' '(1, 9)' '(2, 9)' '(3, 9)' '(4, 9)']\n",
        " ['(0, 8)' '(1, 8)' '(2, 8)' '(3, 8)' '(4, 8)']\n",
        " ['(0, 7)' '(1, 7)' '(2, 7)' '(3, 7)' '(4, 7)']\n",
        " ['(0, 6)' '(1, 6)' '(2, 6)' '(3, 6)' '(4, 6)']\n",
        " ['(0, 5)' '(1, 5)' '(2, 5)' '(3, 5)' '(4, 5)']\n",
        " ['(0, 4)' '(1, 4)' '(2, 4)' '(3, 4)' '(4, 4)']\n",
        " ['(0, 3)' '(1, 3)' '(2, 3)' '(3, 3)' '(4, 3)']\n",
        " ['(0, 2)' '(1, 2)' '(2, 2)' '(3, 2)' '(4, 2)']\n",
        " ['(0, 1)' '(1, 1)' '(2, 1)' '(3, 1)' '(4, 1)']\n",
        " ['(0, 0)' '(1, 0)' '(2, 0)' '(3, 0)' '(4, 0)']]\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L_x = 1.\n",
      "L_y = 2.\n",
      "T = 1.\n",
      "\n",
      "J = 100\n",
      "I = 200\n",
      "N = 1000\n",
      "\n",
      "dx = float(L_x)/float(J-1)\n",
      "dy = float(L_y)/float(I-1)\n",
      "dt = float(T)/float(N-1)\n",
      "\n",
      "x_grid = numpy.array([j*dx for j in range(J)])\n",
      "y_grid = numpy.array([i*dy for i in range(I)])\n",
      "t_grid = numpy.array([n*dt for n in range(N)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "D_U = 0.1\n",
      "D_V = 10."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "tot_protein = 2.26"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "no_high = 10\n",
      "U =  numpy.array([[0.1 if (I-1)-i+1 > no_high else 2.0 for j in range(J)] for i in range(I)])\n",
      "\n",
      "U_protein = sum(sum(U))*dx*dy\n",
      "V_protein_px = float(tot_protein-U_protein)/float(I*J*dx*dy)\n",
      "V = numpy.array([[V_protein_px for i in range(0,J)] for i in range(I)])\n",
      "\n",
      "print 'Initial protein mass', (sum(sum(U))+sum(sum(V)))*(dx*dy)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Initial protein mass 2.26\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = pyplot.subplots()\n",
      "ax.set_xlim(left=0., right=(J-1)*dx)\n",
      "ax.set_ylim(bottom=0., top=(I-1)*dy)\n",
      "heatmap = ax.pcolor(x_grid, y_grid, U, vmin=0., vmax=2.1)\n",
      "colorbar = pyplot.colorbar(heatmap)\n",
      "colorbar.set_label('concentration U')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD9CAYAAACcJ53WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1wVNX9P/D3JY9IAgHNppqNBJ8I2WzCLgQMjyqQAB0T\nJNAQS4DQlEwcpXbU0bEPUxzHkWqtWEuM3wo+QTDCaCPBhbFlVTISMEkNofYJsdmNwI8oQyRrwm5y\nf3/EbFmSvbm7dze79+b9mrkz2XvP2XuOnX728LnnniOIoiiCiIhUY0yoG0BERL5h4CYiUhkGbiIi\nlWHgJiJSGQZuIiKVYeAmIlIZycBts9mwYMECGI1GTJ06Fb/97W+HLLd582YYDAaYzWY0NzcHpaFE\nRNQvUupidHQ0tm/fjoyMDFy6dAlmsxl5eXnIyspyl9m3bx/a2tpw8uRJNDc3o7S0FH/729+C3nAi\notFKcsSdlJSEjIwMAEBcXBwyMzPx1VdfeZQ5cOAASkpKAAAmkwkulwt2uz1IzSUiItk57i+//BLH\njx/HvHnzPM7b7XakpKS4P+v1egZuIqIgkkyVDLh06RJWr16Nbdu2IT4+ftD1q9+aFwRhUJmhzhER\neaNkNY5rBAHf+VB+4sSJ+Oabb/y+30gbNnA7nU4UFhbi3nvvxYoVKwZd1+v1sNlsmD17NoD+Ebhe\nrx/yuw4obGy42QXgx6FuRBCwX+qhxT4BwHKF9b8D8KQP5X954YLCO44syVSJKIr4yU9+gvT0dPz8\n5z8fsszy5cuxa9cuAEBTUxMiIiKQnJwc+JYSEfkgyodDbSRH3PX19XjzzTeRmZkJk8kEAHjqqafQ\n1tYGACgvL0dhYSEOHz4Mg8GAmJgY7Ny5M/itJiIahqw8sEpJ9m3evHno6+sb9ktefPHFgDVITYyh\nbkCQsF/qocU+BcrYUDcgiLT8oxR0maFuQJCwX+qhxT4FihpTIHIxcBORJmk5uHGtEiLSJCUPJ5Uu\n92GxWGA0GpGeno6tW7cGtF+Atn+UiGgUUxLclCz30dPTg4qKChw5cgRJSUnIyclBbm6ue4JHIHDE\nTUSapGTErWS5j4aGBhgMBiQnJyMyMhJFRUWoq6sLaN9GdMS9XHOv4BBRcCh9BUf64eQJAK0yv2dg\nuY+rpzp7W+6jvb190Hmr1Sq73XIwVUJEmiQ1HXDW98eAai/lfF3uY6QwVUJEmhTpwzEUuct9DBgY\ngV993mazeYzAA4GBm4g0SUmOW8lyH9nZ2WhtbUV7ezucTidqamqwbNmygPaNqRIi0iQlwU3Jch+x\nsbGorKxEXl4e+vr6UFJSArPZrLQ7HgRxhJI0/cu68uEkEcmxXFH+WBAENPpQfgZCl6/2B0fcRKRJ\nWg5uWu4bEY1iXKuEiEhluDogEZHKcMRNRKQyWg5uWu4bEY1iUb5EN1fQmhEUDNxEpEmRDNxEROoS\nFRHqFgQPAzcRaZJPI26V0XDXiGg0i4oJdQuCh4GbiLRJw9FNw10jolFNw9FNw10jolFNw9FNw10j\nolGNs0qIiFRGw9FNw10jolGNs0qIiFRGw9FNw10jolFNQXTbuHEj6urqoNPpcOLEiUHXn332Wfd+\nky6XC59//jk6OjqQkJCA1NRUjB8/HhEREYiKisKxY8f8b4gX3LqMiMKQ8q3LxLk+lK/33Lrs448/\nRlxcHNatWzdk4L7S/v378fzzz+ODDz4AAEyZMgWNjY2YNGmSX22Xg7u8E5E2RfpwXGX+/PmYOHGi\nrNvs3r0bxcXFHueCPR5m4CYibVIQuOVyOBw4ePAgCgsL3ecEQcCSJUuQmZmJF1980f8vl8AcNxFp\nk0R0s37Tfyj13nvvYd68eUhISHCfO3r0KHQ6Hc6fP4+lS5ciLS0NixcvVn6zKzBwE5E2SUwHvOP6\n/mPAllP+3WLPnj2D0iQ6nQ4AkJiYiFWrVuH48eMBD9xMlRCRNgU5VXLx4kV89NFHKCgocJ9zOBxw\nOBwAgK6uLlgsFhgMBv/74AVH3ESkTQpeeS8uLsaHH36Ijo4OpKSkYMuWLXA6nQCA8vJyAMC7776L\nvLw8jB37v/3kz507hxUrVkAQBDgcDqxZswb5+fmKujEUTgckojAUgOmAxcOXc5evDv5MkEDiiJuI\ntEnD0U3DXSOiUY2rAxIRqYyGo5uGu0ZEo1psqBsQPAzcRKRNTJUQEamMhqObhrtGRKOahqObhrtG\nRKMaUyVERCqj4eim4a4R0aim4eg27CJTGzduRFJSEoxG45DXrVYrJkyYAJPJBJPJhCeffDLgjSQi\n8lmMD4fKDPubVFpaigceeADr1q3zWmbhwoWora0NaMOIiBQZzSNuOVv4qGlxFiIaJUZgB5xQUdxk\nQRDwySefwGg0QqfT4bnnnkNWVpaX0ruu+NsIIFPp7YlIE1oASG/K6zPOKvFuxowZsNvtiI2NxaFD\nh7BixQqcPn3aS+kfK70dEWlSJjwHcru8FZRPhSNpuRTvgBMXF4fY2P5FAXJzcxEdHY2zZ88qbhgR\nkSIaTpUoDtwdHR3uvxsbG9HV1eXec42IKGQifDhUZtjfmuG28KmursbLL78MAIiOjsbu3bsxZgy3\nsiSiENPw6oDcuoyIwlAAti5724fyqz1nx23cuBF1dXXQ6XQ4cWLwQ1Or1YqCggLcdNNNAIDCwkL8\n8pe/BABYLBY88sgj6O3txfr16/Hoo4/63Q9vVJjdISKSQUEKxN/3V3p6elBRUYEjR44gKSkJOTk5\nyM3Nhclk8r8xQ2BOg4i0ScHDSX/fX2loaIDBYEBycjIiIyNRVFSEuro6Zf0YAkfcRKRNEtHN2gJY\nFUwb9/b+it1uR0pKirucXq+H1Wr1/0ZeMHATkTZJpEruMPUfA7b4OG18qPdXvvjiC//a6QcGbiLS\npiDOKomLi3P/PfD+yrlz55CSkgKbzea+ZrPZPEbggcIcNxFpUxDncV/9/sqlS5eg0+mQnZ2N1tZW\ntLe3w+l0oqamBsuWLVPclatxxE1E2qQguvn6/kp1dTXGjBmD2NhYVFZWIi8vD319fSgpKYHZbA5E\nbzxwHjcRhaEAzOP+1IfyM9W1yilH3ESkTRqObhruGhGNaipcg0QuBm4i0iYNRzcNd42IRjUV7iUp\nFwM3EWlTmEe3u+++2+OzIAiIi4vDnXfeibKysu8ndAyNs0qIKAwFYFbJVz6Uv2HkZ5UM9Sp8Z2cn\nqqurodPpsG3bNq91GbiJKAwFIHD/Px/K68JnOmBfXx+MRiNOnjzptUyY/2OCiMg/okpnlYwZM2bY\nzWgYuIlIk3rDPLp98803g84NpEpuvfVWybph3jUiIv+Ee+A2m80eDyAFQcC4ceOwcOFC7NixQ7Ju\nmHeNiMg/PTHRPpS+HLR2ePPll1/6XZeBm4g0qTdCpUluGRi4iUiTejX8zjsDNxFpkouBm4hIXXpV\nFN5Onz4Nu90OURQhiiIEQcCCBQu8lldPz4iIfKCWVMmDDz6Id955BwaDARFX5OUZuIlo1FFL4K6t\nrcW//vUvxMTIXxWLgZuINKkHvkwHDJ309HT09vb6VIeBm4g0SUmOe+PGjairq4NOp8OJEycGXX/j\njTfwzDPPQBRFxMTEoKqqCjNmzAAApKamYvz48YiIiEBUVBSOHTsmea/o6GgYjUYsWrTIPeoWBAEv\nvPCC1zoM3ESkSUpSJaWlpXjggQewbt26Ia9PnToV9fX1iI+Ph8ViQVlZGZqbmwH0B12r1YpJkybJ\nuld+fj7y8/Pdb1EOPJyUwsBNRJqkJHDPnz9f8s3GWbNmuf+eO3cu2tvbPa77stLghg0b0N3djdbW\nVgiCgIyMjGHz3QzcRKRJUvO4G61daLJ2BeQ+VVVVKCgocH8WBAFLliyBy+XCpk2bcP/990vWP3jw\nIDZs2OBeWOrf//43XnvtNeTm5nqtw8BNRJokleOefscETL9jgvvz/23xYfHuK1itVuzYsQP19fXu\nc0ePHoVOp8P58+exdOlSpKWlYfHixV6/4+GHH8bHH3+MW265BQBw6tQprFixYsjc+gDpRV+JiFSq\nFxGyD3+0tLSgrKwMtbW1mDhxovu8TqcDACQmJmLVqlU4fvy45Pf09fW5gzYA3Hzzzejr65OswxE3\nEWnS5SBOB2xra8PKlSvx5ptvegRdh8MBALjmmmvQ1dUFi8WChx56SPK7MjMzUV5ejuLiYoiiiLfe\neguZmZmSdRi4iUiTlKxVUlxcjA8//BAdHR1ISUnBli1b4HQ6AQDl5eV44okncOHCBVRUVACAe9rf\n2bNncc8990AQBDgcDqxZswb5+fmS99qxYwd+//vf45lnngHQ/2D0wQcflKzDPSeJKAwp33PyffEO\n2eWXCdaw2XNSDo64iUiTwv2V99WrV+Ptt99GRkbGoHnbgiCgpaXFa10GbiLSpHAP3Nu2bQMA1NXV\nDRrtD/cCDmeVEJEmuRAh+wiFG264AQCwfft2pKamehzbt2+XrMvATUSadBkxso9QOnTo0KBz7733\nnmQdpkqISJPCPVVSWVmJ7du349SpUzAaje7zDocD06dPl6zLwE1EmhTuW5fde++9WLZsGR577DFs\n3brVneceO3YskpKSJOsycBORJoX71mUTJkzAhAkTsGfPHoiiiDNnzsDlcqGnpwdtbW248cYbvdZl\njpuINCnYr7wHyttvv42bbroJt956KxYuXIjU1FQsW7ZMsg4DNxFpkloC969+9SscP34ct912G06f\nPg2r1Yrbb79dsg4DNxFpkloC97hx43DdddfB6XRCFEUsWLAAn376qWSd8E4CERH5qSfE0/zkGj9+\nPBwOB+bMmYPi4mLodDpERUVJ1uGIm4g0SS0j7v379yMmJgYvvPAC7rrrLtxyyy04ePCgZJ1hA/fG\njRuRlJTkMc/waps3b4bBYIDZbHbvu0ZEFEpqCNy9vb3Iz89HREQEYmNjsWnTJmzevBnXXnutZL1h\nA3dpaSksFovX6/v27UNbWxtOnjyJV155BaWlpb63nogowML9lXcAiIiIQGRkJL799luf6g2b4x5u\n08wDBw6gpKQEAGAymeByuWC326HX631qCBFRIIX7PO4BMTExSE9PR25uLq655hoA/YtMvfDCC17r\nKO6Z3W5HSkqK+7Ner5cI3Luu+NsIQHqXByIaLVoAeN9j0R+hzl3LVVhYiJUrV7pXBBRFcdjVAQPy\nkyR/ScIfB+J2RKQ5mfAcyO3yVlA2tQTuCxcuDNrx5vnnn5eso3hWiV6vh81mc39mmoSIwkEPomUf\nV1MyKcNiscBoNCI9PR1bt24dtp2vvfbaoHOvvPKKZB3FgXv58uXYtav/17GpqQkRERFITk5W+rVE\nRIr0IlL2cTV/J2X09PSgoqICFosFLS0t2Lt3r9eZdtXV1bj77rtx+vRp3H333e5j0aJFSEhIkOzb\nsKmS4TbNLCwsxOHDh2EwGBATE4OdO3cO95VEREGnJFXi76SML774AgaDwT14LSoqQl1dHUwm06Dv\nmDNnDq6//nqcP38eDz/8sMfqgEOVv9Kwgbu6unq4InjxxReHLUNENJKCmeP2Nimjvb190Hmr1Trk\nd0yePBmTJ0/G0aNHfb6/OubLEBH5SGp+9lfW/+CM9T+Kvj9Qu8Lv3r0bjz/+ODo6OtwTOwRBQGdn\np9c6DNxEpElS87iT7khD0h1p7s+NW7zns4cyMClj9uzZAP43Anc6nR6TNWw2m8cIfCiPPfYYDh48\niGnTpsm+P9cqISJNCuYr794mZWRnZ6O1tRXt7e1wOp2oqakZdm3t1NRUn4I2wBE3EWnU5SGm+cnl\n76SM2NhYVFZWIi8vD319fSgpKYHZbJa8l8lkQnFxMfLz8xEd3d9mQRCwcuVKr3UEMVCJmmH0524O\njMStiEj1livKIQuCgB+Jr8ouXyNsCFjO2lcbNmwAMPjFRakZehxxE5EmqWWtkldffdXnOsxxE5Em\nqWFZVwA4efIk5s2bh7S0/oelf//737FlyxbJOgzcRKRJagncGzduxO9+9zuMHTsWADBt2jTU1NRI\n1lHHvyWIiHwUynW2fdHd3e2eVgj057ojIqTbzsBNRJqklhz3pEmT8J///O9loP379w+7A446ekZE\n5CMl0wFH0ksvvYT169fjH//4B2688UYkJibirbfekqzDwE1EmqSWVMnUqVNRX1+Pjo4OiKKIxMTE\nYevw4SQRaZKSZV1H0qOPPorOzk5cd911SExMxMWLF/H4449L1mHgJiJNUsuskoMHD2L8+PHuzxMm\nTMD7778vWYepEiLSpFAHZLl6enrgdDoRFRUFALh8+TK+++47yToM3ESkSWoJ3GvWrMGdd96J0tJS\niKKIV199FcXFxZJ1uFYJEYUh5WuVTBObZJf/XDCHbK0SAHjnnXfwwQcfQBAELFmyBAUFBZLlGbiJ\nKAwpD9y3iZ/JLv8vISukgdtXfDhJRJqkloeTu3fvRmpqKuLi4hAfH4/4+HiPh5VDYY6biDRJLfO4\n/dkBh4GbiDQp1POz5eIOOERE3wt1CkQuf3bAYeAmIk1SS+C+ePEiYmJicOjQIY/z3LqMiFRG+ayS\nCT1nZJe/GHP9oPtZLBY88sgj6O3txfr16/Hoo496XH/22WfdGwa7XC58/vnn6OjoQEJCAlJTUzF+\n/HhEREQgKioKx44d87svQ2HgJqIwpDxwx3Wdl13+0rhEj/v19PQgLS0NR44cQVJSEnJycvDyyy/D\nZDINWX///v14/vnn8cEHHwAApkyZgsbGRkyaNGnYe3/55ZeoqKhAfX09AGD+/Pn44x//iNTUVK91\nOB2QiDSp1xUh+7haQ0MDDAYDkpOTERkZiaKiItTV1Xm91+7duwe97Sj3h2ft2rUoLi7G119/ja+/\n/hpr1qzB2rVrJeswx01EmjRUQB7Qd+RjiPVHvF632+1ISUlxf9br9bBarUOWdTgcOHjwILZv3+4+\nN/AGpMvlwqZNm3D//fd7vde3336LdevWuT+XlJTgmWee8VoeYOAmIo1yOSUeTs6+o/8Y8NunPS73\np3blee+99zBv3jwkJCS4zx09ehQ6nQ7nz5/H0qVLkZaWhsWLFw9Zf9y4caiursaPfvQjAEBNTQ3i\n4+Ml78lUCRFpUl9vpOzjanq9Hjabzf3ZZrN5jMCvtGfPnkFpEp1OBwBITEzEqlWrcPz4ca/tfOON\nN7Bz504kJCRg4sSJeO211/D6669L9o0PJ4koDCl/OIn/OuVXmBzlcb/u7m6kpaWhvr4eOp0Oc+bM\nQVVVFcxms0e1ixcv4qabboLdbnfv0u5wOAAA11xzDbq6urB8+XI89NBDyM/P97s/V+OIm4i0qTtS\n/nGV2NhYVFZWIi8vD1lZWVi5ciXMZjOqqqpQVVXlLvfuu+8iLy/PHbQB4Ny5c8jJycH06dNhMpmw\ncOFCyaC9du1adHZ2uj9fvHgR69evl+waR9xEFIYCMOI+6UN9gxCy1QFNJhOam5uHPXcljriJSJtc\nPhwh1NPTM2jE3d3dLVmHs0qISJtCHJDl+tnPfoaZM2eiqKgIoiiipqYGDz30kGQdpkqIKAwFIFVy\n1If6t4cuVQIATU1N+Mtf/gJBELBo0SKvb2gOYOAmojAUgMBd70P9uaEN3L5iqoSItEklqRJ/MHAT\nkTZJP99TNQZuItImjriJiFSGgZuISGUYuImIVMaHpUrUhoGbiLSpN9QNCB4GbiLSJqZKiIhUhtMB\niYhUhiNuIiKV0XDgHnZZV4vFAqPRiPT0dGzdunXQdavVigkTJsBkMsFkMuHJJ58MSkOJiHyikmVd\n/SE54u7p6UFFRQWOHDmCpKQk5OTkIDc3d9DKVQsXLkRtbW1QG0pE5BMNTweUHHE3NDTAYDAgOTkZ\nkZGRKCoqQl1d3aByalpVi4hGiV4fjiEoyTYMV1cpyRG33W732NlYr9fDarV6lBEEAZ988gmMRiN0\nOh2ee+45ZGVlefnGXVf8bQSQ6V+riUhjWgCcCOxXKphVoiTbILeuEpKBu38NbWkzZsyA3W5HbGws\nDh06hBUrVuD06dNeSv/YnzYSkeZlwnMgt8tbQfkU5K6vzDYAcGcbrg6+Q2Ub5NZVQjJVotfrYbPZ\n3J9tNpvHCBwA4uLiEBsbCwDIzc1FdHQ0zp49G7AGEhH5xenDcZWhsg12u92jzJXZhkWLFuGzzz6T\nXVcpyRF3dnY2Wltb0d7eDp1Oh5qaGo+t6QGgo6MD1113HQCgsbERXV1d0Ol0AW0kEZHPpF55t1kB\nu9XrZX+zDV988YXPzfSHZOCOjY1FZWUl8vLy0NfXh5KSEpjNZnfwLi8vR3V1NV5++WUAQHR0NHbv\n3o0xY7h5PBGFmFSq5Po7+o8BR7d4XJabbRgwkG04d+4cUlJShq2rFPecJKIwFIA9Jyt8qF/puedk\nd3c30tLSUF9fD51Ohzlz5qCqqgpms9ld5upsQ35+Pmw2Gy5fvjxsXaX45iQRaZOCedz+ZBuqq6sx\nZswYr3UDiSNuIgpDARhxl/hQ/w3u8k5EFHoqfJVdLgZuItImDb/yzsBNRNrEHXCIiFSGqRIiIpVh\n4CYiUhnmuImIVKYn1A0IHgZuItImpkqIiFSGqRIiIpXhdEAiIpVhqoSISGUYuImIVIY5biIileF0\nQCIilWGqhIhIZZgqISJSGU4HJCJSGQ2nSrgdOxFpk8uHYwgWiwVGoxHp6enYunXroOtvvPEGMjMz\nYTQaMXPmTDQ2NrqvpaamIjMzEyaTCbNmzQpsv8ARNxFplYIcd09PDyoqKnDkyBEkJSUhJycHubm5\nMJlM7jJTp05FfX094uPjYbFYUFZWhubmZgD9e15arVZMmjRJaS+GxBE3EWmTghF3Q0MDDAYDkpOT\nERkZiaKiItTV1XmUmTVrFuLj4wEAc+fORXt7u8f1YG4+zBE3EY1C1u+PodntdqSkpLg/6/V6WK3e\ny1dVVaGgoMD9WRAELFmyBC6XC5s2bcL999+vvMlXYOAmolHoju+PAVs8rgqCIPubrFYrduzYgfr6\neve5o0ePQqfT4fz581i6dCnS0tKwePFiRS2+ElMlRERX0ev1sNls7s82m81jBD6gpaUFZWVlqK2t\nxcSJE93ndTodACAxMRGrVq3C8ePHA9o+Bm4i0iinD4en7OxstLa2or29HU6nEzU1NVi2bJlHmba2\nNqxcuRJvvvkmbrnlFvd5h8MBh8MBAOjq6oLFYoHBYAhoz5gqISKN8n8id2xsLCorK5GXl4e+vj6U\nlJTAbDajqqoKAFBeXo4nnngCFy5cQEVFBQAgKioKx44dw9mzZ3HPPfdAEAQ4HA6sWbMG+fn5AenR\nAEEM5qPPK28kCAAOjMStiEj1liualdEfby76UGNCUGeBBBpH3ESkUd+FugFBw8BNRBql3VWmGLiJ\nSKO0u1gJAzcRaRRH3EREKsMRNxGRynDETUSkMpxVQkSkMkyVEBGpDFMlREQqwxE3EZHKcMRNRKQy\nHHETEakMR9xERCrD6YBERCrDETcRkcowx01EpDLaHXEPu+ekxWKB0WhEeno6tm7dOmSZzZs3w2Aw\nwGw2o7m5OeCNDF8toW5AkLBf6qHFPgWKy4djMCWxT05dJSQDd09PDyoqKmCxWNDS0oK9e/cOCsz7\n9u1DW1sbTp48iVdeeQWlpaUBb2T4OhHqBgQJ+6UeWuxToPi/WbCS2CenrlKSgbuhoQEGgwHJycmI\njIxEUVER6urqPMocOHAAJSUlAACTyQSXywW73R7QRhIR+c7/EbeS2CenrlKSgdtutyMlJcX9Wa/X\nDwrKcsoQEY2873w4PCmJfe3t7UGPiZIPJ/t3Sh7e1bsje6+3XNb3qcuuUDcgSNgv9dBinwLhN7JL\nxsXFeXz2N/aNFMnArdfrYbPZ3J9tNpvHL8mVZWbPng2g/1dIr9cP+q5QdZCIRh+l8cbf2JeSkgKn\n0zlsXaUkUyXZ2dlobW1Fe3s7nE4nampqsGzZMo8yy5cvx65d/b/4TU1NiIiIQHJyckAbSUQ0kpTE\nPjl1lZIcccfGxqKyshJ5eXno6+tDSUkJzGYzqqqqAADl5eUoLCzE4cOHYTAYEBMTg507dwa0gURE\nI01J7PNWN6DEAHv//ffFjIwMcdq0aeLTTz89ZJkHHnhATE9PF00mk9jU1BToJgTccH16/fXXRaPR\nKGZkZIgzZswQP/300xC00ndy/rcSRVE8duyYGBERIe7bt28EW+cfOX06fPiwmJ2dLWZlZYkLFiwY\n4Rb6Z7h+nTlzRrzrrrvE9PR08bbbbhNfeumlELTSN6WlpaJOpxMzMjK8llFbrBgpAQ3c3d3dYmpq\nqmi320Wn0ynOnDlz0H/svXv3igUFBaIoimJTU5OYlZUVyCYEnJw+NTQ0iJ2dnaIo9v8fbPr06aFo\nqk/k9EsURdHlcol33nmn+MMf/lDcu3dvCFoqn5w+nTlzRjQYDOK5c+dEURTFr7/+OhRN9Ymcfv3i\nF78QH3vsMVEURfH8+fNiQkKC2N3dHYrmyvbRRx+JTU1NXgO32mLFSBr2zUlfaHHet5w+zZo1C/Hx\n8QCAuXPnor29PRRN9YncuaZ/+MMfsGrVKiQmJoaglb6R06c9e/agqKgIOp0OADBp0qRQNNUncvqV\nkpKCzs5OAEBnZycSExMRExMTiubKNn/+fEycONHrdbXFipEU0MCtxXnfvra3qqoKBQUFI9E0ReT0\nq729HX/+859RUVEBQP4UqVCR06d//vOf+Oqrr5CTk4PMzEz86U9/Gulm+kxOv37605/i5MmTuOGG\nG5CVlYVt27aNdDMDTm2xYiQFdJGpwM/7Dj1f2ma1WrFjxw7U19cHsUWBIadfDz74IJ5++mkIggCx\nP602Ai3zn5w+9fb2orW1FX/961/hcDhw++23IycnBwaDYQRa6B85/Xrqqacwffp0WK1WnDp1CkuW\nLMFnn33m/pegWqkpVoykgI64fZn7OMDbvO9wIadPANDS0oKysjLU1tZK/vMvXMjpV2NjI9asWYMp\nU6Zg3759uO+++1BbWzvSTZVNTp9uvPFG5ObmYuzYsbj22muxcOFCtLSE90JNcvp15MgRrF69GgBw\n8803Y8qDddV8AAABSElEQVSUKfj8889HtJ2BprZYMaICmTD/7rvvxMmTJ4t2u128fPmyOHPmTLGx\nsdGjzN69e8UVK1aIoiiKjY2NYmZmZiCbEHBy+vTf//5XvPnmm8VPPvkkRK30nZx+XWnDhg1hP6tE\nTp+amprERYsWiS6XS+zq6hLT09PF5ubmELVYHjn9uu+++8Tf/OY3oiiK4tmzZ8Uf/OAH7gew4ez0\n6dOSDyfVFCtGUkBTJVqc9y2nT0888QQuXLjgzgVHRUXh2LFjoWz2sOT0S23k9MlkMmHp0qXIzMyE\n0+lEWVkZpk+fHuKWS5PTr1//+tdYu3Yt0tPT0dvbiyeffNL9ADZcFRcX48MPP0RHRwdSUlKwZcsW\nOJ39K/WpMVaMJEEUwzxxSUREHgKa4yYiouBj4CYiUhkGbiIilWHgJiJSGQZuIiKVYeAmIlKZ/w9y\nLjdyR2qrTwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x3252490>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "alpha_x_U = D_U*dt/(2.*dx*dx)\n",
      "alpha_y_U = D_U*dt/(2.*dy*dy)\n",
      "\n",
      "alpha_x_V = D_V*dt/(2.*dx*dx)\n",
      "alpha_y_V = D_V*dt/(2.*dy*dy)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A_U = numpy.diagflat([-alpha_x_U for j in range(J-1)], -1)+\\\n",
      "    numpy.diagflat([1.+alpha_x_U]+[1.+2.*alpha_x_U for j in range(J-2)]+[1.+alpha_x_U], 0)+\\\n",
      "    numpy.diagflat([-alpha_x_U for j in range(J-1)], 1)\n",
      "    \n",
      "C_U = numpy.diagflat([-alpha_y_U for i in range(I-1)], -1)+\\\n",
      "      numpy.diagflat([1.+alpha_y_U]+[1.+2.*alpha_y_U for i in range(I-2)]+[1.+alpha_y_U], 0)+\\\n",
      "      numpy.diagflat([-alpha_y_U for i in range(I-1)], 1)\n",
      "        \n",
      "A_V = numpy.diagflat([-alpha_x_V for j in range(J-1)], -1)+\\\n",
      "    numpy.diagflat([1.+alpha_x_V]+[1.+2.*alpha_x_V for j in range(J-2)]+[1.+alpha_x_V], 0)+\\\n",
      "    numpy.diagflat([-alpha_x_V for j in range(J-1)], 1)\n",
      "    \n",
      "C_V = numpy.diagflat([-alpha_y_V for i in range(I-1)], -1)+\\\n",
      "      numpy.diagflat([1.+alpha_y_V]+[1.+2.*alpha_y_V for i in range(I-2)]+[1.+alpha_y_V], 0)+\\\n",
      "      numpy.diagflat([-alpha_y_V for i in range(I-1)], 1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f_vec = lambda U, V: numpy.multiply(float(dt)/float(2.), numpy.subtract(numpy.multiply(V, \n",
      "                     numpy.add(k0, numpy.divide(numpy.multiply(U,U), numpy.add(1., numpy.multiply(U,U))))), U))\n",
      "\n",
      "k0 = 0.067"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "b_t_stencil_U = numpy.array([[(1.-alpha_y_U) for j in range(J)],\n",
      "                             [alpha_y_U for j in range(J)]])\n",
      "b_c_stencil_U = numpy.array([[alpha_y_U for j in range(J)],\n",
      "                             [1.-2.*alpha_y_U for j in range(J)],\n",
      "                             [alpha_y_U for j in range(J)]])\n",
      "b_b_stencil_U = numpy.array([[alpha_y_U for j in range(J)],\n",
      "                             [(1.-alpha_y_U) for j in range(J)]])\n",
      "\n",
      "b_t_stencil_V = numpy.array([[(1.-alpha_y_V) for j in range(J)],\n",
      "                             [alpha_y_V for j in range(J)]])\n",
      "b_c_stencil_V = numpy.array([[alpha_y_V for j in range(J)],\n",
      "                             [1.-2.*alpha_y_V for j in range(J)],\n",
      "                             [alpha_y_V for j in range(J)]])\n",
      "b_b_stencil_V = numpy.array([[alpha_y_V for j in range(J)],\n",
      "                             [(1.-alpha_y_V) for j in range(J)]])\n",
      "\n",
      "d_l_stencil_U = numpy.array([[1.-alpha_x_U, alpha_x_U] for i in range(I)])\n",
      "d_c_stencil_U = numpy.array([[alpha_x_U, 1.-2.*alpha_x_U, alpha_x_U] for i in range(I)])\n",
      "d_r_stencil_U = numpy.array([[alpha_x_U,1.-alpha_x_U] for i in range(I)])\n",
      "                             \n",
      "d_l_stencil_V = numpy.array([[1.-alpha_x_V, alpha_x_V] for i in range(I)])\n",
      "d_c_stencil_V = numpy.array([[alpha_x_V, 1.-2.*alpha_x_V, alpha_x_V] for i in range(I)])\n",
      "d_r_stencil_V = numpy.array([[alpha_x_V, 1.-alpha_x_V] for i in range(I)])\n",
      "\n",
      "f_curr = f_vec(U,V)\n",
      "\n",
      "def b_U(i):\n",
      "    if i <= I-2 and i >= 1:\n",
      "        U_y = U[[i+1, i, i-1], :]\n",
      "        return numpy.add(U_y+b_c_stencil_U, f_curr[[i+1, i, i-1], :]).sum(axis=0)\n",
      "    elif i == I-1:\n",
      "        U_y = U[[I-1, I-2], :]\n",
      "        return numpy.add(U_y+b_t_stencil_U, f_curr[[I-1, I-2], :]).sum(axis=0)\n",
      "    elif i == 0:\n",
      "        U_y = U[[1, 0], :]\n",
      "        return numpy.add(U_y+b_b_stencil_U, f_curr[[1, 0], :]).sum(axis=0)\n",
      "    \n",
      "def b_V(i):\n",
      "    if i <= I-2 and i >= 1:\n",
      "        V_y = V[[i+1, i, i-1], :]\n",
      "        return numpy.subtract(V_y+b_c_stencil_V, f_curr[[i+1, i, i-1], :]).sum(axis=0)\n",
      "    elif i == I-1:\n",
      "        V_y = V[[I-1, I-2], :]\n",
      "        return numpy.subtract(V_y+b_t_stencil_V, f_curr[[I-1, I-2], :]).sum(axis=0)\n",
      "    elif i == 0:\n",
      "        V_y = V[[1, 0], :]\n",
      "        return numpy.subtract(V_y+b_b_stencil_V, f_curr[[1, 0], :]).sum(axis=0)\n",
      "    \n",
      "def d_U(j):\n",
      "    if j <= J-2 and j >= 1:\n",
      "        U_x = U[:, [j-1, j, j+1]]\n",
      "        return numpy.add(U_x+d_c_stencil_U, f_curr[:, [j-1, j, j+1]]).sum(axis=1)\n",
      "    if j == 0:\n",
      "        U_x = U[:, [0, 1]]\n",
      "        return numpy.add(U_x+d_l_stencil_U, f_curr[:, [0, 1]]).sum(axis=1)\n",
      "    if j == J-1:\n",
      "        U_x = U[:, [J-2, J-1]]\n",
      "        return numpy.add(U_x+d_r_stencil_U, f_curr[:, [J-2, J-1]]).sum(axis=1)\n",
      "    \n",
      "def d_V(j):\n",
      "    if j <= J-2 and j >= 1:\n",
      "        V_x = V[:, [j-1, j, j+1]]\n",
      "        return numpy.subtract(V_x+d_c_stencil_V, f_curr[:, [j-1, j, j+1]]).sum(axis=1)\n",
      "    if j == 0:\n",
      "        V_x = V[:, [0, 1]]\n",
      "        return numpy.subtract(V_x+d_l_stencil_V, f_curr[:, [0, 1]]).sum(axis=1)\n",
      "    if j == J-1:\n",
      "        V_x = V[:, [J-2, J-1]]\n",
      "        return numpy.subtract(V_x+d_r_stencil_V, f_curr[:, [J-2, J-1]]).sum(axis=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for n in range(N):\n",
      "    print 'time = %.6f' % (n*dt)\n",
      "    print 'total protein = %.6f' % ((dx*dy)*(sum(sum(U))+sum(sum(V))))\n",
      "    f_curr = f_vec(U,V)\n",
      "    \n",
      "    for i in range(I):\n",
      "        U[i, :] = numpy.linalg.solve(A_U, b_U(i))\n",
      "        V[i, :] = numpy.linalg.solve(A_V, b_V(i))\n",
      "    for j in range(J):\n",
      "        U[:, j] = numpy.linalg.solve(C_U, d_U(j))\n",
      "        V[:, j] = numpy.linalg.solve(C_V, d_V(j))\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = pyplot.subplots()\n",
      "ax.set_xlim(left=0., right=(J-1)*dx)\n",
      "ax.set_ylim(bottom=0., top=(I-1)*dy)\n",
      "heatmap = ax.pcolor(x_grid, y_grid, U, vmin=0., vmax=2.1)\n",
      "colorbar = pyplot.colorbar(heatmap)\n",
      "colorbar.set_label('concentration U')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}