{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.2 Incompressible Navier-Stokes equations\n",
    "\n",
    "As in [the whetting the appetite demo](../wta/navierstokes.ipynb) we consider the unsteady Navier Stokes equations\n",
    "\n",
    "Find $(u,p):[0,T] \\to (H_{0,D}^1)^d \\times L^2$, s.t.\n",
    "\n",
    "\\begin{align*}\n",
    "\\int_{\\Omega} \\partial_t u \\cdot v + \\int_{\\Omega} \\nu \\nabla u \\nabla v + u \\cdot \\nabla u v - \\int_{\\Omega} \\operatorname{div}(v) p &= \\int f v  && \\forall v \\in (H_{0,D}^1)^d, \\\\ \n",
    "- \\int_{\\Omega} \\operatorname{div}(u) q &= 0 && \\forall q \\in L^2, \\\\\n",
    "\\quad u(t=0) & = u_0\n",
    "\\end{align*}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Schäfer-Turek benchmark\n",
    "We consider the benchmark setup of from http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html . The geometry is a 2D channel with a circular obstacle which is positioned (only slightly) off the center of the channel. The geometry:\n",
    "\n",
    "![title](geometry.png) \n",
    "\n",
    "The viscosity is set to $\\nu = 10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# import libarries, define geometry and generate mesh\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *\n",
    "from netgen.webgui import Draw as DrawGeo\n",
    "shape = Rectangle(2,0.41).Circle(0.2,0.2,0.05).Reverse().Face()\n",
    "shape.edges.name=\"cyl\"\n",
    "shape.edges.Min(X).name=\"inlet\"\n",
    "shape.edges.Max(X).name=\"outlet\"\n",
    "shape.edges.Min(Y).name=\"wall\"\n",
    "shape.edges.Max(Y).name=\"wall\"\n",
    "mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.07)).Curve(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "DrawGeo (shape)\n",
    "Draw (mesh)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Viscosity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# viscosity\n",
    "nu = 0.001"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Taylor-Hood velocity-pressure pair\n",
    "\n",
    "We use a Taylor-Hood discretization of degree $k$. The finite element space for the (vectorial) velocity $\\mathbf{u}$ and the pressure $p$ are:\n",
    "\n",
    "\\begin{align}\n",
    "\\mathbf{u} \\in \\mathbf{V} &= \\{ v \\in H^1(\\Omega) | v|_T \\in \\mathcal{P}^k(T) \\}^2 \\\\\n",
    "p \\in Q &= \\{ q \\in H^1(\\Omega) | q|_T \\in \\mathcal{P}^{k-1}(T) \\}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3\n",
    "V = VectorH1(mesh,order=k, dirichlet=\"wall|cyl|inlet\")\n",
    "Q = H1(mesh,order=k-1)\n",
    "X = V*Q"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The `VectorH1` is the product space of two `H1` spaces with convenience operators for identity, `grad`(`Grad`) and `div`. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Stokes problem for initial values\n",
    "Find $\\mathbf{u} \\in \\mathbf{V}$, $p \\in Q$ so that\n",
    "\n",
    "\\begin{align}\n",
    "\\int_{\\Omega} \\nu \\nabla \\mathbf{u} : \\nabla \\mathbf{v} - \\int_{\\Omega} \\operatorname{div}(\\mathbf{v}) p &= \\int \\mathbf{f}  \\cdot \\mathbf{v}  && \\forall \\mathbf{v} \\in \\mathbf{V}, \\\\- \\int_{\\Omega} \\operatorname{div}(\\mathbf{u}) q &= 0 && \\forall q \\in Q,\n",
    "\\end{align}\n",
    "\n",
    "with boundary conditions:\n",
    "\n",
    " * the inflow boundary data ('inlet') with mean value $1$\n",
    " \n",
    "\\begin{align}\n",
    "u_x(0,y) &= 6 y(0.41-y)/(0.41)^2, \\quad \\int_{0}^{0.41} u_x(0,y) dy = 1 \\\\\n",
    "u_y(0,y) &= 0\n",
    "\\end{align}\n",
    " \n",
    " * \"do-nothing\" boundary conditions on the ouflow ('outlet') and\n",
    " * homogeneous Dirichlet conditions on all other boundaries ('wall')."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "gfu = GridFunction(X)\n",
    "velocity, pressure = gfu.components\n",
    "# parabolic inflow at bc=1:\n",
    "uin = CoefficientFunction((1.5*4*y*(0.41-y)/(0.41*0.41),0))\n",
    "velocity.Set(uin, definedon=mesh.Boundaries(\"inlet\"))\n",
    "\n",
    "Draw(velocity,mesh,\"u\")\n",
    "Draw(pressure,mesh,\"p\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### inflow profile plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib; import numpy as np; import matplotlib.pyplot as plt; \n",
    "%matplotlib inline  \n",
    "s = np.arange(0.0, 0.42, 0.01)\n",
    "bvs = 6*s*(0.41-s)/(0.41)**2 \n",
    "plt.plot(s, bvs)\n",
    "plt.xlabel('y'); plt.ylabel('$u_x(0,y)$'); plt.title('inflow profile')\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To solve the Stokes (and later the Navier-Stokes) problem, we introduce the bilinear form:\n",
    "\n",
    "$$\n",
    "  a((u,p),(v,q)) := \\int_{\\Omega} \\nu \\nabla u : \\nabla v - \\operatorname{div}(v) p - \\operatorname{div}(u) q\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "and directly use it to solve the Stokes problem for suitable initial values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "(u,p), (v,q) = X.TnT()\n",
    "\n",
    "a = BilinearForm(X)\n",
    "stokes = (nu*InnerProduct(grad(u),grad(v))-div(u)*q-div(v)*p)*dx\n",
    "a += stokes\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(X)   \n",
    "f.Assemble()\n",
    "\n",
    "res = f.vec - a.mat*gfu.vec\n",
    "inv_stokes = a.mat.Inverse(X.FreeDofs())\n",
    "gfu.vec.data += inv_stokes * res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "sceneu = Draw(velocity,mesh,\"u\")\n",
    "scenep = Draw(pressure,mesh,\"p\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## IMEX time discretization\n",
    "For the time integration we consider an semi-implicit Euler method (IMEX) where the convection is treated only explicitly and the remaining part implicitly:\n",
    "\n",
    "Find $(\\mathbf{u}^{n+1},p^{n+1}) \\in X = \\mathbf{V} \\times Q$, s.t. for all $(\\mathbf{v},q) \\in X = \\mathbf{V} \\times Q$\n",
    "\n",
    "\\begin{align}\n",
    "\\underbrace{m(\\mathbf{u}^{n+1},\\mathbf{v}) + \\Delta t ~\\cdot~a((\\mathbf{u}^{n+1},p^{n+1}),(\\mathbf{v},q))}_{ \\to M^\\ast} ~=~m(\\mathbf{u}^{n},\\mathbf{v}) - \\Delta t ~\\cdot~c(\\mathbf{u}^{n}; \\mathbf{u}^{n},\\mathbf{v}) \n",
    "\\end{align}\n",
    "\n",
    "with \n",
    "\\begin{align}\n",
    "m(\\mathbf{u},\\mathbf{v}) = \\int \\mathbf{u} \\cdot \\mathbf{v}\n",
    "\\end{align}\n",
    "\n",
    "and \n",
    "\\begin{align}\n",
    "c(\\mathbf{w},\\mathbf{u},\\mathbf{v}) = \\int \\mathbf{w} \\cdot \\nabla \\mathbf{u} \\cdot \\mathbf{v}\n",
    "\\end{align}\n",
    "\n",
    "We prefer the incremental form (as it homogenizes the boundary conditions):\n",
    "\\begin{align}\n",
    "m(\\delta \\mathbf{u}^{n+1},\\mathbf{v}) + \\Delta t ~\\cdot~a((\\delta \\mathbf{u}^{n+1},\\delta p^{n+1}),(\\mathbf{v},q)) ~=~ \\Delta t (-a((\\mathbf{u}^{n},p^n),(\\mathbf{v},q)) -c(\\mathbf{u}^{n}; \\mathbf{u}^{n},\\mathbf{v}))\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dt = 0.001\n",
    "# matrix for implicit part of IMEX(1) scheme:\n",
    "mstar = BilinearForm(X)\n",
    "mstar += InnerProduct(u,v)*dx + dt*stokes\n",
    "mstar.Assemble()\n",
    "inv = mstar.mat.Inverse(X.FreeDofs())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Next, we want to evaluate the convection term $c(u^n;u^n,v)$ for given $u^n$, i.e. a linear form. We can do this by setting up a `LinearForm`\n",
    "```\n",
    "fc = LinearForm(X)\n",
    "fc += (Grad(gfu) * gfu) * v * dx\n",
    "```\n",
    "Then, for $u^n =$`gfu` a call `fc.Assemble()` yields a vector `fc.vec`representing the linear form. \n",
    "\n",
    "Alternatively, we can use the semi-linear form structure of the convection form and put it into a `BilinearForm` (which indeed is able to represent non-bilinear forms as we will discuss in more detail in [unit-3.3](../unit-3.3-nonlinear/nonlinear.ipynb)). We will not need the assembly (or the storage of a sparse matrix) though:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "conv = BilinearForm(X, nonassemble = True)\n",
    "conv += (Grad(u) * u) * v * dx"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Instead of assembling a `LinearForm` in every time step, we use an operator application of a `BilinearForm`, see [unit-2.11](../unit-2.11-matrixfree/matrixfree.ipynb) or the next units for more details.\n",
    "\n",
    "The `BilinearForm`(which indeed is rather a semi-linear form) can compute the linear form at a given state `u`:\n",
    "\n",
    " * either by a call `conv.Apply(gfu.vec,res)` where the result ends up in the vector `res`\n",
    " * alternatively we can use the notation of an operator application `conv.mat * gfu.vec` yielding the same.\n",
    "\n",
    "Below, we will use the latter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "tend = 0\n",
    "gfut = GridFunction(gfu.space,multidim=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# implicit Euler/explicit Euler splitting method:\n",
    "tend += 1\n",
    "gfut.AddMultiDimComponent(gfu.vec)\n",
    "t = 0; cnt = 0\n",
    "while t < tend-0.5*dt:\n",
    "    print (\"\\rt=\", t, end=\"\")\n",
    "\n",
    "    conv.Assemble()\n",
    "    res = (a.mat + conv.mat)* gfu.vec\n",
    "    gfu.vec.data -= dt * inv * res    \n",
    "\n",
    "    t = t + dt; cnt += 1\n",
    "    sceneu.Redraw()\n",
    "    scenep.Redraw()\n",
    "    if cnt % 50 == 0:\n",
    "        gfut.AddMultiDimComponent(gfu.vec)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "Draw (gfut.components[0], mesh, interpolate_multidim=True, animate=True,\n",
    "      min=0, max=1.9, autoscale=False, vectors = True)\n",
    "Draw (gfut.components[1], mesh, interpolate_multidim=True, animate=True,\n",
    "      min=-0.5, max=1, autoscale=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In this benchmark scenario of specific interests are the forces acting on the disk (cylinder). We will discuss how to compute them next."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Supplementary 1: Computing stresses"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "In the previously used benchmark the stresses on the obstacle are evaluated, the so-called drag and lift coefficients. \n",
    "\n",
    "The force acting on the obstacle is\n",
    "\n",
    "$$\n",
    "  F_{\\circ} = (F_D,F_L) = \\int_{\\Gamma_{cyl}} \\sigma_n = \\int_{\\Gamma_{cyl}} (-\\nu \\frac{\\partial u}{\\partial n} + p I) \\cdot n\n",
    "$$\n",
    "\n",
    "The drag/lift coefficients are \n",
    "\n",
    "$$\n",
    "  c_D = \\frac{2 }{\\bar{u} L} F_D, \\quad c_L = \\frac{2 }{\\bar{u} L} F_L.\n",
    "$$\n",
    "\n",
    "where $\\bar{u} = 1$ is the mean inflow velocity and $L = 0.41$ is the channel width."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We use the residual of our discretization to compute the forces. On the boundary degrees of freedoms of the disk we overwrote the equation by prescribing the (homogeneous) boundary conditions. The equations related to these dofs describe the force (im)balance at this boundary. \n",
    "\n",
    "Testing the residual (functional) with the characteristic function on that boundary (in the x- or y-direction) we obtain the integrated stresses (in the x- or y-direction):\n",
    "\n",
    "We define the functions which are characteristic functions (in terms of boundary evaluations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "drag_x_test = GridFunction(X)\n",
    "drag_x_test.components[0].Set(CoefficientFunction((-20.0,0)), definedon=mesh.Boundaries(\"cyl\"))\n",
    "drag_y_test = GridFunction(X)\n",
    "drag_y_test.components[0].Set(CoefficientFunction((0,-20.0)), definedon=mesh.Boundaries(\"cyl\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We will collect drag and lift forces over time and therefore create empty arrays\n",
    "\n",
    "and reset initial data (and time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "time_vals, drag_x_vals, drag_y_vals = [],[],[]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# restoring initial data\n",
    "res = f.vec - a.mat*gfu.vec\n",
    "gfu.vec.data += inv_stokes * res"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "With the same discretization as before."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Remarks:\n",
    " * you can call the following block several times to continue a simulation and collect more data\n",
    " * note that you can reset the array without rerunning the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "t = 0\n",
    "tend = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tend += 1\n",
    "while t < tend-0.5*dt:\n",
    "    print (\"\\rt=\", t, end=\"\")\n",
    "    res = (a.mat + conv.mat) * gfu.vec\n",
    "    gfu.vec.data -= dt * inv * res    \n",
    "    t += dt\n",
    "    time_vals.append( t )\n",
    "    drag_x_vals.append(InnerProduct(res, drag_x_test.vec) )\n",
    "    drag_y_vals.append(InnerProduct(res, drag_y_test.vec) )"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    " Below you can plot the gathered drag and lift values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Plot drag force over time\n",
    "plt.plot(time_vals, drag_x_vals)\n",
    "plt.xlabel('time'); plt.ylabel('drag'); plt.title('drag'); plt.grid(True)\n",
    "plt.show()   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Plot lift force over time\n",
    "plt.plot(time_vals, drag_y_vals)\n",
    "plt.xlabel('time'); plt.ylabel('lift'); plt.title('lift'); plt.grid(True)\n",
    "plt.show()    "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Tasks**\n",
    "\n",
    "If you understood the tutorial so far, you may want to play around with it further. Here are some suggestions on what you could do:\n",
    "\n",
    "* Compute drag and lift forces directly through the surface integrals and compare to the residual based evaluation\n",
    "* Change the geometry and solve a driven cavity example\n",
    "* Use a fully implicit time integration scheme \n",
    "* Use a higher order implict or IMEX scheme for time integration"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}