{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stokes flow around a sphere\n",
    "***"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Formulation\n",
    "\n",
    "In this example, we numerically reproduce the famous Stokes law (see Landau, Lifschitz: Fluid mechanics, Pergamon Press, Oxford, 1966)\n",
    "\n",
    "$$ F_\\mathrm{drag} = 6\\pi \\mu R v, $$\n",
    "\n",
    "which evaluates the drag on a sphere of radius $R$ moving at speed $v$. We load an external mesh which was prepared using [Gmsh](https://gmsh.info/) software.\n",
    "\n",
    "<div align=center>\n",
    "<img src=\"fig/geometry.png\" width=\"350\"/>\n",
    "</div>\n",
    "\n",
    "Instead of moving the sphere in a time-dependant geometry, it is much simpler to fix the sphere and let the fluid move around it. Also, since FEM works on bounded domains, we confine the flow in a finite cylinder of radius $R_\\mathrm{cylinder} = 100R$ and length $L_\\mathrm{cylinder} = 100R$. At the walls and the inflow, we define a Dirichlet boundary condition\n",
    "\n",
    "$$ \\mathbf{v} = \\begin{pmatrix}\n",
    "                 1 \\\\\n",
    "                 0 \\\\\n",
    "                 0 \\\\\n",
    "                \\end{pmatrix},\n",
    "$$\n",
    "\n",
    "homogenous Dirichlet condition on sphere and natural boundary condition at the outflow. For $R = 1$ and $\\mu = 1$, the drag predicted by theory is $6\\pi$.\n",
    "\n",
    "This mesh has coarse resolution of the cylinder but is significantly refined near the sphere. Ability to work with non-uniform resolution is a significant advantage of FEM.\n",
    "\n",
    "## Implementation\n",
    "\n",
    "We begin by importing few things. Since we will run this on multiple cores using MPI, we need to access the rank of threads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfin as df\n",
    "from time import time\n",
    "from numpy import pi\n",
    "\n",
    "comm = df.MPI.comm_world\n",
    "rank = df.MPI.rank(comm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the mesh from a file called *mesh_big.h5*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = df.Mesh()\n",
    "hdf = df.HDF5File(mesh.mpi_comm(), \"mesh_big.h5\", \"r\")\n",
    "hdf.read(mesh, \"/mesh\", False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read boundary marker from the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "bdary = df.MeshFunction(\"size_t\", mesh, mesh.topology().dim()-1, 0)\n",
    "hdf.read(bdary, \"/boundaries\")\n",
    "#1 = INFLOW\n",
    "#2 = WALLS\n",
    "#3 = OUTFLOW\n",
    "#4 = SPHERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define function spaces. For extra speed, let us use Crouzeix-Raviart elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ev = df.VectorElement(\"CR\", mesh.ufl_cell(), 1)\n",
    "Ep = df.FiniteElement(\"DG\", mesh.ufl_cell(), 0)\n",
    "W = df.FunctionSpace(mesh, df.MixedElement([Ev, Ep]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we define boundary conditions as stated above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "V0 = df.Constant((0.0, 0.0, 0.0))\n",
    "Vin = df.Constant((1.0, 0.0, 0.0))\n",
    "bc_inflow = df.DirichletBC(W.sub(0), Vin, bdary, 1)\n",
    "bc_walls = df.DirichletBC(W.sub(0), Vin, bdary, 2)\n",
    "bc_sphere = df.DirichletBC(W.sub(0), V0, bdary, 4)\n",
    "bc = [bc_inflow, bc_walls, bc_sphere]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variational form and solver will be same as in two dimensions. When printing output, we use ```if rank == 0``` condition. (Otherwise, it would be printed by every thread.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "nu = df.Constant(1.0)\n",
    "u, p = df.TrialFunctions(W)\n",
    "v, q = df.TestFunctions(W)\n",
    "F = nu*df.inner(df.grad(u), df.grad(v))*df.dx - p*df.div(v)*df.dx - q*df.div(u)*df.dx\n",
    "\n",
    "w = df.Function(W)\n",
    "problem = df.LinearVariationalProblem(df.lhs(F), df.rhs(F), w, bc)\n",
    "solver = df.LinearVariationalSolver(problem)\n",
    "solver.parameters[\"linear_solver\"] = \"mumps\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve the linear system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ellapsed = Solving linear variational problem.\n",
      " 111.6617522239685 s\n"
     ]
    }
   ],
   "source": [
    "tick = time()\n",
    "solver.solve()\n",
    "\n",
    "if rank == 0:\n",
    "    print(\"ellapsed = \", time() - tick, \"s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the result for inspection in [Paraview](https://www.paraview.org/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "ufile = df.XDMFFile(\"results/u.xdmf\")\n",
    "pfile = df.XDMFFile(\"results/p.xdmf\")\n",
    "u, p = w.split()\n",
    "u.rename(\"u\", \"velocity\")\n",
    "p.rename(\"p\", \"pressure\")\n",
    "ufile.write(u)\n",
    "pfile.write(p)\n",
    "ufile.close()\n",
    "pfile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Report the drag. Compare it to the Stokes' law."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "drag by surface integral = 18.92\n",
      "drag by volume integral = 18.92\n",
      "drag according to Stokes' law = 18.85\n",
      "relative error = 0.3954%\n"
     ]
    }
   ],
   "source": [
    "n = df.FacetNormal(mesh)\n",
    "I = df.Identity(mesh.geometry().dim())\n",
    "ds = df.Measure(\"ds\", subdomain_data=bdary)\n",
    "T = -p*I + 2.0*nu*df.sym(df.grad(u))\n",
    "force = df.dot(T, n)\n",
    "drag = df.assemble(-force[0]*ds(4))\n",
    "\n",
    "#alt force computation\n",
    "bc_ = df.DirichletBC(W.sub(0), (1.0, 0.0, 0.0), bdary, 4)\n",
    "w_ = df.Function(W)\n",
    "bc_.apply(w_.vector())\n",
    "drag_ = -df.assemble(df.inner(T, df.grad(w_.sub(0)))*df.dx)\n",
    "\n",
    "if rank == 0:\n",
    "        print(\"drag by surface integral = {:.4}\".format(drag))\n",
    "        print(\"drag by volume integral = {:.4}\".format(drag_))\n",
    "        print(\"drag according to Stokes' law = {:.4}\".format(6*pi))\n",
    "        print(\"relative error = {:.4%}\".format(drag/(6*pi) - 1.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remark\n",
    "\n",
    "In *mesh_big.h5*, the ratio of cylinder radius to sphere radius is 100:1. This is not an overkill. You can try to recompute this problem with *mesh_small.h5*, where the ratio is 10:1. You will find that the boundary effect is strong and leads to much higher drag."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complete Code \n",
    "\n",
    "Run in parallel using the command\n",
    "\n",
    "```$ mpirun -n 4 python3 Stokes3D.py ```\n",
    "\n",
    "where 4 is number of threads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfin as df\n",
    "from time import time\n",
    "from numpy import pi\n",
    "\n",
    "comm = df.MPI.comm_world\n",
    "rank = df.MPI.rank(comm)\n",
    "\n",
    "mesh = df.Mesh()\n",
    "hdf = df.HDF5File(mesh.mpi_comm(), \"mesh_big.h5\", \"r\")\n",
    "hdf.read(mesh, \"/mesh\", False)\n",
    "\n",
    "bdary = df.MeshFunction(\"size_t\", mesh, mesh.topology().dim()-1, 0)\n",
    "hdf.read(bdary, \"/boundaries\")\n",
    "#1 = INFLOW\n",
    "#2 = WALLS\n",
    "#3 = OUTFLOW\n",
    "#4 = SPHERE\n",
    "\n",
    "Ev = df.VectorElement(\"CR\", mesh.ufl_cell(), 1)\n",
    "Ep = df.FiniteElement(\"DG\", mesh.ufl_cell(), 0)\n",
    "W = df.FunctionSpace(mesh, df.MixedElement([Ev, Ep]))\n",
    "\n",
    "V0 = df.Constant((0.0, 0.0, 0.0))\n",
    "Vin = df.Constant((1.0, 0.0, 0.0))\n",
    "bc_inflow = df.DirichletBC(W.sub(0), Vin, bdary, 1)\n",
    "bc_walls = df.DirichletBC(W.sub(0), Vin, bdary, 2)\n",
    "bc_sphere = df.DirichletBC(W.sub(0), V0, bdary, 4)\n",
    "bc = [bc_inflow, bc_walls, bc_sphere]\n",
    "\n",
    "nu = df.Constant(1.0)\n",
    "u, p = df.TrialFunctions(W)\n",
    "v, q = df.TestFunctions(W)\n",
    "F = nu*df.inner(df.grad(u), df.grad(v))*df.dx - p*df.div(v)*df.dx - q*df.div(u)*df.dx\n",
    "\n",
    "w = df.Function(W)\n",
    "problem = df.LinearVariationalProblem(df.lhs(F), df.rhs(F), w, bc)\n",
    "solver = df.LinearVariationalSolver(problem)\n",
    "solver.parameters[\"linear_solver\"] = \"mumps\"\n",
    "\n",
    "tick = time()\n",
    "solver.solve()\n",
    "\n",
    "if rank == 0:\n",
    "    print(\"ellapsed = \", time() - tick, \"s\")\n",
    "\n",
    "ufile = df.XDMFFile(\"results/u.xdmf\")\n",
    "pfile = df.XDMFFile(\"results/p.xdmf\")\n",
    "u, p = w.split()\n",
    "u.rename(\"u\", \"velocity\")\n",
    "p.rename(\"p\", \"pressure\")\n",
    "ufile.write(u)\n",
    "pfile.write(p)\n",
    "\n",
    "n = df.FacetNormal(mesh)\n",
    "I = df.Identity(mesh.geometry().dim())\n",
    "ds = df.Measure(\"ds\", subdomain_data=bdary)\n",
    "T = -p*I + 2.0*nu*df.sym(df.grad(u))\n",
    "force = df.dot(T, n)\n",
    "drag = df.assemble(-force[0]*ds(4))\n",
    "\n",
    "if rank == 0:\n",
    "        print(\"drag according to simulation = {:.4}\".format(drag))\n",
    "        print(\"drag according to Stokes' law = {:.4}\".format(6*pi))\n",
    "        print(\"relative error = {:.4%}\".format(drag/(6*pi) - 1.0))"
   ]
  }
 ],
 "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.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}