{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Navier Stokes Equations\n",
    "===\n",
    "Find velocity $u : \\Omega \\times [0,T] \\rightarrow R^d$ and pressure $p : \\Omega \\times [0,T] \\rightarrow R$ such that\n",
    "\n",
    "$$ \n",
    "\\begin{array}{ccccl}\n",
    "\\frac{\\partial u}{\\partial t} - \\nu \\Delta u + u \\nabla u & + & \\nabla p & = & f \\\\\n",
    "\\operatorname{div} u & & & = & 0\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "import ipywidgets as widgets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Schäfer-Turek benchmark geometry:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "\n",
    "shape = Rectangle(2,0.41).Circle(0.2,0.2,0.05).Reverse().Face()\n",
    "shape.edges.name=\"wall\"\n",
    "shape.edges.Min(X).name=\"inlet\"\n",
    "shape.edges.Max(X).name=\"outlet\"\n",
    "Draw (shape);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.07)).Curve(3)\n",
    "Draw (mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Higher order Taylor-Hood element pairing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = VectorH1(mesh,order=3, dirichlet=\"wall|cyl|inlet\")\n",
    "Q = H1(mesh,order=2)\n",
    "X = V*Q\n",
    "\n",
    "u,p = X.TrialFunction()\n",
    "v,q = X.TestFunction()\n",
    "\n",
    "nu = 0.001  # viscosity\n",
    "stokes = (nu*InnerProduct(grad(u), grad(v))+ \\\n",
    "    div(u)*q+div(v)*p - 1e-10*p*q)*dx\n",
    "\n",
    "a = BilinearForm(stokes).Assemble()\n",
    "\n",
    "# nothing here ...\n",
    "f = LinearForm(X).Assemble()\n",
    "\n",
    "# gridfunction for the solution\n",
    "gfu = GridFunction(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "parabolic inflow at inlet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uin = CoefficientFunction( (1.5*4*y*(0.41-y)/(0.41*0.41), 0) )\n",
    "gfu.components[0].Set(uin, definedon=mesh.Boundaries(\"inlet\"))\n",
    "Draw (gfu.components[0], mesh, \"vel\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "solve Stokes problem for initial conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inv_stokes = a.mat.Inverse(X.FreeDofs())\n",
    "\n",
    "res = f.vec - a.mat*gfu.vec\n",
    "gfu.vec.data += inv_stokes * res\n",
    "\n",
    "Draw (gfu.components[0], mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "implicit/explicit time-stepping:\n",
    "\n",
    "$$\n",
    "\\frac{u_{n+1}-u_n}{\\tau} - \\nu \\Delta u_{n+1} + \\nabla p_{n+1} = f - u_n \\nabla u_n\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\operatorname{div} u_{n+1} = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = 0.001 # timestep\n",
    "\n",
    "mstar = BilinearForm(u*v*dx+tau*stokes).Assemble()\n",
    "inv = mstar.mat.Inverse(X.FreeDofs(), inverse=\"sparsecholesky\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the non-linear convective term $\\int u \\nabla u v$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conv = BilinearForm(X, nonassemble=True)\n",
    "conv += (Grad(u) * u) * v * dx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "implicit Euler/explicit Euler splitting method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = 0; i = 0\n",
    "tend = 10\n",
    "gfut = GridFunction(V, multidim=0)\n",
    "vel = gfu.components[0]\n",
    "\n",
    "scene = Draw (gfu.components[0], mesh, min=0, max=2, autoscale=False)\n",
    "tw = widgets.Text(value='t = 0')\n",
    "display(tw)\n",
    "\n",
    "with TaskManager():\n",
    "    while t < tend:\n",
    "        res = conv.Apply(gfu.vec) + a.mat*gfu.vec\n",
    "        gfu.vec.data -= tau * inv * res    \n",
    "\n",
    "        t = t + tau; i = i + 1\n",
    "        if i%10 == 0: scene.Redraw()\n",
    "        if i%50 == 0: gfut.AddMultiDimComponent(vel.vec)\n",
    "        if i%100 == 0: tw.value = \"t = \" + str(t)\n",
    "    tw.value = \"t = \"+str(tend)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the multidim - GridFunction gfut we have collected several time-steps, which can be animated now by the visualization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfut, mesh, interpolate_multidim=True, animate=True, min=0, max=2, autoscale=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}