notebooks/advection-diffusion-reaction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advection-diffusion-reaction equations #\n",
    "\n",
    "The problems we have encountered so far---with the notable exception\n",
    "of the Navier-Stokes equations---all share a common feature: they all\n",
    "involve models expressed by a *single* scalar or vector PDE. In many\n",
    "situations the model is instead expressed as a system of PDEs,\n",
    "describing different quantities possibly governed by (very) different\n",
    "physics. As we saw for the Navier-Stokes equations, one way to solve\n",
    "a system of PDEs in FEniCS is to use a splitting method where we solve\n",
    "one equation at a time and feed the solution from one equation into\n",
    "the next. However, one of the strengths with FEniCS is the ease by\n",
    "which one can instead define variational problems that couple several\n",
    "PDEs into one compound system. In this section, we will look at how to use\n",
    "FEniCS to write solvers for such systems of coupled PDEs.\n",
    "The goal is to demonstrate how easy it is to implement fully implicit,\n",
    "also known as monolithic, solvers in FEniCS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PDE problem ##\n",
    "\n",
    "Our model problem is the following system of\n",
    "advection-diffusion-reaction equations:\n",
    "\\begin{align}\n",
    "  \\label{ftut1-reactionsystem-system-1}\n",
    "  \\frac{\\partial u_1}{\\partial t} +\n",
    "  \\boldsymbol{w} \\cdot \\boldsymbol{\\nabla} u_1 - \\boldsymbol{\\nabla}\\cdot(\\epsilon\\boldsymbol{\\nabla} u_1)\n",
    "    &= f_1 -K u_1 u_2, \\\\\n",
    "  \\label{ftut1-reactionsystem-system-2}\n",
    "  \\frac{\\partial u_2}{\\partial t} +\n",
    "  \\boldsymbol{w} \\cdot \\boldsymbol{\\nabla} u_2 - \\boldsymbol{\\nabla}\\cdot(\\epsilon\\boldsymbol{\\nabla} u_2)\n",
    "    &= f_2 -K u_1 u_2, \\\\\n",
    "  \\label{ftut1-reactionsystem-system-3}\n",
    "  \\frac{\\partial u_3}{\\partial t} +\n",
    "  \\boldsymbol{w} \\cdot \\boldsymbol{\\nabla} u_3 - \\boldsymbol{\\nabla}\\cdot(\\epsilon\\boldsymbol{\\nabla} u_3)\n",
    "    &= f_3 + K u_1 u_2 - K u_3.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This system models the chemical reaction between two species $A$ and\n",
    "$B$ in some domain $\\Omega$:\n",
    "\n",
    "\\[\n",
    "  A + B \\rightarrow C.\n",
    "\\]\n",
    "\n",
    "The variables $u_1$, $u_2$, and $u_3$ denote the\n",
    "concentrations of the three species:\n",
    "\\[\n",
    "  u_1 = [A], \\quad u_2 = [B], \\quad u_3 = [C].\n",
    "\\]\n",
    "The $f_i$ are the source terms of the three species."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The chemical reactions are accounted for in the\n",
    "right-hand sides of the PDE system (\\ref{ftut1-reactionsystem-system-1})-(\\ref{ftut1-reactionsystem-system-3}). We assume that the reaction is *first-order*, meaning that the\n",
    "reaction rate is proportional to the concentrations $[A]$ and $[B]$ of\n",
    "the two species $A$ and $B$:\n",
    "\n",
    "\\[\n",
    "  \\frac{\\mathrm{d}}{\\mathrm{d}t} [C] = K [A] [B].\n",
    "\\]\n",
    "\n",
    "We also assume that the formed species $C$ spontaneously decays with a\n",
    "rate proportional to the concentration $[C]$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The chemical reactions take part at each point in the domain\n",
    "$\\Omega$. In addition, we assume that the species $A$, $B$, and $C$\n",
    "diffuse throughout the domain with diffusivity $\\epsilon$ (the terms\n",
    "$-\\boldsymbol{\\nabla}\\cdot(\\epsilon\\boldsymbol{\\nabla} u_i)$) and are advected with velocity\n",
    "$\\boldsymbol{w}$ (the terms $\\boldsymbol{w}\\cdot\\boldsymbol{\\nabla} u_i$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make things visually and physically interesting, we shall let the\n",
    "chemical reaction take place in the velocity field computed from the\n",
    "solution of the incompressible Navier-Stokes equations around a\n",
    "cylinder. In summary, we will thus be\n",
    "solving the following coupled system of nonlinear PDEs:\n",
    "\n",
    "\\begin{align}\n",
    "  \\label{ftut1-reactionsystem-full}\n",
    "  \\varrho\\left(\\frac{\\partial \\boldsymbol{w}}{\\partial t} +\n",
    "      \\boldsymbol{w} \\cdot \\boldsymbol{\\nabla} \\boldsymbol{w}\\right) &= \\boldsymbol{\\nabla}\\cdot\\boldsymbol{\\sigma}(\\boldsymbol{w}, p) + \\boldsymbol{f}, \\\\\n",
    "  \\boldsymbol{\\nabla} \\cdot \\boldsymbol{w} &= 0, \\\\\n",
    "  \\frac{\\partial u_1}{\\partial t} +\n",
    "  \\boldsymbol{w} \\cdot \\boldsymbol{\\nabla} u_1 - \\boldsymbol{\\nabla}\\cdot(\\epsilon\\boldsymbol{\\nabla} u_1)\n",
    "    &= f_1 - K u_1 u_2, \\\\\n",
    "  \\frac{\\partial u_2}{\\partial t} +\n",
    "  \\boldsymbol{w} \\cdot \\boldsymbol{\\nabla} u_2 - \\boldsymbol{\\nabla}\\cdot(\\epsilon\\boldsymbol{\\nabla} u_2)\n",
    "    &= f_2 - K u_1 u_2, \\\\\n",
    "  \\frac{\\partial u_3}{\\partial t} +\n",
    "  \\boldsymbol{w} \\cdot \\boldsymbol{\\nabla} u_3 - \\boldsymbol{\\nabla}\\cdot(\\epsilon\\boldsymbol{\\nabla} u_3)\n",
    "    &= f_3 + K u_1 u_2 - K u_3.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This system is one-way coupled from the Navier--Stokes subsystem\n",
    "to the advection-diffusion-reaction subsystem.\n",
    "\n",
    "We assume that $u_1 = u_2 = u_3 = 0$ at $t = 0$ and inject the species\n",
    "$A$ and $B$ into the system by specifying nonzero source terms $f_1$\n",
    "and $f_2$ close to the corners at the inflow, and take $f_3 = 0$. The\n",
    "result will be that $A$ and $B$ are convected by advection and\n",
    "diffusion throughout the channel, and when they mix the species $C$\n",
    "will be formed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variational formulation ##\n",
    "\n",
    "We obtain the variational formulation of our system by multiplying\n",
    "each equation by a test function, integrating the second-order terms\n",
    "$-\\boldsymbol{\\nabla}\\cdot(\\epsilon\\boldsymbol{\\nabla} u_i)$ by parts, and summing up the\n",
    "equations. When working with FEniCS it is convenient to think of the\n",
    "PDE system as a vector of equations. The test functions are collected in\n",
    "a vector too, and the variational formulation is the inner product of\n",
    "the vector PDE and the vector test function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need introduce some discretization in time. We will use the\n",
    "backward Euler method as before when we solved the heat equation and\n",
    "approximate the time derivatives by $(u_i^{n+1}-u_i^n) / \\delta t$. Let\n",
    "$v_1$, $v_2$, and $v_3$ be the test functions, or the components of\n",
    "the test vector function. The inner product results in\n",
    "\\begin{align}\n",
    "  \\label{ftut1-reactionsystem-varproblem}\n",
    "  & \\int_{\\Omega}\n",
    "  (\\delta t^{-1} (u_1^{n+1} - u_1^n) v_1 + \\boldsymbol{w} \\cdot \\boldsymbol{\\nabla} u^{n+1}_1 \\, v_1\n",
    "  + \\epsilon \\boldsymbol{\\nabla} u^{n+1}_1 \\cdot \\boldsymbol{\\nabla} v_1) \\rm{d}\\boldsymbol{x} \\\\\n",
    "  + & \\int_{\\Omega} (\\delta t^{-1} (u_2^{n+1} - u_2^n) v_2\n",
    "  + \\boldsymbol{w} \\cdot \\boldsymbol{\\nabla} u^{n+1}_2 \\, v_2\n",
    "  + \\epsilon \\boldsymbol{\\nabla} u^{n+1}_2 \\cdot \\boldsymbol{\\nabla} v_2) \\rm{d}\\boldsymbol{x} \\nonumber \\\\\n",
    "  + & \\int_{\\Omega} (\\delta t^{-1} (u_3^{n+1} - u_3^n) v_3\n",
    "  + \\boldsymbol{w} \\cdot \\boldsymbol{\\nabla} u^{n+1}_3 \\, v_3\n",
    "  + \\epsilon \\boldsymbol{\\nabla} u^{n+1}_3 \\cdot \\boldsymbol{\\nabla} v_3) \\rm{d}\\boldsymbol{x} \\nonumber \\\\\n",
    "  - & \\int_{\\Omega} (f_1 v_1 + f_2 v_2 + f_3 v_3) \\rm{d}\\boldsymbol{x} \\nonumber \\\\\n",
    "  - & \\int_{\\Omega} (-K u^{n+1}_1 u^{n+1}_2 v_1 - K u^{n+1}_1\n",
    "  u^{n+1}_2 v_2 + K u^{n+1}_1 u^{n+1}_2 v_3 - K u^{n+1}_3 v_3) \\rm{d}\\boldsymbol{x} = 0.\n",
    "  \\nonumber\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this problem it is natural to assume homogeneous Neumann boundary\n",
    "conditions on the entire boundary for $u_1$, $u_2$, and $u_3$; that\n",
    "is, $\\partial u_i/\\partial n = 0$ for $i = 1, 2, 3$. This means that\n",
    "the boundary terms vanish when we integrate by parts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FEniCS implementation ##\n",
    "\n",
    "The first step is to read the mesh from file. If you have written out the mesh previously, you can read it in using the following statement:\n",
    "```\n",
    "mesh = Mesh('navier_stokes_cylinder/cylinder.xml.gz')\n",
    "```\n",
    "The mesh is stored in the native FEniCS XML format (with additional\n",
    "gzipping to decrease the file size).\n",
    "\n",
    "Here, we will generate the mesh from scratch. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f2f8b4843c8>,\n",
       " <matplotlib.lines.Line2D at 0x7f2fc23b42e8>]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAABhCAYAAAAtOOHMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH/lJREFUeJztnXlsHNd9x7+/5S4p3rdWEnWQOinqti7Lsh07li8ZsWPY\nCBSgh4sUQVoETeG2QNwCcRPAqGsBqes2QGqkAdqgaOLKjSE3ChxfsSFDckzbkq3DkiVZJ21JlGRa\nS0okl3z9g3ybt7NzvDfzZnZ2+T6AoOXM7MybN29/83u/6xFjDAaDwWAoLxLFboDBYDAY9GOEu8Fg\nMJQhRrgbDAZDGWKEu8FgMJQhRrgbDAZDGWKEu8FgMJQhRrgbDAZDGWKEu8FgMJQhRrgbDAZDGZIs\n1oXb2tpYZ2dnsS5vMBgMJcm7777bzxhr9zpOSrgT0T0A/hlABYCfMMaedDjuIQA7AKxnjPW6nbOz\nsxO9va6HGAwGg8ECEZ2SOc5TuBNRBYAfAbgTwFkA7xDRTsbYIctx9QC+A+Bt9ebKs3v3brz66qtY\nsWIFPvzww7x9ra2tuHTpEhYtWoSPP/44b9+CBQtw/Phx1NfX4+rVq3n7Vq1ahf379+dtSyaTyGaz\n6OjowLlz5/L2LVmyBEeOHMldT2Tp0qU4fPhw3raamhoMDQ1h+vTpuHDhgu3xzc3NuHLlCgCgrq4O\nmUwG3d3d+Oijj/KOnzNnDs6cOYPZs2fj7NmzefuWL1+OAwcOgIjAawY1NTXh888/x+LFi3H06NG8\n42fNmoW+vj60tbWhv78/b19PTw8OHcp7xKisrMTIyEju/kXmzp2L06dP564HABUVFRgbG7N9Vu3t\n7bh48aLts1q4cCGOHTuGxsZGDAwM2LYrlUphdHQUAHLP1K5P+H23tLTg8uXLefuWLVuGgwcP5m2b\nOXMmPv3009x4EeFt5fdldx2RdDqN8+fP5/4X4WNOfFb8fu3GEL83uzHEn3sikcD4+Hjevvnz5+PE\niRN52/gYEp8Vh/evXbvs+oQ/K7sxxNslUlVVheHh4Vw/252rtrYWg4ODtu0S4ePX7VmJ9+F2Lt4X\ndn1iN975b8Fu/PI+t5MPN9xwA9577z1s3boV69evR5iQV+EwItoE4O8ZY3dP/v0YADDG/sFy3NMA\nXgbwNwD+2ktzX7duHfOjuW/fvh1DQ0PK3zMYDIa4kE6n8a1vfcvXd4noXcbYOq/jZByqHQDOCH+f\nndwmXuwGAHMYY79SaqUPbrrpJtvtqVQq7EsbyoRksmiuphxVVVXFboIhZGprax33rVvnKZsDEzha\nhogSAH4I4K8kjv0mEfUSUe/Fixd9Xc9pptHc3OzrfJw5c+bk/d3d3V1wTFNTE5qamvIemnhce7u7\nj4OICrY1NDTkfb+trS1v/6ZNm2yPVb2O6jniSnV1de6z6gv9xhtvBABks1mtbeLce++9nsds3rwZ\nLS0tGB4exqxZsxyP6+rqQn19fcF4CEJFRQUqKiq0nc/K5s2bQzt3GNTV1YV6fusLvKenJ/fZaj4L\ng8BmGSJqBHAcQGbyKzMAXAZwv5tpxq9Z5vvf/77t9oaGBnzxxRfK5zMYDIZi8Pjjj/v6nk6zzDsA\nFhFRFxFVAtgGYCffyRgbYIy1McY6GWOdAPbCQ7CHgaxgv/nmm32d30lLFLX4adOmoaWlJff3woUL\nAUw4w6KksrIy1POraP+y2pHbbEMHbtPgRMLfBLa+vl7p+Nra2rzxUQpMmzZN+li7mYg404qaMMaU\n22wrbniOasZYFsC3AbwE4DCA5xhjB4noB0R0f9gNtLJly5ZA39+9e7ev7/HIDCvcq79q1Spcv34d\nly9fxoIFC9Dc3Ixjx44BQEEER9iMjIxIHdfU1IR0Oq18fpUZUiaT8T4I7i9dHfkQvb29jj9Mv1Nk\na9SViN0LdnBwsCBiJ+5cv37ddX8qlcKCBQsAAH19fQDyzRHXrl3zvMZdd92l3VxIRGhublZ6OcnA\n79EvXBm87777dDTHFU+zTFj4Ncs89dRTUgOmmMybNw9DQ0MQ/Qq33nor3nzzzSK2yrBmzRq8//77\ngc6xYcMG/O53v9PUoujhYbmlgl3YqRc8jDlu8BBnYGLW9+ijj/o6j06zTKzQJdjFKVtFRYXym1TU\neK2O1FOnTsHqMDaCvRA7p3WYBBXsAEpasAMITbBbHbW33nprYOdtKpVSFuxAeA7zoIizWLdZny5K\nTrjfcsstAIIJhrlz5+ZF3YyNjWHXrl1K5xCTUvxG/oio2m9LHSIqSNDyYvr06SG1prx56KGHAn1f\nJmzTKoTfeustV8G8ZMkSz3M6mUI5GzdulLLp6ww7JSLX3+rcuXPz/nbyt915553a2uREyQl3/qCs\nGWMqnD59umDbxo0bfZ8vCDU1NQCieZOrEqYw5fetgjUzM240NTVpP6c1RNcPzz//fKDvDw8PK3/H\nKtitjvUgv1/O22+/jWvXrqGxsdHxmEQi4dj+22+/XfmajDFXP5JVtlj9ba2trQDCD3oASlC4v/LK\nKwAmOnnZsmWBz8c7ee/evb6+HzQO+bbbbpM6rhhRB2EKU2t6OWfevHmhXVMWv1EW1rR1WZwcvYlE\nAh0dHbb7okLHuCMitLa25qKSdLywRKwlKkTcnOWvv/56wbaampqcg9iJIH5KXo7gV78KPd+z9IS7\niLUuiB3iFIp/FqdpspElTljraaiya9cuKQ25WE5ktyy7IMc6TVdPnZqoibR69Wrpc6nCNWz+Yhdn\nETU1NYF+vG44RYQ4RWCMj49j7969WLx4MdauXWt7TNimqqDjrqenB4wxnDp1KtevZ86c8fhW8Rga\nGiqoUyPD/Pnz8+SK22wiKkpOuKuGQormDv7ZzzQzTLw0BUBPOKAqjY2Njhq2HfxYmfAzPl11yuq1\nFpzSCdewR0ZG0NjYmDd7CjOSxG+S3dGjR/Huu+/a7ou7qUos0BXkpZlMJkMr2eC3HIWY+3LixIk8\nueI2mwCiCYUsOeFurS5YDuzZs8dxHx9AJ0+e1HpNGfuw1wB1wis2WsTJGR1mxENXV1fu88DAQIEz\nXTVRJWz7qXEkT4wHq1Kmy1QpjjWVhDYvh68bUcixkhPu1tKpTuhOXtBNS0sL6urqCmzMS5cuzfs7\nyAByw699WCQOBbj88Mknn7juV01UCWrac6O9vT2nnYc91ffrEE4mk3l1U3Qg4/cIw1QZRc0XwD6o\nQzclJ9xXrFghdRzXHq2hSbKE/XK4fPkyMplMzsbMq1361ZajpqqqKrB2HUVlxLCLQ4WNOLMJe2x4\nvfCdtNpsNltQHz0oxUquDBPROc6L2IVJyQl3njwkTq3d8POGTKVSuOeee5S/54UozKxaLzfNBE1v\njoLKykotfosofB980ROvWYZdzZdi1kXhiC+n1tZWz8qjIrpnVuPj45g5c2agc8j2qahc+a39Exa8\nAq1qX5w7dy73/HhIZJjEq9ck4KGQXlPrIIyOjuKFF17I26YjyUgUZlbNxBrWGWezUphmiDD46KOP\nPGcZdjVfZBLlwhY8mUwG9fX1aGxsxKVLl3Dx4kXU1tZKlbiWnVnNmDEj97m5udk1HNW6epIqsqYU\n0W8ThanE6zny/YlEIrdimp++4DMxEwoZI3QnGVmTPKzRIXxwW6NkUqmU1hrfcSHMant+BbCYQetU\nFTQKwXP16tU8k8zg4GBOwLghY6NPJBK4dOkSUqkU7rzzTly5ciVnKgwLJ0VJNb+AV13Vcaz1OVpn\nGHx/sfMOVJiSwt3JDvulL30p4pZ4Y42SGR0dVZpF+F2hyk+1SD/MnDkTzc3NoZqj/AjgqqqqPLso\nd2wvXrxYW7tU8JPRK74QnPwb4+PjGB0dxejoKF5++WXf7ZOltrbWUVFys7Pb3T+vuspxG+vWY71w\nmmE45bXEMbhgSgp3p/ThN954A+3t7XnT1DiiYpLyG20jG5UUlE8//VRKC9WBbCGrBx98EMPDw7kM\nRvF71kWwo0Im/t5N85Wx8fK6TWFilzch81xkSgIHiSzr6elBIpHwtKM7ZcXHsVjZlBTuYlxyXV1d\n3o+CMeaoWUQZeVGsNWGLkSwVFbIVBn/5y1/m/R3m0nQ6cdN8ZWZGbvkWYSLzXD777DNf506lUlK/\npUOHDmF8fNzTju6nTHmxKDnhLlNNzgvuEEwmk8hkMnk/iv7+fsesTNmFJ3QQVny7F+fOnSvKdcPC\nmmAk468Q7f9EhKamplAcp2GvPqWKjPaZTCZ9rSZVrIiXuXPnKpXFcENn9FSY5TU4JSfcdRYdEgez\n3bR17dq1+PKXv4w77rgDy5cv13ZdN/zYVnWi86UShxhz8X4aGxulagFxLbe7uxuMMVy4cCEUx2mp\nxHKLs7lsNutrNSmn/rOOd91j5vjx47bx+zU1NcqZxTqSprgfJwrHbMmtxOS0QLZOvvrVr2LlypUF\nmhVjDLt378Zrr72m7Vpr1651rBsSZ2pra5XqzhSTadOmSZdE6Ozs1F7qoVSor6+3NUmKKwgBE+UQ\ndNa0IaLQXnTJZBLpdDrSGWlDQ4NUHaE4LJA9pdi4cSNWrVplO2UmItxyyy1ay9JaBXupLNrBBbuO\ntS+9zBOrVq0KdH4nwW4XKhhFWrgudJkbOE6+Jm6G4REhXLDX1dV5hlvalTRYs2ZN3t9+zFOypqFs\nNusq2FV/y/Pnz3fcx802fgvE6UZKuBPRPUR0hIiOEdF3bfY/SkSHiOgDInqViEIryh10gWwvZDJT\nH3nkkUDXcBuYMvH0MgksYSL+ON0GspudVaxFwrU2u35pamrC/v37/TTTE7vY+jiUapXFa+Yklt7w\nYycHJnwU/IVntclnMhnPkgjcJLJp06bctvfffz9PoKuYvHgkm66FxlVj+k+cOOG4T8Vss3XrVqXr\n+sFTuBNRBYAfAbgXQA+ArxORtUrQ+wDWMcZWAtgB4CndDeXoKHjlxAMPPCB9bBDHrtfAdLK780SM\nqEIHnZBdi9QtSsGuFoldv7j5IFRCVu20w8OHDxdsK3bfOuGV5GWtj9/a2poTyps2bXIcc14zxaDr\nFXCskTjcoagaiaSrPcUminEmo7lvAHCMMXaCMTYC4OcA8qQgY+x1xhgPxN0LwH4lBg3oCkWyG9Qr\nV66U/v6DDz6opR12OMU0qyZi6MQqHGVS83XUjrEL4eMlcFXC43TYdHVHt1RWVkoHCIj9kEwmC4Si\ndTk3vuIPMLEcnRPFWt6RKwiqC2DHMZ7cD1GEncoI9w4A4tIpZye3OfENAL+220FE3ySiXiLq9buo\ntNMKPqrwQS2eT0UAxK2YkSqq03RrZIHq4tZOPPzwwwXbvMLEZJ2jQYWxtb6PdXwEreM+MjLia1Wi\nbDZbIBRXr17tWPc9qjK2UwEZn1hbW5tnbL1s4cMgaJVQRPQHANYB2G63nzH2LGNsHWNsnUp1OxGu\nMcp8XyYlWNR4VLSCKGPew0DVZhlWBccdO3YUvCj37dsHALj55pttv+PlsOJOXqvzS/WF5vUSiVMB\ntX379vmKYNHhEI8COyVAxKoQhJWMJzPT6e/vz4XgOuVV6K5/b4eMcD8HQJw7zp7clgcRbQHwdwDu\nZ4yFVsuVV4WU0fxVKys++eST0sc+88wzSueOI36dbLpx0ix3794t9X2rFsSFv3UtzFKJRIqSoJEd\nYZfq4GvH7tixw/U4rhBwZMJZoygz4uQjiEtVyHcALCKiLiKqBLANwE7xACJaA+DfMCHYI1vUcdas\nWa7aucxiEFavtYz2Hrc1WP2iK+IgLG6//Xap42Rr7YRd7TBKnBzNMr4QHVRWVqKiosJ3WQBZwoxP\nD7vtQL4CFXU1V0/hzhjLAvg2gJcAHAbwHGPsIBH9gIjunzxsO4A6AP9DRPuIaKfD6QLDQyGTyST6\n+vryhHEqlUIikcg5m0SnkpOgP3jwYN7fTzzxhKvtfWRkREnDn0rwEEldU31euMsPcS/+FhQnp7su\nX4gXIyMjys5QJ9wWiw5bAMtE6/jx3XCbu6hAiVp8LEIhAYAxtosxtpgxtoAx9sTktu8xxnZOft7C\nGEszxlZP/rvf/Yz+4XZOO81ldHS0YLUYPmV3coDaaXPbt2/HG2+8USDk9+zZg6effjpvWxxW64kL\nPAIiyiQOpx+nk1CIY2nWqU4UJgo7qqurPV9QK1asyJMD1dXVWL9+fe5v0XEqlk5weiFweRFF7ShT\nfiAAd999N1566aViNyNSenp68mLUOzo6cO7cOSxatAgff/yx9HlSqVTRiqMVk6amplBzNUqNRCJR\nMtE8M2fOtK0aOXv27IJQVGBCkLslNpnyAxZknKTccaZahEi1zK6fzEndmqPVjhf28nzW5AtuE1UR\n7EAwzSWdTiORSEhHRMjaoaNYoCSKtTNLiTgKdi4HrNnKTuWA7QQ7MGFdcDJRRpFlXnLCnYfHOdXV\n6O7uzoUrZTKZnICXSRbhAke23oQfe6DoI9ARK2/1xodde7xYSS8i58+fx/j4uHSBLyc79OzZs/Om\nz0EWKJFNgOMRPEHr5bghxru7jYegvhHdFUzjYjLjcoCXVuCLmDj1l5NSmEgkCkyUX/nKVwAAN910\nk5a2ulFywp2HQjrV1bD+kLkjVSVZRLTDW4sc6cSv1nLXXXfZbq+trZVasScIpR7fL1JTU6OtlswH\nH3ygdHxY9XIA5MW7u9mUrYJHVdnQPdacItUaGxu1F0kD5LXnd955B4CzL8lpFmoXjfbiiy8CAH7z\nm99IXTsIJSfcveAlPoHfL/7LcdOunKbkdnVUirVKEjBxf04DY3BwMM/5Y1eRLwqCaGAVFRWRxd8f\nPXp0ytq/7aKJgppIdGryYub4wMCAUnlpt3EvzmSuXLmCqqqqgnVxxQTJ5cuX+05WGxsbc6zbHoW/\nqeSE+x133OG6P5vN5qbX1sHqpF0lEgnbVPKWlhbbt3sUD8YptV0mDp/7HD7//HNtPzivBCDRvxGk\n/sfY2Jh0/L3fLOepgpsp0q2P/YaRBtHkeVE8jpMdWwa3F7Z1JjM8PFywLq6YIHngwIFAL71irmxW\ncsJdtJHqsNFNnz4d1dXVtmaby5cv+67eFtRx5ldbmDFjRp5d3O0HpzIN97K181R9vtKMLmpqahyF\nlFOWchxWgOJYtUI3dM8InQR4Mpl0XQjayZfEK1PKxH3LJuz09PQglUr5Loon0xa3Gux+HZsyv51i\nj8OSE+7c5g741xDFB3rhwoVQVhQSzUGy6DCjqDh5ubYUVLvv6enJPYu9e/cGOpeVoaEh5eJaqn4B\nnUs3WrFqhW7onhE6jetsNiudrStGX/HKlDL2b9nSvIcOHQp03zKh3G412P0qbzLafCaTKeraCyUn\n3EV48R03LcQO6wPt7u4umAXYvZkbGhrQ3NyMhoaGUBZ10GX/lbVZc8EjavdElOtXN6Enmo2stdk3\nb94sdf2wIntUQxr9VGaUQbYf4oxd8TS/TnU+ZsJeGLyxsRHbtm0DgEBrH4uLnfjF6eURRUhsySUx\nvfXWW3jllVccEwecqKio8MxGk137MGxqampCj3pZtmxZQemFdDqNixcvhhZ73NXVJV0HJiqSyWTZ\n1AiXRWVN2VKlvb1dqrigG5WVlaFV/rzvvvuwbp1nHpItZZvE9OabbwJQd7jI1MGwCnY3DSPMN6+T\nYHfTdFUjTKyCnZ8jzKQSN8HOy0TI2GpVZ2puTDXBDsjXw48r06ZN8wxjtQp2P873MMIvOb/97W9D\nOzen5IS7+CaV/ZH7TRYSZzVWZ5eTTd0uQzToog4ctxeUjgqPhw8fjnwN0a6uLiSTyZzgl8k69ePP\nCIqK8161tHAxQ2ujQlSGvOoxufmA6uvrcf36dQwNDXmu3yriR4sPcym8MPx8VkpOuPOFdjs6OhzT\nga3o0EZlnD5tbW22WpGfqR3XNJYuXar8XTeNo7a21tXJ09zcHGmI4SeffJKnPff29jq2j2deFmOR\nDLGN4gyqrq6uIIJLNsOZw8dW2NnFYZemcEN8OV67ds1xjDY1NbmaJK1RW16rdonMmDFDe1Yt4G8W\nL1vOOgglJ9z5oPBaMLgYiBECTssByoaIcRORuIjzhg0bpJw8bj+OwcFBV83p5MmT0lqO07JuQeGa\nr9UhqbLSUJi1s8fGxpBOp1FdXY1MJpM3w8tmszhw4ICv8+pos5t5TiZ+XVd0h3WmY/17y5YttrMh\nlVk2ESk5dz/77DMtvixr1rrdTNKrHEUYLxkrJedQFatCLly4sKiLRrsR1KHT1tambaX3JUuW4MiR\nI7m/dVfis6vwWFVVVRaLmnDH/dKlS/NetFGgqw955c4ormWFV8Gsr6/XXpdINahCN7xksN+ZpKkK\n6QIX7EHrv+go4JVKpbBx40YAwPr1631H3SxbtgyAfJywDFywc3u2rkXGOaJg51PUqAS7DlODmxa1\nZMkSACgQ7LNmzQp9/VGxD4OED3oJ9s7OztCeFw/v9SvY3e5bRrCHqSGnUilXwc5DihctWhRaG9wo\nOeHOV2ISO8yu/osKTlqsSlTG+Ph4zgHzxRdfeEbncPuq1bxkF8WiC15F8fTp04HP5WTrdHN29vT0\nSAljFdOAjsgPt6n6q6++aru9r68v0rBZxljuxd/Z2Vlg5+W+KI6K/V62uqYd6XS6wEfDSyzLxpi7\nKVfieKmoqFB+oYYZUiw+f/5sRHj+h91vwm31KV2UnHDn9WFU64f7QaUEbEdHRy4p6MiRI7Z2bXFg\ncuHPs/44TppKOp1WdtSFARGhuro6tyCxNdLDGhkkCqFDhw5JCeMwoxS8mD9/fu4ZzJkzJzTt3Oqv\nkImY4S/lkydPFkRH7dmzJ+9v1Wxnv2Pr/PnzBeZHPlPkvgevtYxF5co6/q9du4bGxkZUVVUhnU77\nnqk5RWGtXr1aKRLKqRCYm1JmF8mmWkXUD1LCnYjuIaIjRHSMiL5rs7+KiH4xuf9tIurU3VCOilNN\nBjcnk4pd2qoN201DZTQ9Jx/I+fPnY7HAM2Msb3UZq62dT1N5dmsxwhaDcOLEidwzOHPmTGjauXUc\ny0RjiWPKy1em2u+6xlZ9fX1B21RMPnb3NTAwgOHhYfT19fn+/TvNTqzRWl7oKgQWVla0iKdwJ6IK\nAD8CcC+AHgBfJ6Iey2HfAHCFMbYQwD8B+EfdDeWohD7ZYZ0+6Uzo4FqFaqx4R0eHcmx0EHRpo262\nRKfBW6wyxFasGmCU/W/H1772tcDnEM02dnCt0218emnZXly9ejXvN6ojqi3MEtBBbPKq1TPF+4jL\nYh0bABxjjJ1gjI0A+DmAByzHPADgPyY/7wBwB4VUQCJouJh1+mTVsINknvIXhUpyBeC8fJcbQZzA\nuuKd3UxjTgLErn5OW1ubrxjvuro6qfHAs19FrC/1MEMnRZye23PPPRf43HxsO2mpXOsUx2dbW1ve\ns9LhWBXPbzU7+kFHgp4TqnVyxL5SHbPifcRlmb0OAKIadnZym+0xjLEsgAEAoeTni1UhdWB1fBbD\njDA+Pq4cTRAklFG3acsOlRdcf3+/VHkIK5lMRiqqSKaeTVQ1b+K2Zmh/f7+yMuJF3OoHuaH6uxP7\nKoiJJoqw2kgdqkT0TSLqJaJevzHg3GGnM5uvXNK/i5mB6Iabw0p1gqcjbNXvtaNCdZ0CN9NCMe9R\n55qoqvfh9ltQlR0667LzJMwoVhuT6f1zAMTar7Mnt9kdc5aIkgAaARSowIyxZwE8C0wkMflp8GOP\nPQbGGIgo97/lGo7bZI7nDh2nweT32irHyF5H3Od2H7LXdEPmPmTvh7dLpq2y9yhzbdVtdtdWfR5+\nx5zq85G5D6++d7t20PFibYfK+XWM7TCeX9D7CBsZNegdAIuIqIuIKgFsA7DTcsxOAH88+flhAK+x\nEFNfeQfZdZTbNpnjicj1Afi9tsoxstcR97ndh+w13ZC5D7d9YltU2ip7jzLXVt1md23V5+F3zAUR\nAk7n8up7t2sHHS/WdqicX8fYDuP5Bb2PsJEqP0BEWwE8DaACwE8ZY08Q0Q8A9DLGdhLRNAA/A7AG\nwGUA2xhjzsufTJzzIoBTPtvdBkBfCmd5YPrEHtMvhZg+KaSU+mQeY8yzul/RassEgYh6ZWorTCVM\nn9hj+qUQ0yeFlGOflFyGqsFgMBi8McLdYDAYypBSFe7PFrsBMcT0iT2mXwoxfVJI2fVJSdrcDQaD\nweBOqWruBoPBYHAh1sKdYlSNMi5I9MkjRHSRiPZN/vvTYrQzSojop0R0gYhs17ejCZ6Z7LMPiOiG\nqNsYNRJ9chsRDQjj5HtRtzFqiGgOEb1ORIeI6CARfcfmmPIZK4yxWP7DREz9cQDzAVQC2A+gx3LM\nnwP48eTnbQB+Uex2x6BPHgHwr8Vua8T9ciuAGwAccNi/FcCvARCAGwG8Xew2x6BPbgPwf8VuZ8R9\nMhPADZOf6wEctfn9lM1YibPmHqtqlDFBpk+mHIyxNzGRPOfEAwD+k02wF0ATEckvs1WCSPTJlIMx\n9ilj7L3Jz1cBHEZhEcSyGStxFu6xqkYZE2T6BAAempxS7iCiOTb7pxqy/TbV2ERE+4no10RUuE5c\nGTNpwl0D4G3LrrIZK3EW7gZ/vAigkzG2EsDL+P3MxmAQeQ8TaeyrAPwLgBeK3J7IIKI6AM8D+EvG\nWHQL4UZMnIW7SjVKkEs1yjLCs08YY5cYY3zFhZ8AWBtR2+KMzFiaUjDGvmCMZSY/7wKQIqJoViwp\nIkSUwoRg/y/G2P/aHFI2YyXOwj121ShjgGefWOyD92PCrjjV2QngjyYjIW4EMMAYU1/+qowgohnc\nP0VEGzAhC8pZMcLk/f47gMOMsR86HFY2Y0VfNX3NMMayRPRtAC/h99UoD4rVKDHxoH5GRMcwWY2y\neC0OH8k++Qsiuh9AFhN98kjRGhwRRPTfmIj+aCOiswAeB5ACAMbYjwHswkQUxDEAQwD+pDgtjQ6J\nPnkYwJ8RURbANUxUci1nxQgANgP4QwAfEtG+yW1/C2AuUH5jxWSoGgwGQxkSZ7OMwWAwGHxihLvB\nYDCUIUa4GwwGQxlihLvBYDCUIUa4GwwGQxlihLvBYDCUIUa4GwwGQxlihLvBYDCUIf8Pkf64pudq\ncAEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2f8b343f98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from fenics import *\n",
    "%matplotlib inline\n",
    "\n",
    "from mshr import *\n",
    "\n",
    "channel = Rectangle(Point(0, 0), Point(2.2, 0.41))\n",
    "cylinder = Circle(Point(0.2, 0.2), 0.05, 16)\n",
    "domain = channel - cylinder\n",
    "\n",
    "mesh = generate_mesh(domain, 64)\n",
    "plot(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to define the finite element function space. For this\n",
    "problem, we need to define several spaces. The first space we create\n",
    "is the space for the velocity field $w$ from the Navier--Stokes\n",
    "simulation. We call this space $W$ and define the space by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define function space for velocity\n",
    "W = VectorFunctionSpace(mesh, 'P', 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important that this space is exactly the same as the space we\n",
    "used for the velocity field in the Navier-Stokes solver. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the three concentrations $u_1$, $u_2$, and $u_3$, we want to\n",
    "create a *mixed space* using `MixedElement` with three scalar functions representing the full system $(u_1, u_2, u_3)$ as a single entity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define function space for system of concentrations\n",
    "P1 = FiniteElement('P', triangle, 1)\n",
    "element = MixedElement([P1, P1, P1])\n",
    "V = FunctionSpace(mesh, element)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: We can also define the mixed finite elements with:\n",
    "```\n",
    "element = VectorElement('P', triangle, 1, dim=3)\n",
    "V = FunctionSpace(mesh, element)\n",
    "```\n",
    "\n",
    "Note that for Navier-Stoke's equation, we used the product of elements to define the mixed Taylor-Hood element:\n",
    "```\n",
    "P2 = VectorElement('P', triangle, 2)\n",
    "P1 = FiniteElement('P', triangle, 1)\n",
    "TH = P2 * P1\n",
    "```\n",
    "However, this syntax does not work well with three or more elements because `element = P1 * P1 * P1` is interpreted as `element = (P1 * P1) * P1`, and the mixed element would be consisting of two subsystems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the space has been created, we need to define our test functions\n",
    "and finite element functions. Test functions for a mixed function\n",
    "space can be created by replacing `TestFunction` by `TestFunctions`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define test functions\n",
    "v_1, v_2, v_3 = TestFunctions(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the problem is nonlinear, we need to work with functions rather\n",
    "than trial functions for the unknowns. This can be done by using the\n",
    "corresponding `Functions` construction in FEniCS. However, as we will\n",
    "need to access the `Function` for the entire system itself, we first\n",
    "need to create that function and then access its components:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define functions for velocity and concentrations\n",
    "w = Function(W)\n",
    "# w = Constant((1,0))\n",
    "u = Function(V)\n",
    "u_n = Function(V)\n",
    "\n",
    "# Split system functions to access components\n",
    "u_1, u_2, u_3 = split(u)\n",
    "u_n1, u_n2, u_n3 = split(u_n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These functions will be used to represent the unknowns $u_1$, $u_2$, and $u_3$\n",
    "at the new time level $n+1$. The corresponding values at the previous\n",
    "time level $n$ are denoted by `u_n1`, `u_n2`, and `u_n3` in our program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define source terms\n",
    "f_1 = Expression('pow(x[0]-0.1,2)+pow(x[1]-0.1,2)<0.05*0.05 ? 0.1 : 0',\n",
    "                 degree=1)\n",
    "f_2 = Expression('pow(x[0]-0.1,2)+pow(x[1]-0.3,2)<0.05*0.05 ? 0.1 : 0',\n",
    "                 degree=1)\n",
    "f_3 = Constant(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define some parameters and expressions used in variational forms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "T = 0.5            # final time\n",
    "num_steps = 500    # number of time steps\n",
    "dt = T / num_steps # time step size\n",
    "eps = 0.01         # diffusion coefficient\n",
    "K = 10.0           # reaction rate\n",
    "\n",
    "# Define expressions used in variational forms\n",
    "k = Constant(dt)\n",
    "K = Constant(K)\n",
    "eps = Constant(eps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now all functions and test functions have been defined. We can\n",
    "express the nonlinear variational problem\n",
    "(\\ref{ftut1-reactionsystem-varproblem}):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "F = ((u_1 - u_n1) / k)*v_1*dx + dot(w, grad(u_1))*v_1*dx \\\n",
    "  + eps*dot(grad(u_1), grad(v_1))*dx + K*u_1*u_2*v_1*dx  \\\n",
    "  + ((u_2 - u_n2) / k)*v_2*dx + dot(w, grad(u_2))*v_2*dx \\\n",
    "  + eps*dot(grad(u_2), grad(v_2))*dx + K*u_1*u_2*v_2*dx  \\\n",
    "  + ((u_3 - u_n3) / k)*v_3*dx + dot(w, grad(u_3))*v_3*dx \\\n",
    "  + eps*dot(grad(u_3), grad(v_3))*dx - K*u_1*u_2*v_3*dx + K*u_3*v_3*dx \\\n",
    "  - f_1*v_1*dx - f_2*v_2*dx - f_3*v_3*dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One could write out the velocity field into a `TimeSeries` and then read it back in, such as \n",
    "```\n",
    "timeseries_w = TimeSeries('navier_stokes_cylinder/velocity_series')\n",
    "```\n",
    "However, file I/O is too expensive compared to computation. Hence, we solve the Navier-Stokes equation directly.\n",
    "\n",
    "To avoid name conflicts between the Navier-Stokes solver and the advection-diffusion-reaction solver, we encapsulate the Navier-Stokes solver into a class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class NavierStokesSolver:\n",
    "    # Define function spaces\n",
    "    inflow   = 'near(x[0], 0)'\n",
    "    outflow  = 'near(x[0], 2.2)'\n",
    "    walls    = 'near(x[1], 0) || near(x[1], 0.41)'\n",
    "    cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3'\n",
    "\n",
    "    # Define inflow profile\n",
    "    inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')\n",
    "\n",
    "    mu = 0.001         # dynamic viscosity\n",
    "    rho = 1            # density\n",
    "\n",
    "    f  = Constant((0, 0))\n",
    "    mu = Constant(mu)\n",
    "\n",
    "    def epsilon(u):\n",
    "        return sym(nabla_grad(u))\n",
    "\n",
    "    # Define stress tensor\n",
    "    def sigma(self, u, p):\n",
    "        return 2*self.mu*NavierStokesSolver.epsilon(u) - p*Identity(len(u))\n",
    "\n",
    "    # Define symmetric gradient\n",
    "    def __init__(self, mesh, dt):\n",
    "        self.V = VectorFunctionSpace(mesh, 'P', 2)\n",
    "        self.Q = FunctionSpace(mesh, 'P', 1)\n",
    "\n",
    "        # Define boundary conditions\n",
    "        bcu_inflow = DirichletBC(self.V, Expression(self.inflow_profile, degree=2), self.inflow)\n",
    "        bcu_walls = DirichletBC(self.V, Constant((0, 0)), self.walls)\n",
    "        bcu_cylinder = DirichletBC(self.V, Constant((0, 0)), self.cylinder)\n",
    "        bcp_outflow = DirichletBC(self.Q, Constant(0), self.outflow)\n",
    "        self.bcu = [bcu_inflow, bcu_walls, bcu_cylinder]\n",
    "        self.bcp = [bcp_outflow]\n",
    "\n",
    "        # Define trial and test functions\n",
    "        self.v = TestFunction(self.V)\n",
    "        self.q = TestFunction(self.Q)\n",
    "\n",
    "        self.u = TrialFunction(self.V)\n",
    "        self.p = TrialFunction(self.Q)\n",
    "\n",
    "        # Define functions for solutions at previous and current time steps\n",
    "        self.u_k = Function(self.V)\n",
    "        self.u_  = Function(self.V)\n",
    "\n",
    "        self.p_k = Function(self.Q)\n",
    "        self.p_  = Function(self.Q)\n",
    "\n",
    "        # Define expressions used in variational forms\n",
    "        self.DT = Constant(dt)\n",
    "\n",
    "        self.u_mid  = 0.5*(self.u_k + self.u)\n",
    "        self.n  = FacetNormal(mesh)\n",
    "\n",
    "        # Define variational problem for step 1\n",
    "        self.F1 = self.rho*dot((self.u - self.u_k) / self.DT, self.v)*dx \\\n",
    "           + self.rho*dot(dot(self.u_k, nabla_grad(self.u_k)), self.v)*dx \\\n",
    "           + inner(self.sigma(self.u_mid, self.p_k), NavierStokesSolver.epsilon(self.v))*dx \\\n",
    "           + dot(self.p_k*self.n, self.v)*ds - dot(self.mu*nabla_grad(self.u_mid)*self.n, self.v)*ds \\\n",
    "           - dot(self.f, self.v)*dx\n",
    "\n",
    "        self.a1 = lhs(self.F1)\n",
    "        self.L1 = rhs(self.F1)\n",
    "        self.A1 = assemble(self.a1)\n",
    "        [bc.apply(self.A1) for bc in self.bcu]\n",
    "\n",
    "        # Define variational problem for step 2\n",
    "        self.a2 = dot(nabla_grad(self.p), nabla_grad(self.q))*dx\n",
    "        self.L2 = dot(nabla_grad(self.p_k), nabla_grad(self.q))*dx - (1/self.DT)*div(self.u_)*self.q*dx\n",
    "        self.A2 = assemble(self.a2)\n",
    "        [bc.apply(self.A2) for bc in self.bcp]\n",
    "\n",
    "        # Define variational problem for step 3\n",
    "        self.a3 = dot(self.u, self.v)*dx\n",
    "        self.L3 = dot(self.u_, self.v)*dx - self.DT*dot(nabla_grad(self.p_ - self.p_k), self.v)*dx\n",
    "        self.A3 = assemble(self.a3)\n",
    "        \n",
    "    def advance(self):\n",
    "        # Update current time\n",
    "        # Step 1: Tentative velocity step\n",
    "        b1 = assemble(self.L1)\n",
    "        [bc.apply(b1) for bc in self.bcu]\n",
    "        solve(self.A1, self.u_.vector(), b1, 'bicgstab', 'ilu')\n",
    "\n",
    "        # Step 2: Pressure correction step\n",
    "        b2 = assemble(self.L2)\n",
    "        [bc.apply(b2) for bc in self.bcp]\n",
    "        solve(self.A2, self.p_.vector(), b2, 'bicgstab', 'ilu')\n",
    "\n",
    "        # Step 3: Velocity correction step\n",
    "        b3 = assemble(self.L3)\n",
    "        solve(self.A3, self.u_.vector(), b3, 'cg', 'sor')\n",
    "        \n",
    "        # Update previous solution\n",
    "        self.u_k.assign(self.u_)\n",
    "        self.p_k.assign(self.p_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create VTK files for visualization output\n",
    "vtkfile_u_1 = File('reaction_system/u_1.pvd')\n",
    "vtkfile_u_2 = File('reaction_system/u_2.pvd')\n",
    "vtkfile_u_3 = File('reaction_system/u_3.pvd')\n",
    "\n",
    "# Create progress bar\n",
    "progress = Progress('Time-stepping')\n",
    "set_log_level(PROGRESS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The time-stepping simply consists of solving this variational problem\n",
    "in each time step by a call to the `solve` function.\n",
    "In each time step, we first read the current value for the velocity\n",
    "field from the time series we have previously stored. We then solve\n",
    "the nonlinear system, and assign the computed values to the left-hand\n",
    "side values for the next time interval. When retrieving values from a\n",
    "`TimeSeries`, the values will by default be interpolated (linearly) to\n",
    "the given time `t` if the time does not exactly match a sample in the\n",
    "series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u max:  0.000102884795044\n",
      "u max:  0.000622918813026\n",
      "u max:  0.00113870053516\n",
      "u max:  0.00164590541457\n",
      "u max:  0.00212818577667\n",
      "u max:  0.00257295714994\n",
      "u max:  0.00297380436068\n",
      "u max:  0.00332800481701\n",
      "u max:  0.0036619628641\n",
      "u max:  0.00398779981288\n",
      "u max:  0.00427665295716\n",
      "u max:  0.00452881375789\n",
      "u max:  0.00474523501681\n",
      "u max:  0.00492735408589\n",
      "u max:  0.00507703708594\n",
      "u max:  0.00519655891911\n",
      "u max:  0.00534775055116\n",
      "u max:  0.00548224944723\n",
      "u max:  0.00559843226931\n",
      "u max:  0.00569764378726\n",
      "u max:  0.00578085980865\n",
      "u max:  0.00584882756265\n",
      "u max:  0.00590206962203\n",
      "u max:  0.0059409527439\n",
      "u max:  0.00596623186284\n",
      "u max:  0.00597962627691\n",
      "u max:  0.00598372775757\n",
      "u max:  0.00598138331324\n",
      "u max:  0.0059751073556\n",
      "u max:  0.00596679047832\n",
      "u max:  0.00595771662041\n",
      "u max:  0.00594874369799\n",
      "u max:  0.00594045956532\n",
      "u max:  0.00593325515397\n",
      "u max:  0.00592735924945\n",
      "u max:  0.00592287123065\n",
      "u max:  0.00591979369416\n",
      "u max:  0.0059180594189\n",
      "u max:  0.00591755248402\n",
      "u max:  0.00591812377219\n",
      "u max:  0.00591960462297\n",
      "u max:  0.00592182262882\n",
      "u max:  0.00592461524073\n",
      "u max:  0.00592783603494\n",
      "u max:  0.0059313555554\n",
      "u max:  0.00593506100857\n",
      "u max:  0.00593885592105\n",
      "u max:  0.00594265934764\n",
      "u max:  0.00594640500035\n",
      "u max:  0.0059500405091\n",
      "u max:  0.00595352655672\n",
      "u max:  0.00595683548114\n",
      "u max:  0.00595994952699\n",
      "u max:  0.00596285921274\n",
      "u max:  0.00596556173652\n",
      "u max:  0.00596805936626\n",
      "u max:  0.00597035810908\n",
      "u max:  0.00597246676366\n",
      "u max:  0.00597439614858\n",
      "u max:  0.00597615845571\n",
      "u max:  0.00597776670834\n",
      "u max:  0.00597923431638\n",
      "u max:  0.00598057472804\n",
      "u max:  0.00598180116252\n",
      "u max:  0.00598292638088\n",
      "u max:  0.00598396254329\n",
      "u max:  0.00598492110836\n",
      "u max:  0.00598581275664\n",
      "u max:  0.00598664739526\n",
      "u max:  0.00598743417151\n",
      "u max:  0.00598818148231\n",
      "u max:  0.00598889700579\n",
      "u max:  0.00598958775769\n",
      "u max:  0.00599026013218\n",
      "u max:  0.00599091995438\n",
      "u max:  0.0059915725329\n",
      "u max:  0.00599222270723\n",
      "u max:  0.00599287488682\n",
      "u max:  0.00599353310137\n",
      "u max:  0.00599420104643\n",
      "u max:  0.00599488210589\n",
      "u max:  0.00599557939974\n",
      "u max:  0.0059962957981\n",
      "u max:  0.00599703393962\n",
      "u max:  0.00599779625158\n",
      "u max:  0.00599858496138\n",
      "u max:  0.0059994021033\n",
      "u max:  0.00600024952458\n",
      "u max:  0.00600112889\n",
      "u max:  0.00600204168514\n",
      "u max:  0.00600298920756\n",
      "u max:  0.00600397257801\n",
      "u max:  0.00600499272572\n",
      "u max:  0.00600605039145\n",
      "u max:  0.00600714612159\n",
      "u max:  0.00600828025481\n",
      "u max:  0.00600945293116\n",
      "u max:  0.00601066407143\n",
      "u max:  0.00601191337938\n",
      "u max:  0.00601320032301\n"
     ]
    }
   ],
   "source": [
    "# Time-stepping\n",
    "t = 0\n",
    "out_interval = num_steps / 100;\n",
    "\n",
    "nss = NavierStokesSolver(mesh, dt)\n",
    "\n",
    "for k in range(num_steps):\n",
    "    # Update current time\n",
    "    t += dt\n",
    "\n",
    "    # Advance the Navier-Stokes solver in time\n",
    "    nss.advance()\n",
    "\n",
    "    # Copy velocities from Navier-Stokes solver to a-d-r solver\n",
    "    w.assign(nss.u_k)\n",
    "\n",
    "    # Solve variational problem for time step\n",
    "    solve(F == 0, u)\n",
    "\n",
    "    if k%out_interval ==0 or k==num_steps:\n",
    "        # Save solution to file (VTK)\n",
    "        _u_1, _u_2, _u_3 = u.split()\n",
    "        vtkfile_u_1 << (_u_1, t)\n",
    "        vtkfile_u_2 << (_u_2, t)\n",
    "        vtkfile_u_3 << (_u_3, t)\n",
    "\n",
    "        print('u max: ', u.vector().array().max())\n",
    "        \n",
    "    # Update previous solution\n",
    "    u_n.assign(u)\n",
    "\n",
    "    # Update progress bar\n",
    "    progress.update(t / T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we comment on three important techniques that are very useful\n",
    "when working with systems of PDEs: setting initial conditions, setting\n",
    "boundary conditions, and extracting components of the system for\n",
    "plotting or postprocessing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting initial conditions for mixed systems ##\n",
    "\n",
    "In our example, we did not need to worry about setting an initial\n",
    "condition, since we start with $u_1 = u_2 = u_3 = 0$. This happens\n",
    "automatically in the code when we set `u_n = Function(V)`. This\n",
    "creates a `Function` for the whole system and all degrees of freedom\n",
    "are set to zero.\n",
    "\n",
    "If we want to set initial conditions for the components of the system\n",
    "separately, the easiest solution is to define the initial conditions\n",
    "as a vector-valued `Expression` and then project (or interpolate) this\n",
    "to the `Function` representing the whole system. For example,\n",
    "\n",
    "```\n",
    "u_0 = Expression(('sin(x[0])', 'cos(x[0]*x[1])', 'exp(x[1])'), degree=1)\n",
    "u_n = project(u_0, V)\n",
    "```\n",
    "This defines $u_1$, $u_2$, and $u_2$ to be the projections of $\\sin\n",
    "x$, $\\cos (xy)$, and $\\exp(y)$, respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting boundary conditions for mixed systems ##\n",
    "\n",
    "In our example, we also did not need to worry about setting boundary\n",
    "conditions since we used a natural Neumann condition. If we want to set\n",
    "Dirichlet conditions for individual components of the system, this can\n",
    "be done as usual by the class `DirichletBC`, but we must specify for\n",
    "which subsystem we set the boundary condition. For example, to specify\n",
    "that $u_2$ should be equal to $xy$ on the boundary defined by\n",
    "`boundary`, we do\n",
    "\n",
    "```\n",
    "u_D = Expression('x[0]*x[1]', degree=1)\n",
    "bc = DirichletBC(V.sub(1), u_D, boundary)\n",
    "```\n",
    "The object `bc` or a list of such objects containing different\n",
    "boundary conditions, can then be passed to the `solve` function as usual.\n",
    "Note that numbering starts at $0$ in FEniCS so the subspace\n",
    "corresponding to $u_2$ is `V.sub(1)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing components of mixed systems ##\n",
    "\n",
    "If `u` is a `Function` defined on a mixed function space in FEniCS,\n",
    "there are several ways in which `u` can be *split* into components.\n",
    "Above we already saw an example of the first of these:\n",
    "\n",
    "```\n",
    "u_1, u_2, u_3 = split(u)\n",
    "```\n",
    "This extracts the components of `u` as *symbols* that can be used in a\n",
    "variational problem. The above statement is in fact equivalent to\n",
    "\n",
    "```\n",
    "u_1 = u[0]\n",
    "u_2 = u[1]\n",
    "u_3 = u[2]\n",
    "```\n",
    "Note that `u[0]` is not really a `Function` object, but merely a\n",
    "symbolic expression, just like `grad(u)` in FEniCS is a symbolic\n",
    "expression and not a `Function` representing the gradient.  This means\n",
    "that `u_1`, `u_2`, `u_3` can be used in a variational problem, but\n",
    "cannot be used for plotting or postprocessing.\n",
    "\n",
    "To access the components of `u` for plotting and saving the solution\n",
    "to file, we need to use a different variant of the `split` function:\n",
    "\n",
    "```\n",
    "u_1_, u_2_, u_3_ = u.split()\n",
    "```\n",
    "This returns three subfunctions as actual objects with access to the\n",
    "common underlying data stored in `u`, which makes plotting and saving\n",
    "to file possible. Alternatively, we can do\n",
    "\n",
    "```\n",
    "u_1_, u_2_, u_3_ = u.split(deepcopy=True)\n",
    "```\n",
    "which will create `u_1_`, `u_2_`, and `u_3_` as stand-alone `Function`\n",
    "objects, each holding a copy of the subfunction data extracted from\n",
    "`u`. This is useful in many situations but is not necessary for\n",
    "plotting and saving solutions to file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Acknowledgment**: This notebook was adapted from [The FEniCS Tutorial Volume I](https://fenicsproject.org/pub/tutorial/sphinx1/) by Hans Petter Langtangen and Anders Logg, released under CC Attribution 4.0 license."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}