notebooks/navier_stokes_cylinder.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Case Study: Flow past a cylinder\n",
    "\n",
    "We solve a relatively more challenging problem: flow past a \n",
    "circular cylinder. The geometry is as shown below:\n",
    "![Geoetry](http://www.featflow.de/media/dfg_bench2_2d/geometry.png)\n",
    "\n",
    "This problem was taken from the [\"FEATFLOW/1995-DFG benchmark suite\"](\n",
    "http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html).\n",
    "Notice the slightly perturbed and unsymmetric geometry.\n",
    "\n",
    "The kinematic viscosity is\n",
    "given by $\\nu = 0.001 = \\mu/\\rho$ and the inflow velocity profile is\n",
    "specified as\n",
    "$$  u(x, y, t) = \\left(1.5 \\cdot \\frac{4y(0.41 - y)}{0.41^2}, 0\\right),$$\n",
    "\n",
    "which has a maximum magnitude of $1.5$ at $y = 0.41/2$. We do not\n",
    "use any scaling for this problem since all exact parameters are known."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FEniCS implementation\n",
    "\n",
    "### Generate mesh and function space\n",
    "FEniCS supports the creation of relatively complex meshes via a technique\n",
    "called *constructive solid geometry* (CSG), which lets us define\n",
    "geometries in terms of simple shapes (primitives) and set operations:\n",
    "union, intersection, and set difference. The set operations are\n",
    "encoded in FEniCS using the operators `+` (union), `*` (intersection),\n",
    "and `-` (set difference). To access the CSG functionality in FEniCS,\n",
    "one must import the FEniCS module `mshr` which provides the\n",
    "extended meshing functionality of FEniCS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fenics import *\n",
    "from mshr import *\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The geometry for the cylinder flow test problem can be defined easily\n",
    "by first defining the rectangular channel and then subtracting the\n",
    "circle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in the above, we passed 16 for the third argument when calling `Circle` to request explicitly using 16 segments to represent the circle. This is necessary to circumvent a bug in `mshr` as of release 2016.02, which would generate a triangle to represent the circle if the third argument was missing.\n",
    "\n",
    "We may then generate the mesh by calling the function `generate_mesh`, plot the mesh, and save it to a file for inspection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "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 0x7ff6ecc71198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mesh = generate_mesh(domain, 64)\n",
    "plot(mesh)\n",
    "File('navier_stokes_cylinder/cylinder.pvd') << mesh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the argument `64` indicates that we want to resolve the geometry\n",
    "with 64 cells across its diameter (the channel length)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define the function spaces based on the mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define function spaces\n",
    "V = VectorFunctionSpace(mesh, 'P', 2)\n",
    "Q = FunctionSpace(mesh, 'P', 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Specify boundary conditions\n",
    "To solve the cylinder test problem, we only need to make a few minor\n",
    "changes to the code we wrote for the channel flow test\n",
    "case. Besides defining the new mesh, the only change we need to make\n",
    "is to modify the boundary conditions and the time step size. The\n",
    "boundaries are specified as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "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'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last line may seem cryptic before you catch the idea: we want to pick\n",
    "out all boundary points (`on_boundary`) that also lie within the 2D\n",
    "domain $[0.1,0.3]\\times [0.1,0.3]$. The only possible points are then the \n",
    "points on the circular boundary!\n",
    "\n",
    "We then collect the boundary conditions into lists."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "--- Instant: compiling ---\n",
      "--- Instant: compiling ---\n",
      "--- Instant: compiling ---\n",
      "--- Instant: compiling ---\n"
     ]
    }
   ],
   "source": [
    "# Define inflow profile\n",
    "inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0')\n",
    "\n",
    "# Define boundary conditions\n",
    "bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)\n",
    "bcu_walls = DirichletBC(V, Constant((0, 0)), walls)\n",
    "bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)\n",
    "bcp_outflow = DirichletBC(Q, Constant(0), outflow)\n",
    "bcu = [bcu_inflow, bcu_walls, bcu_cylinder]\n",
    "bcp = [bcp_outflow]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specify parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We specify density, viscosity, and timestep. We will set the time steps. \n",
    "relatively small time step to compute the solution (a time step that\n",
    "is too large will make the solution blow up)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu = 0.001         # dynamic viscosity\n",
    "rho = 1            # density\n",
    "\n",
    "T = 5.0            # final time\n",
    "num_steps = 5000   # number of time steps\n",
    "dt = T / num_steps # time step size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the variational problems\n",
    "\n",
    "The definition of the variational problem is the same as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n"
     ]
    }
   ],
   "source": [
    "# Define trial and test functions\n",
    "u = TrialFunction(V)\n",
    "v = TestFunction(V)\n",
    "p = TrialFunction(Q)\n",
    "q = TestFunction(Q)\n",
    "\n",
    "# Define symmetric gradient\n",
    "def epsilon(u):\n",
    "    return sym(nabla_grad(u))\n",
    "\n",
    "# Define stress tensor\n",
    "def sigma(u, p):\n",
    "    return 2*mu*epsilon(u) - p*Identity(len(u))\n",
    "\n",
    "# Define functions for solutions at previous and current time steps\n",
    "u_k = Function(V)\n",
    "u_  = Function(V)\n",
    "\n",
    "p_k = Function(Q)\n",
    "p_  = Function(Q)\n",
    "\n",
    "# Define expressions used in variational forms\n",
    "f  = Constant((0, 0))\n",
    "DT = Constant(dt)\n",
    "mu = Constant(mu)\n",
    "\n",
    "u_mid  = 0.5*(u_k + u)\n",
    "n  = FacetNormal(mesh)\n",
    "\n",
    "# Define variational problem for step 1\n",
    "F1 = rho*dot((u - u_k) / DT, v)*dx \\\n",
    "   + rho*dot(dot(u_k, nabla_grad(u_k)), v)*dx \\\n",
    "   + inner(sigma(u_mid, p_k), epsilon(v))*dx \\\n",
    "   + dot(p_k*n, v)*ds - dot(mu*nabla_grad(u_mid)*n, v)*ds \\\n",
    "   - dot(f, v)*dx\n",
    "\n",
    "a1 = lhs(F1)\n",
    "L1 = rhs(F1)\n",
    "A1 = assemble(a1)\n",
    "[bc.apply(A1) for bc in bcu]\n",
    "\n",
    "# Define variational problem for step 2\n",
    "a2 = dot(nabla_grad(p), nabla_grad(q))*dx\n",
    "L2 = dot(nabla_grad(p_k), nabla_grad(q))*dx - (1/DT)*div(u_)*q*dx\n",
    "A2 = assemble(a2)\n",
    "[bc.apply(A2) for bc in bcp]\n",
    "\n",
    "# Define variational problem for step 3\n",
    "a3 = dot(u, v)*dx\n",
    "L3 = dot(u_, v)*dx - DT*dot(nabla_grad(p_ - p_k), v)*dx\n",
    "A3 = assemble(a3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time-stepping\n",
    "\n",
    "We can now perform time-stepping. We will also add a *progress bar* so that \n",
    "we can follow the progress of our computation. This can be done as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "progress = Progress('Time-stepping')\n",
    "set_log_level(PROGRESS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Log levels and printing in FEniCS**. Notice the call to `set_log_level(PROGRESS)` which is essential to\n",
    "make FEniCS actually display the progress bar. FEniCS is actually\n",
    "quite informative about what is going on during a computation but the\n",
    "amount of information printed to screen depends on the current log\n",
    "level. Only messages with a priority higher than or equal to the\n",
    "current log level will be displayed. The predefined log levels in\n",
    "FEniCS are\n",
    "`DBG`,\n",
    "`TRACE`,\n",
    "`PROGRESS`,\n",
    "`INFO`,\n",
    "`WARNING`,\n",
    "`ERROR`, and\n",
    "`CRITICAL`. By default, the log level is set to `INFO` which means\n",
    "that messages at level `DBG`, `TRACE`, and `PROGRESS` will not be\n",
    "printed. Users may print messages using the FEniCS functions `info`,\n",
    "`warning`, and `error` which will print messages at the obvious log\n",
    "level (and in the case of `error` also throw an exception and\n",
    "exit). One may also use the call `log(level, message)` to print a\n",
    "message at a specific log level."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also make another change: since the system(s) of linear equations are \n",
    "significantly larger than for the simple channel flow test problem, we choose \n",
    "to use an iterative method instead of the default direct (sparse) solver used by\n",
    "FEniCS when calling `solve`. Efficient solution of linear systems\n",
    "arising from the discretization of PDEs requires the choice of both a\n",
    "good iterative (Krylov subspace) method and a good\n",
    "preconditioner. For this problem, we will simply use the biconjugate\n",
    "gradient stabilized method (BiCGSTAB) and the conjugate gradient method. This can be done by adding the\n",
    "keywords `bicgstab` or `cg` in the call to `solve`. You can determine which solvers\n",
    "are available by running the command\n",
    "```\n",
    "list_krylov_solver_methods()\n",
    "```\n",
    "Note that the output may not be displaced in the Jupyter Notebook but will be shown in the terminal where\n",
    "you started Jupyter Notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also specify suitable preconditioners to speed up the \n",
    "computations. You can determine which solvers\n",
    "are available by running the command\n",
    "```\n",
    "list_krylov_solver_preconditioners()\n",
    "```\n",
    "\n",
    "For example, we may use `petsc_amg` or `ilu`.\n",
    "```\n",
    "solve(A1, u1.vector(), b1, 'bicgstab', 'petsc_amg')\n",
    "solve(A2, p1.vector(), b2, 'bicgstab', 'petsc_amg')\n",
    "solve(A3, u1.vector(), b3, 'cg', 'sor')\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition, to be able to postprocess the computed solution in ParaView,\n",
    "we store the solution to a file in each time step. \n",
    "If you have HDF5 installed, you may use the `xdmf` file format, \n",
    "which has a few advantages. First, the storage is\n",
    "much more efficient, both in terms of speed and file sizes. Second,\n",
    "`.xdmf` files work in parallell, both for writing and reading\n",
    "(postprocessing). Much like `.pvd` files, the actual data will not be\n",
    "stored in the `.xdmf` file itself, but will instead be stored in a\n",
    "(single) separate data file with the suffix `.hdf5` which is an\n",
    "advanced file format designed for high-performance computing.\n",
    "However, if HDF5 is not available, we will still use `pvd` files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if 'HDF5File' in globals():\n",
    "    # Use XDMF files\n",
    "    xdmffile_u = XDMFFile('navier_stokes_cylinder/velocity.xdmf')\n",
    "    xdmffile_p = XDMFFile('navier_stokes_cylinder/pressure.xdmf')\n",
    "else:\n",
    "    vtkfile_u = File('navier_stokes_cylinder/velocity.pvd')\n",
    "    vtkfile_p = File('navier_stokes_cylinder/pressure.pvd')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, we may also store the solution using a FEniCS `TimeSeries`. This allows us\n",
    "to store the solution not for visualization, but for later reuse in a\n",
    "computation as we will see in the next section. Using a `TimeSeries`\n",
    "it is easy and efficient to read in solutions from certain points in\n",
    "time during a simulation. The `TimeSeries` class also uses the HDF5\n",
    "file format for efficient storage and access to data, so we will \n",
    "enable it only if HDF5 is available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if 'HDF5File' in globals():\n",
    "    # Create time series (for use in reaction_system.py)\n",
    "    timeseries_u = TimeSeries('navier_stokes_cylinder/velocity_series')\n",
    "    timeseries_p = TimeSeries('navier_stokes_cylinder/pressure_series')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is the complete code for time-stepping. We will print `max(u)` at each step. Please be advised that this example program is considerably more demanding than our previous examples in terms of CPU time and memory, but it should be possible to run the program on a reasonably modern laptop.\n",
    "\n",
    "Note that the progress log may not be displayed inside the Jupyter Notebook but in the terminal where JN was started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
      "u max: 1.92762010483\n",
      "u max: 2.01520053123\n",
      "u max: 2.22347152812\n",
      "u max: 2.21720346441\n",
      "u max: 2.19623569704\n",
      "u max: 2.17779478429\n",
      "u max: 2.16216516992\n",
      "u max: 2.15205491193\n",
      "u max: 2.14885171706\n",
      "u max: 2.15174188487\n",
      "u max: 2.15719648979\n",
      "u max: 2.15987132177\n",
      "u max: 2.15539813761\n",
      "u max: 2.14362397166\n",
      "u max: 2.1298462418\n",
      "u max: 2.12234017814\n",
      "u max: 2.12691842862\n",
      "u max: 2.14213090558\n",
      "u max: 2.15892963584\n",
      "u max: 2.16571572717\n",
      "u max: 2.15619093737\n",
      "u max: 2.13504630758\n",
      "u max: 2.11711902332\n",
      "u max: 2.1184334594\n",
      "u max: 2.14314658182\n",
      "u max: 2.17741113218\n",
      "u max: 2.19719858855\n",
      "u max: 2.18603504752\n",
      "u max: 2.15038678476\n",
      "u max: 2.13082446614\n",
      "u max: 2.13928040317\n",
      "u max: 2.17309063317\n",
      "u max: 2.22360359133\n",
      "u max: 2.2389764105\n",
      "u max: 2.20471231223\n",
      "u max: 2.14761011222\n",
      "u max: 2.16956636838\n",
      "u max: 2.15418814745\n",
      "u max: 2.21741381849\n",
      "u max: 2.26650722432\n",
      "u max: 2.25707086732\n",
      "u max: 2.19309667928\n",
      "u max: 2.17574477406\n",
      "u max: 2.19118494917\n",
      "u max: 2.18136739773\n",
      "u max: 2.25941686882\n",
      "u max: 2.2931865044\n",
      "u max: 2.25285947733\n",
      "u max: 2.16830516194\n",
      "u max: 2.20324673275\n",
      "u max: 2.19442489366\n",
      "u max: 2.21800039369\n",
      "u max: 2.29249656694\n",
      "u max: 2.299671007\n",
      "u max: 2.23086052242\n",
      "u max: 2.18074015283\n",
      "u max: 2.21663527608\n",
      "u max: 2.18093649966\n",
      "u max: 2.2541591662\n",
      "u max: 2.31168518668\n",
      "u max: 2.28791523367\n",
      "u max: 2.19924424663\n",
      "u max: 2.20413866041\n",
      "u max: 2.21616070757\n",
      "u max: 2.19528719142\n",
      "u max: 2.28464430209\n",
      "u max: 2.3164514034\n",
      "u max: 2.26327156863\n",
      "u max: 2.16607462113\n",
      "u max: 2.21851075243\n",
      "u max: 2.20378767477\n",
      "u max: 2.22982341625\n",
      "u max: 2.30649499312\n",
      "u max: 2.30813665535\n",
      "u max: 2.23097726456\n",
      "u max: 2.18946584119\n",
      "u max: 2.22306074058\n",
      "u max: 2.18187486751\n",
      "u max: 2.26251296061\n",
      "u max: 2.31792556373\n",
      "u max: 2.28865450824\n",
      "u max: 2.19574831098\n",
      "u max: 2.20904714363\n",
      "u max: 2.21751343484\n",
      "u max: 2.2007657204\n",
      "u max: 2.2900122321\n",
      "u max: 2.31804784242\n",
      "u max: 2.26048449981\n",
      "u max: 2.16837498971\n",
      "u max: 2.22078045821\n",
      "u max: 2.2022957774\n",
      "u max: 2.23486705256\n",
      "u max: 2.30950382057\n",
      "u max: 2.30693322457\n",
      "u max: 2.22676873628\n",
      "u max: 2.19282179646\n",
      "u max: 2.22337187192\n",
      "u max: 2.17876422034\n",
      "u max: 2.26677036378\n",
      "u max: 2.31897191745\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "\n",
    "out_interval = num_steps / 100;\n",
    "\n",
    "# Time-stepping\n",
    "t = 0\n",
    "for k in range(1,num_steps+1):\n",
    "    # Update current time\n",
    "    t += dt\n",
    "\n",
    "    # Step 1: Tentative velocity step\n",
    "    b1 = assemble(L1)\n",
    "    [bc.apply(b1) for bc in bcu]\n",
    "    solve(A1, u_.vector(), b1, 'bicgstab', 'ilu')\n",
    "\n",
    "    # Step 2: Pressure correction step\n",
    "    b2 = assemble(L2)\n",
    "    [bc.apply(b2) for bc in bcp]\n",
    "    solve(A2, p_.vector(), b2, 'bicgstab', 'ilu')\n",
    "\n",
    "    # Step 3: Velocity correction step\n",
    "    b3 = assemble(L3)\n",
    "    solve(A3, u_.vector(), b3, 'cg', 'sor')\n",
    "\n",
    "    # Plot solution\n",
    "    # plot(u_, title='Velocity')\n",
    "    # plot(p_, title='Pressure')\n",
    "\n",
    "    if k%out_interval ==0 or k==num_steps:\n",
    "        print('u max:', u_.vector().array().max())\n",
    "\n",
    "        # Save solution to files\n",
    "        if 'HDF5File' in globals():\n",
    "            xdmffile_u.write(u_, t)\n",
    "            xdmffile_p.write(p_, t)\n",
    "        else:\n",
    "            vtkfile_u << (u_, t)\n",
    "            vtkfile_p << (p_, t)\n",
    "\n",
    "    # Update previous solution\n",
    "    u_k.assign(u_)\n",
    "    p_k.assign(p_)\n",
    "\n",
    "    # Update progress bar\n",
    "    progress.update(t / T)\n",
    "    \n",
    "interactive()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following figure shows the velocity and pressure at final time visualized in ParaView. \n",
    "It uses the *Glyph* filter to visualize the vector velocity field."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://fenicsproject.org/pub/tutorial/html/mov/navier_stokes_cylinder.ogv\" target=\"_blank\"><img src=\"https://fenicsproject.org/pub/tutorial/html/fig/navier_stokes_cylinder_velocity.png\" alt=\"Movie\" width=\"800\" border=\"10\" /></a>"
   ]
  }
 ],
 "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
}