{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Navier Stokes - flow around cylinder\n",
    "***"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "11c94af0",
   "metadata": {},
   "source": [
    "## Formulation\n",
    "\n",
    "In this [benchmark example](https://wwwold.mathematik.tu-dortmund.de/~featflow/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark3_re100.html) we will show how to compute the flow of a Newtonian fluid around a cylinder. Even though the inlet velocity profile is constant, the resulting velocity field may not be stationary. We will see a pattern called [Kármán vortex street](https://en.wikipedia.org/wiki/K%C3%A1rm%C3%A1n_vortex_street).\n",
    "\n",
    "The following figure shows the computational mesh used. \n",
    "\n",
    "<div align=center>\n",
    "<img src=\"fig/hron_turek_bench.png\" width=\"500\"/>\n",
    "</div>\n",
    "\n",
    "The strong formulation of the time dependent Navier Stokes equations for this problem is as follows\n",
    "\\begin{align}\n",
    "\\frac{\\partial\\mathbf{u}}{\\partial t} + \\mathbf{u}\\cdot\\nabla\\mathbf{u} + \\nabla p - \\nu\\Delta\\mathbf{u} &= 0\\quad\\text{in }\\,\\Omega\\times (0,T)\\\\\n",
    "\\nabla\\cdot\\mathbf{u}&=0\\quad\\text{in }\\,\\Omega\\times (0,T)\\\\\n",
    "\\mathbf{u}(t=0)&=\\mathbf{0}\\quad\\text{in }\\,\\Omega\\\\\n",
    "\\mathbf{u} &= \\mathbf{u}_{\\text{in}}\\quad\\text{on }\\,\\Gamma_{\\text{in}}\\times (0,T)\\\\\n",
    "\\mathbf{u} &= \\mathbf{0}\\quad\\text{on }\\,\\Gamma_{\\text{wall}}\\times (0,T)\\\\\n",
    "\\mathbf{u} &= \\mathbf{0}\\quad\\text{on }\\,\\Gamma_{\\text{cylinder}}\\times (0,T)\\\\\n",
    "(-p \\mathbb{I} + \\nu \\nabla \\mathbf{u}) \\mathbf{n} &= \\mathbf{0}\\quad\\text{on }\\,\\Gamma_{\\text{out}}\\times (0,T),\n",
    "\\end{align}\n",
    "where $\\mathbf{u}_{\\text{in}}$ is defined as \n",
    "$$\n",
    "\\mathbf{u}_{\\text{in}} = v_{\\text{avg}}\\begin{pmatrix}\n",
    "\\frac{4\\,y\\,(D-y)}{D^2}\\\\\n",
    "0\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Taking test functions $q$ and $\\mathbf{v}$ for pressure and velocity respectively, we multiply the first equation by $\\mathbf{v}$ and the second equation by $q$ and integrate over $\\Omega$. Then we use the integration by parts to obtain the weak formulation.\n",
    "\n",
    "$$\n",
    "\\int_\\Omega\\frac{\\partial\\mathbf{u}}{\\partial t}\\cdot\\mathbf{v}\\,\\text{d}x+\\int_\\Omega(\\mathbf{u}\\cdot\\nabla\\mathbf{u})\\cdot\\mathbf{v}\\,\\text{d}x+\\nu\\int_\\Omega\\nabla\\mathbf{u}\\cdot\\nabla\\mathbf{v}\\,\\text{d}x-\\int_\\Omega p\\,\\text{div}\\,\\mathbf{v}\\,\\text{d}x-\\int_\\Omega q\\,\\text{div}\\,\\mathbf{u}\\,\\text{d}x = 0\n",
    "$$\n",
    "\n",
    "Compared to the Stokes equations, the Navier Stokes equations include a convection term $\\int_\\Omega(\\mathbf{u}\\cdot\\nabla\\mathbf{u})\\cdot\\mathbf{v}\\,\\text{d}x$. Therefore, a nonlinear solver has to be used instead of a linear solver. \n",
    "\n",
    "Another new term is the time derivative $\\frac{\\partial \\mathbf{u}}{\\partial t}$ which we will treat using a Crank-Nicholson time stepping scheme.\n",
    "First we discretize the time derivative.\n",
    "\\begin{equation}\n",
    "\\frac{\\partial \\mathbf{u}}{\\partial t}(t) \\approx \\frac{\\mathbf{u}(t)-\\mathbf{u}(t-\\delta t)}{\\delta t}\n",
    "\\end{equation}\n",
    "In order to use the Crank-Nicolson scheme, we make a linear combination of explicit and implicit scheme using a factor $\\theta$ in part of the equation (here we use simpler notation $\\mathbf{u}^{i} = \\mathbf{u}(t)$, $\\mathbf{u}^{i-1} = \\mathbf{u}(t-\\delta t)$). Notice that the last two terms are not part of the linear combination since we need to use the implicit scheme on them to have stability.\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\int_\\Omega\\frac{\\mathbf{u}^i-\\mathbf{u}^{i-1}}{\\delta t}\\cdot\\mathbf{v}\\,\\text{d}x\\\\\n",
    "+\\theta\\left(\\int_\\Omega(\\mathbf{u}^i\\cdot\\nabla\\mathbf{u}^i)\\cdot\\mathbf{v}\\,\\text{d}x+\\nu\\int_\\Omega\\nabla\\mathbf{u}^i\\cdot\\nabla\\mathbf{v}\\,\\text{d}x\\right)\\\\+(1-\\theta)\\left(\\int_\\Omega(\\mathbf{u}^{i-1}\\cdot\\nabla\\mathbf{u}^{i-1})\\cdot\\mathbf{v}\\,\\text{d}x\\\\+\\nu\\int_\\Omega\\nabla\\mathbf{u}^{i-1}\\cdot\\nabla\\mathbf{v}\\,\\text{d}x\\right)\\\\\n",
    "-\\int_\\Omega p^i\\,\\text{div}\\,\\mathbf{v}\\,\\text{d}x-\\int_\\Omega q\\,\\text{div}\\,\\mathbf{u}^i\\,\\text{d}x = 0\n",
    "\\end{split}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29785573",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "First, we need to import all neccessary packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a514d189",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfin as df\n",
    "from time import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c7d7b9",
   "metadata": {},
   "source": [
    "Then we read mesh and boundary markers from the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a1a3e872",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = df.Mesh()\n",
    "hdf = df.HDF5File(mesh.mpi_comm(), \"bench_csg.h5\", \"r\")\n",
    "hdf.read(mesh, \"/mesh\", False)\n",
    "\n",
    "bndry = df.MeshFunction(\"size_t\", mesh, mesh.topology().dim()-1, 0)\n",
    "hdf.read(bndry, \"/boundaries\")\n",
    "\n",
    "with df.XDMFFile(\"results/boundaries.xdmf\") as xdmf:\n",
    "    xdmf.write(bndry)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb420be3",
   "metadata": {},
   "source": [
    "Next, we define the function spaces V, P and create the mixed function space W. The [MINI](https://defelement.com/elements/mini.html) element is used in this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b07279b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = mesh.geometry().dim()\n",
    "U = df.FiniteElement(\"CG\", mesh.ufl_cell(), 1)\n",
    "B = df.FiniteElement(\"Bubble\", mesh.ufl_cell(), dim + 1)\n",
    "P = df.FiniteElement(\"CG\", mesh.ufl_cell(), 1)\n",
    "V = df.VectorElement(df.NodalEnrichedElement(U, B))\n",
    "W = df.FunctionSpace(mesh, df.MixedElement([V, P]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3622069",
   "metadata": {},
   "source": [
    "Let us define some parameters and constants such are the average velocity at the inlet, kinematic viscosity, timestep size and time T."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "84be18b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity = 1.5\n",
    "nu = df.Constant(0.001)\n",
    "dt = 0.1\n",
    "t_end = 15"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50b86bf7",
   "metadata": {},
   "source": [
    "Let us define the boundary conditions. The boundary markings are inlet: 1, outlet: 2, wall: 3, cylinder: 5. We prescribe zero vector on the wall and the cylinder using Dirichlet boundary conditions. Inlet velocity can be defined using an expression object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "81448128",
   "metadata": {},
   "outputs": [],
   "source": [
    "inlet_expression = df.Expression((\"v * 4.0 * x[1] * (D - x[1]) / pow(D,2)\", \"0.0\"), D=0.41, v=velocity, degree=2)\n",
    "bc_in = df.DirichletBC(W.sub(0), inlet_expression, bndry, 1)\n",
    "\n",
    "zero_vector = df.Constant((0.0, 0.0))\n",
    "bc_walls = df.DirichletBC(W.sub(0), zero_vector, bndry, 3)\n",
    "bc_sphere = df.DirichletBC(W.sub(0), zero_vector, bndry, 5)\n",
    "\n",
    "bcs = [bc_in, bc_walls, bc_sphere]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed209f00",
   "metadata": {},
   "source": [
    "Now we define the variational form for functions $\\mathbf{u}, p$ with help of test functions $\\mathbf{v}, q$ and functions $\\mathbf{u}_0, p_0$ representing the velocity and pressure fields from previous timestep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f4a9a974",
   "metadata": {},
   "outputs": [],
   "source": [
    "v, q = df.TestFunctions(W)\n",
    "w = df.Function(W)\n",
    "w0 = df.Function(W)\n",
    "u, p = df.split(w)\n",
    "u0, p0 = df.split(w0)\n",
    "\n",
    "a = lambda u, v: df.inner(df.grad(u)*u, v)*df.dx + nu*df.inner(df.grad(u), df.grad(v))*df.dx\n",
    "b = lambda q, v: q*df.div(v)*df.dx\n",
    "\n",
    "F1 = a(u, v) - b(p, v) - b(q, u)\n",
    "F0 = a(u0, v) - b(p, v) - b(q, u)\n",
    "\n",
    "theta = df.Constant(0.5)\n",
    "F = df.Constant(1.0/dt)*df.inner((u-u0), v)*df.dx + theta*F1 + (1-theta)*F0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16148c04",
   "metadata": {},
   "source": [
    "We create a nonlinear variational problem object using our variational form and a nonlinear variational solver. Optionally, we can set up some of the solver parameters such as absolute and relative tolerance for convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0264d286",
   "metadata": {
    "scrolled": true
   },
   "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"
     ]
    }
   ],
   "source": [
    "problem = df.NonlinearVariationalProblem(F, w, bcs, J=df.derivative(F, w))\n",
    "solver = df.NonlinearVariationalSolver(problem)\n",
    "solver.parameters['newton_solver']['linear_solver'] = 'mumps'\n",
    "solver.parameters['newton_solver']['absolute_tolerance'] = 1e-12\n",
    "solver.parameters['newton_solver']['relative_tolerance'] = 1e-12"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41500b6b",
   "metadata": {},
   "source": [
    "Setup XDMF files to save results for viewing in Paraview."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "da238b27",
   "metadata": {},
   "outputs": [],
   "source": [
    "ufile = df.XDMFFile(\"results/u.xdmf\")\n",
    "pfile = df.XDMFFile(\"results/p.xdmf\")\n",
    "ufile.parameters[\"flush_output\"] = True\n",
    "pfile.parameters[\"flush_output\"] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b653819f",
   "metadata": {},
   "source": [
    "Since we have time dependent problem, we need to solve the problem for each timestep using a loop. We start at time $t=0$, then for each timestep, we assign the previous solution $w$ to $w_0$ and solve for the new $w$. It is also a good idea to save the solution at each timestep into the XDMF file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9e3687b",
   "metadata": {},
   "outputs": [],
   "source": [
    "tick = time()\n",
    "t = 0.0\n",
    "(u, p) = w.split(deepcopy=True)\n",
    "u.rename(\"v\", \"velocity\")\n",
    "p.rename(\"p\", \"pressure\")\n",
    "df.assign(u, w.sub(0))\n",
    "df.assign(p, w.sub(1))\n",
    "ufile.write(u, t)\n",
    "pfile.write(p, t)\n",
    "while t < t_end:\n",
    "    w0.assign(w)\n",
    "    t += dt\n",
    "    print(\"t={:.1f}s\".format(t))\n",
    "    solver.solve()\n",
    "    df.assign(u, w.sub(0))\n",
    "    df.assign(p, w.sub(1))\n",
    "    ufile.write(u, t)\n",
    "    pfile.write(p, t)\n",
    "print(\"ellapsed = \", time() - tick, \"s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a92137bc",
   "metadata": {},
   "source": [
    "## Complete Code\n",
    "\n",
    "Run the code using the command ```python3 NavierStokes2D.py```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "530f6620",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfin as df\n",
    "from time import time\n",
    "\n",
    "mesh = df.Mesh()\n",
    "hdf = df.HDF5File(mesh.mpi_comm(), \"bench_csg.h5\", \"r\")\n",
    "hdf.read(mesh, \"/mesh\", False)\n",
    "\n",
    "bndry = df.MeshFunction(\"size_t\", mesh, mesh.topology().dim()-1, 0)\n",
    "hdf.read(bndry, \"/boundaries\")\n",
    "\n",
    "with df.XDMFFile(\"results/boundaries.xdmf\") as xdmf:\n",
    "    xdmf.write(bndry)\n",
    "    \n",
    "dim = mesh.geometry().dim()\n",
    "U = df.FiniteElement(\"CG\", mesh.ufl_cell(), 1)\n",
    "B = df.FiniteElement(\"Bubble\", mesh.ufl_cell(), dim + 1)\n",
    "P = df.FiniteElement(\"CG\", mesh.ufl_cell(), 1)\n",
    "V = df.VectorElement(df.NodalEnrichedElement(U, B))\n",
    "W = df.FunctionSpace(mesh, df.MixedElement([V, P]))\n",
    "\n",
    "velocity = 1.5\n",
    "nu = df.Constant(0.001)\n",
    "dt = 0.1\n",
    "t_end = 15\n",
    "\n",
    "inlet_expression = df.Expression((\"v * 4.0 * x[1] * (D - x[1]) / pow(D,2)\", \"0.0\"), D=0.41, v=velocity, degree=2)\n",
    "bc_in = df.DirichletBC(W.sub(0), inlet_expression, bndry, 1)\n",
    "\n",
    "zero_vector = df.Constant((0.0, 0.0))\n",
    "bc_walls = df.DirichletBC(W.sub(0), zero_vector, bndry, 3)\n",
    "bc_sphere = df.DirichletBC(W.sub(0), zero_vector, bndry, 5)\n",
    "\n",
    "bcs = [bc_in, bc_walls, bc_sphere]\n",
    "\n",
    "v, q = df.TestFunctions(W)\n",
    "w = df.Function(W)\n",
    "w0 = df.Function(W)\n",
    "u, p = df.split(w)\n",
    "u0, p0 = df.split(w0)\n",
    "\n",
    "a = lambda u, v: df.inner(df.grad(u)*u, v)*df.dx + nu*df.inner(df.grad(u), df.grad(v))*df.dx\n",
    "b = lambda q, v: q*df.div(v)*df.dx\n",
    "\n",
    "F1 = a(u, v) - b(p, v) - b(q, u)\n",
    "F0 = a(u0, v) - b(p, v) - b(q, u)\n",
    "\n",
    "theta = df.Constant(0.5)\n",
    "F = df.Constant(1.0/dt)*df.inner((u-u0), v)*df.dx + theta*F1 + (1-theta)*F0\n",
    "J = df.derivative(F,w)\n",
    "\n",
    "problem = df.NonlinearVariationalProblem(F, w, bcs, J)\n",
    "solver = df.NonlinearVariationalSolver(problem)\n",
    "solver.parameters['newton_solver']['linear_solver'] = 'mumps'\n",
    "solver.parameters['newton_solver']['absolute_tolerance'] = 1e-12\n",
    "solver.parameters['newton_solver']['relative_tolerance'] = 1e-12\n",
    "\n",
    "ufile = df.XDMFFile(\"results/u.xdmf\")\n",
    "pfile = df.XDMFFile(\"results/p.xdmf\")\n",
    "ufile.parameters[\"flush_output\"] = True\n",
    "pfile.parameters[\"flush_output\"] = True\n",
    "\n",
    "tick = time()\n",
    "t = 0.0\n",
    "(u, p) = w.split(True)\n",
    "u.rename(\"v\", \"velocity\")\n",
    "p.rename(\"p\", \"pressure\")\n",
    "df.assign(u, w.sub(0))\n",
    "df.assign(p, w.sub(1))\n",
    "ufile.write(u, t)\n",
    "pfile.write(p, t)\n",
    "while t < t_end:\n",
    "    w0.assign(w)\n",
    "    t += dt\n",
    "    print(\"t={:.1f}s\".format(t))\n",
    "    solver.solve()\n",
    "    df.assign(u, w.sub(0))\n",
    "    df.assign(p, w.sub(1))\n",
    "    ufile.write(u, t)\n",
    "    pfile.write(p, t)\n",
    "print(\"ellapsed = \", time() - tick, \"s\")"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}