{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Navier Stokes - blood flow simulation\n",
    "***"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7a25bc65",
   "metadata": {},
   "source": [
    "## Formulation\n",
    "\n",
    "In this example we present a bloodflow simulation in a 3D geometry. The setting is very similar to the previous example. However, we will demonstrate how to prescribe a time dependent expression as a Dirichlet boundary condition. The inlet velocity is prescribed as\n",
    "\\begin{equation*}\n",
    "    \\mathbf{u}_{\\text{in}} = v_{\\text{avg}}(t)\\begin{pmatrix}\n",
    "    2\\,\\frac{R^2-r^2}{R^2}\\\\\n",
    "    0\\\\\n",
    "    0\n",
    "    \\end{pmatrix},\n",
    "\\end{equation*}\n",
    "where $r$ is the distance of the given point from the center point of the inlet. The average velocity at the inlet is changing over time according to the following equation\n",
    "\\begin{equation*}\n",
    "v_\\text{avg}(t) = 0.6\\max\\left(\\frac{t(0.3-t)}{0.15^2}, 0\\right)\n",
    "\\end{equation*}\n",
    "![Mesh](fig/aneurysm.png \"Mesh\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e36dd55",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "Again, we need to import the neccessary packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "12a4d4d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfin as df\n",
    "import numpy as np\n",
    "from time import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dd0280a",
   "metadata": {},
   "source": [
    "Next, we read the mesh from the file and create the boundary markings since ```.xml``` files do not include them. To do that, we create a mesh function and then loop over all the mesh facets and mark them in case they are a part of the desired boundary. We can save the marked mesh and open in ParaView."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0e023291",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = df.Mesh(\"bifurcation.xml\")\n",
    "dim = mesh.topology().dim()\n",
    "bndry = df.MeshFunction(\"size_t\", mesh, dim-1, 0)\n",
    "mesh.init()\n",
    "\n",
    "points = [np.array([-0.9, 0.0, 0.0]), np.array([0.6, 0.0, 0.6]), np.array([0.6, 0.0, -0.6])]\n",
    "normals = [np.array([-1.0, 0.0, 0.0]), np.array([np.sqrt(2), 0.0, np.sqrt(2)]), np.array([np.sqrt(2), 0.0, -np.sqrt(2)])]\n",
    "for f in df.facets(mesh):\n",
    "    bndry[f] = 0\n",
    "    if f.exterior():\n",
    "        bndry[f] = 1\n",
    "        x = f.midpoint().array() #[:dim]\n",
    "        for mark, point, normal in zip([2,3,4], points, normals):\n",
    "            if abs(np.dot(x-point, normal)) <= 0.001:\n",
    "                bndry[f] = mark\n",
    "\n",
    "with df.XDMFFile(\"marked_mesh.xdmf\") as xdmf:\n",
    "    xdmf.write(bndry)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "707092b7",
   "metadata": {},
   "source": [
    "Now we define the function spaces for MINI element as done in the previous example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "25e26cb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "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": "5db01bc1",
   "metadata": {},
   "source": [
    "Let us define some parameters and constants such as density, dynamic viscosity, timestep size and time $T$. The density and dynamic viscosity of blood are approximately $\\rho=1.05\\,\\frac{\\text{g}}{\\text{cm}^3}$ and $\\mu = 0.04\\,\\text{P (Poise)} = 0.004\\,\\text{Pa⋅s}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "30436bb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = df.Constant(0.04) # Poise\n",
    "rho = df.Constant(1.05) # g/cm^3\n",
    "dt = 0.01\n",
    "t_end = 0.35"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4484b3ec",
   "metadata": {},
   "source": [
    "Then we can define the boundary conditions as no slip on the wall and an expression at the inlet. Notice that the average inflow velocity is now prescribed using a function of time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0805d33e",
   "metadata": {},
   "outputs": [],
   "source": [
    "zero_vector = df.Constant((0.0, 0.0, 0.0))\n",
    "bc_walls = df.DirichletBC(W.sub(0), zero_vector, bndry, 1)\n",
    "\n",
    "velocity = lambda t: 0.6*max((t*(0.3-t))/pow(0.15,2), 0.0)\n",
    "inlet_expression = df.Expression((\"v * 2.0 * (pow(R,2) - pow(x[0]-sx,2) - pow(x[1],2) - pow(x[2],2)) / pow(R,2)\", \"0.0\", \"0.0\"), R=0.15, sx=points[0][0], v=velocity(0.0), degree=2)\n",
    "bc_in = df.DirichletBC(W.sub(0), inlet_expression, bndry, 2)\n",
    "\n",
    "bcs = [bc_in, bc_walls]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d7f8f8a",
   "metadata": {},
   "source": [
    "Then we can define the variational problem the same way as in the previous example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ade86e66",
   "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: rho*df.inner(df.grad(u)*u, v)*df.dx + mu*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(1.0) # Backward Euler\n",
    "F = df.Constant(1.0/dt)*rho*df.inner((u-u0), v)*df.dx + theta*F1 + (1-theta)*F0\n",
    "J = df.derivative(F,w)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8152bafd",
   "metadata": {},
   "source": [
    "We set up the nonlinear variational problem and solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a9683cf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "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-8\n",
    "solver.parameters['newton_solver']['relative_tolerance'] = 1e-6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b811c89",
   "metadata": {},
   "source": [
    "Setup XDMF files to save results for viewing in Paraview."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a47eb8e0",
   "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": "7aa1ea63",
   "metadata": {},
   "source": [
    "Again, we go over all the timesteps and solve for velocity and pressure at each of them using the solution from the previous timestep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1ced0da",
   "metadata": {},
   "outputs": [],
   "source": [
    "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={:.2f}s\".format(t))\n",
    "    inlet_expression.v = velocity(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": "803af718",
   "metadata": {},
   "source": [
    "## Complete Code\n",
    "\n",
    "Run the code using the command ```python3 NavierStokes3D.py```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53794f48",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfin as df\n",
    "import numpy as np\n",
    "from time import time\n",
    "\n",
    "mesh = df.Mesh(\"bifurcation.xml\")\n",
    "dim = mesh.topology().dim()\n",
    "bndry = df.MeshFunction(\"size_t\", mesh, dim-1, 0)\n",
    "mesh.init()\n",
    "\n",
    "points = [np.array([-0.9, 0.0, 0.0]), np.array([0.6, 0.0, 0.6]), np.array([0.6, 0.0, -0.6])]\n",
    "normals = [np.array([-1.0, 0.0, 0.0]), np.array([np.sqrt(2), 0.0, np.sqrt(2)]), np.array([np.sqrt(2), 0.0, -np.sqrt(2)])]\n",
    "for f in df.facets(mesh):\n",
    "    bndry[f] = 0\n",
    "    if f.exterior():\n",
    "        bndry[f] = 1\n",
    "        x = f.midpoint().array()\n",
    "        for mark, point, normal in zip([2,3,4], points, normals):\n",
    "            if abs(np.dot(x-point, normal)) <= 0.001:\n",
    "                bndry[f] = mark\n",
    "\n",
    "with df.XDMFFile(\"marked_mesh.xdmf\") as xdmf:\n",
    "    xdmf.write(bndry)\n",
    "\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",
    "mu = df.Constant(0.04) # Poise\n",
    "rho = df.Constant(1.05) # g/cm^3\n",
    "dt = 0.01\n",
    "t_end = 0.3\n",
    "\n",
    "zero_vector = df.Constant((0.0, 0.0, 0.0))\n",
    "bc_walls = df.DirichletBC(W.sub(0), zero_vector, bndry, 1)\n",
    "\n",
    "velocity = lambda t: 0.6*max((t*(0.3-t))/pow(0.15,2), 0.0)\n",
    "inlet_expression = df.Expression((\"v * 2.0 * (pow(R,2) - pow(x[0]-sx,2) - pow(x[1],2) - pow(x[2],2)) / pow(R,2)\", \"0.0\", \"0.0\"), R=0.15, sx=points[0][0], v=velocity(0.0), degree=2)\n",
    "bc_in = df.DirichletBC(W.sub(0), inlet_expression, bndry, 2)\n",
    "\n",
    "bcs = [bc_in, bc_walls]\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: rho*df.inner(df.grad(u)*u, v)*df.dx + mu*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(1.0) # Backward Euler\n",
    "F = df.Constant(1.0/dt)*rho*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-8\n",
    "solver.parameters['newton_solver']['relative_tolerance'] = 1e-6\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={:.2f}s\".format(t))\n",
    "    inlet_expression.v = velocity(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
}