SteadyNavierStokes_Newton-StreamStab.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Lab9_Ex3.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "04NkQfW4nElA"
      },
      "source": [
        "# FE approximation of the Navier-Stokes problem with Newton linearization"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jkGbxL1FoHAq"
      },
      "source": [
        "We consider again the steady Navier-Stokes problem defined in the domain $\\Omega := {[0,1]}^2$ by the following system of PDEs and Dirichlet BCs: \n",
        "$$\n",
        "\\begin{cases}\n",
        "  ({\\boldsymbol{u}}\\cdot\\boldsymbol{\\nabla})\\ \\boldsymbol{u}\n",
        "  -\\dfrac{1}{\\text{Re}}\\boldsymbol{\\Delta} \\boldsymbol{u} + \\boldsymbol{\\nabla} p = 0, \n",
        "  \\qquad&\\text{in }\\Omega, \\\\\n",
        "  \\boldsymbol{\\nabla} \\cdot \\boldsymbol{u} = 0, \n",
        "  \\qquad &\\text{in }\\Omega, \\\\\n",
        "  \\boldsymbol{u} = 1\\boldsymbol{i} + 0\\boldsymbol{j}, \\qquad&\\text{on }\\Gamma^{\\text{up}} = \\{ 0 \\le x \\le 1, \\, y=1 \\}, \\\\\n",
        "  \\boldsymbol{u} = \\boldsymbol{0}, \\qquad&\\text{on }\\partial\\Omega \\backslash \\Gamma^{\\text{up}}.\n",
        "\\end{cases}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8kz6v8yVo9Ww"
      },
      "source": [
        "The linearized problem for the Netwon's method is given by\n",
        "$$\n",
        "  a(\\boldsymbol{u}_k,\\boldsymbol{v}) \n",
        "  + c(\\boldsymbol{u}_k,\\boldsymbol{u}_{k-1},\\boldsymbol{v})\n",
        "  + c(\\boldsymbol{u}_{k-1},\\boldsymbol{u}_k,\\boldsymbol{v}) \n",
        "  - b(\\boldsymbol{v},p_k) + b(\\boldsymbol{u}_k,q)\n",
        "  =\n",
        "c(\\boldsymbol{u}_{k-1},\\boldsymbol{u}_{k-1},\\boldsymbol{v}),\n",
        "$$\n",
        "where\n",
        "$$\n",
        "  a(\\boldsymbol{u},\\boldsymbol{v}):= \\frac{1}{\\text{Re}}\n",
        "  \\int_\\Omega\\boldsymbol{\\nabla}{\\boldsymbol{u}}:\\boldsymbol{\\nabla}{\\boldsymbol{v}}\\, d\\boldsymbol{x},\\qquad\n",
        "  b(\\boldsymbol{u},p) := -\\int_\\Omega \\boldsymbol{\\nabla}\\cdot{\\boldsymbol{u}}\\ p \\,d\\boldsymbol{x},\\qquad\n",
        "c(\\boldsymbol{u},\\boldsymbol{v},\\boldsymbol{w}) := \\int_\\Omega (\\boldsymbol{u}\\cdot\\boldsymbol{\\nabla})\\boldsymbol{v}\\ \\boldsymbol{w} \\, d\\vec{x}.\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nglPBz_HoyNw"
      },
      "source": [
        "## Numerical approximation with FEniCS"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ILHUaw3Jo1Eh"
      },
      "source": [
        "%%capture\n",
        "try:\n",
        "    import dolfin\n",
        "except ImportError:\n",
        "    !wget \"https://fem-on-colab.github.io/releases/fenics-install.sh\" -O \"/tmp/fenics-install.sh\" && bash \"/tmp/fenics-install.sh\"\n",
        "    import dolfin\n",
        "\n",
        "from fenics import *"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_9T8IIgqofsr"
      },
      "source": [
        "Using the stable pair of spaces $\\mathbb{P}^2/\\mathbb{P}^1$ solve the problem using Netwon's method. Consider a stopping criterion based on the relative increment\n",
        "$$\n",
        "\\frac{||\\boldsymbol{u}_k-\\boldsymbol{u}_{k-1}||_{H^1}}{||\\boldsymbol{u}_{k-1}||_{H^1}}\n",
        "+\\frac{||\\boldsymbol{p}_k-\\boldsymbol{p}_{k-1}||_{L^2}}{||\\boldsymbol{p}_{k-1}||_{L^2}}\n",
        "\\le 10^{-6}.\n",
        "$$\n",
        "\n",
        "*   Try to solve the problem with different values of Reynolds number ($\\text{Re}=10,100,1000$). \n",
        "*   Assess the convergence performance of the iterative method."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KiXSp65XnA20",
        "outputId": "265e276c-a5c4-45fb-da50-7d699b60c0a6",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 114
        }
      },
      "source": [
        "# 1. mesh generation\n",
        "n = 15\n",
        "mesh = UnitSquareMesh(n, n, 'crossed')\n",
        "\n",
        "# 2. finite element space\n",
        "V = VectorElement('CG', mesh.ufl_cell(), 2)\n",
        "Q = FiniteElement('CG', mesh.ufl_cell(), 1)\n",
        "X = FunctionSpace(mesh, V*Q)\n",
        "\n",
        "u_boundary = Expression((\n",
        "        'near(x[1], 1) ? 1.0 : 0.0',\n",
        "        '0'\n",
        "    ), degree=0)\n",
        "\n",
        "def boundary(x, on_boundary):\n",
        "  return on_boundary\n",
        "\n",
        "def origin(x, on_boundary):\n",
        "  return near(x[0], 0) and near(x[1], 0)\n",
        "\n",
        "bc = [DirichletBC(X.sub(0), u_boundary, boundary),\n",
        "      DirichletBC(X.sub(1), Constant(0), origin, 'pointwise')]\n",
        "\n",
        "# 3. problem definition\n",
        "def solve_linear(X, advection, Re, f, bc):\n",
        "  u, p = TrialFunctions(X)\n",
        "  v, q = TestFunctions(X)\n",
        "\n",
        "  a = (inner(grad(u), grad(v)) / Re - p * div(v) + div(u) * q) * dx\n",
        "  L = dot(f, v) * dx\n",
        "  if advection:\n",
        "    # only the first term was present in the fixed point iterative method\n",
        "    a += (dot(grad(u) * advection, v) + dot(grad(advection) * u, v)) * dx\n",
        "    L += dot(grad(advection) * advection, v) * dx\n",
        "  \n",
        "  x = Function(X)\n",
        "  solve(a == L, x, bc)\n",
        "  \n",
        "  return x.split()\n",
        "\n",
        "# 4. solution\n",
        "Re = Constant(100)\n",
        "f = Constant((0, 0))\n",
        "\n",
        "# compute initial guess\n",
        "uh, ph = solve_linear(X, None, Re, f, bc)\n",
        "\n",
        "niter = 50\n",
        "tolerance = 1e-6\n",
        "for i in range(niter):\n",
        "  uh_old, ph_old = uh, ph\n",
        "  uh, ph = solve_linear(X, uh_old, Re, f, bc)\n",
        "  \n",
        "  error = (errornorm(uh, uh_old, 'H1') / norm(uh_old, 'H1') +\n",
        "           errornorm(ph, ph_old, 'L2') / norm(ph_old, 'L2'))\n",
        "  print('step {}: {:.3e}'.format(i, error))\n",
        "  if error < tolerance:\n",
        "    break\n",
        "\n",
        "uh.rename('velocity', 'velocity')\n",
        "ph.rename('pressure', 'pressure')\n",
        "\n",
        "File('velocity.pvd') << uh\n",
        "File('pressure.pvd') << ph"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "step 0: 9.746e-01\n",
            "step 1: 1.199e-01\n",
            "step 2: 1.329e-03\n",
            "step 3: 2.000e-07\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JsqH6RE8sunp"
      },
      "source": [
        "### Streamline diffusion stabilizition of the $\\mathbb{P}^1-\\mathbb{P}^1$ FEM"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Bu7u1MrPtJcV"
      },
      "source": [
        "If the pair of spaces $\\mathbb{P}^1/\\mathbb{P}^1$ is considered, then the linearized problem with the stabilized terms is given by\n",
        "\\begin{multline*}\n",
        "  a(\\boldsymbol{u}_k,\\boldsymbol{v}) \n",
        "  + c(\\boldsymbol{u}_k,\\boldsymbol{u}_{k-1},\\boldsymbol{v})\n",
        "  + c(\\boldsymbol{u}_{k-1},\\boldsymbol{u}_k,\\boldsymbol{v}) \n",
        "  + b(\\boldsymbol{v},p_k) - b(\\boldsymbol{u}_k,q) + \\\\\n",
        "  \\sum_{K\\in\\mathcal{T}_h} \\delta_K \\left[{\\left(\n",
        "    (\\boldsymbol{u}_k\\cdot\\boldsymbol{\\nabla})\n",
        "    \\boldsymbol{u}_{k-1}\n",
        "    +(\\boldsymbol{u}_{k-1}\\cdot\\boldsymbol{\\nabla})\n",
        "    \\boldsymbol{u}_{k}\n",
        "    + \\boldsymbol{\\nabla}{p},\n",
        "    (\\boldsymbol{v}\\cdot\\boldsymbol{\\nabla})\n",
        "    \\boldsymbol{u}_{k-1}\n",
        "    +(\\boldsymbol{u}_{k-1}\\cdot\\boldsymbol{\\nabla})\n",
        "    \\boldsymbol{v}\n",
        "    + \\boldsymbol{\\nabla}{q}\\right)}_K\n",
        "    + (\\boldsymbol{\\nabla}\\cdot{\\boldsymbol{u}_k}, \n",
        "    \\boldsymbol{\\nabla}\\cdot{\\boldsymbol{v}})_K\\right]= \\\\ \n",
        "c(\\boldsymbol{u}_{k-1},\\boldsymbol{u}_{k-1},\\boldsymbol{v})\n",
        "+ \\sum_{K\\in\\mathcal{T}_h} \\delta_K {\\left(\n",
        "  (\\boldsymbol{u}_{k-1}\\cdot\\boldsymbol{\\nabla})\n",
        "  \\boldsymbol{u}_{k-1},\n",
        "  (\\boldsymbol{v}\\cdot\\boldsymbol{\\nabla})\n",
        "  \\boldsymbol{u}_{k-1}\n",
        "  +(\\boldsymbol{u}_{k-1}\\cdot\\boldsymbol{\\nabla})\n",
        "  \\boldsymbol{v}\n",
        "  + \\boldsymbol{\\nabla}{q}\\right)}_K,\n",
        "\\end{multline*}\n",
        "where ${(\\cdot, \\cdot)}_K$ denotes the $L^2$ inner product over the element $K$. The stabilization parameters are defined by\n",
        "$$\n",
        "\\delta_K = \n",
        "\\delta \\min{\\left\\{ \\frac{h_K}{2 ||\\boldsymbol{u}||}_K, \\frac{h^2_K \\text{Re}}{24} \\right\\}},\\qquad \\delta > 0.\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qlUIUTG7rrsN",
        "outputId": "c14f5716-79f9-4751-f48a-2f7d2beed457",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 146
        }
      },
      "source": [
        "# 1. mesh generation\n",
        "n = 15\n",
        "mesh = UnitSquareMesh(n, n, 'crossed')\n",
        "\n",
        "# 2. finite element space\n",
        "V = VectorElement('CG', mesh.ufl_cell(), 1)\n",
        "Q = FiniteElement('CG', mesh.ufl_cell(), 1)\n",
        "X = FunctionSpace(mesh, V*Q)\n",
        "\n",
        "u_boundary = Expression((\n",
        "        'near(x[1], 1) ? 1.0 : 0.0',\n",
        "        '0'\n",
        "    ), degree=0)\n",
        "\n",
        "def boundary(x, on_boundary):\n",
        "  return on_boundary\n",
        "\n",
        "def origin(x, on_boundary):\n",
        "  return near(x[0], 0) and near(x[1], 0)\n",
        "\n",
        "bc = [DirichletBC(X.sub(0), u_boundary, boundary),\n",
        "      DirichletBC(X.sub(1), Constant(0), origin, 'pointwise')]\n",
        "\n",
        "# 3. problem definition\n",
        "def solve_linear(X, advection, Re, f, bc):\n",
        "  u, p = TrialFunctions(X)\n",
        "  v, q = TestFunctions(X)\n",
        "\n",
        "  a = (inner(grad(u), grad(v)) / Re - p * div(v) + div(u) * q) * dx\n",
        "  L = dot(f, v) * dx\n",
        "  if advection:\n",
        "    a += (dot(grad(u) * advection, v) + dot(grad(advection) * u, v)) * dx\n",
        "    L += dot(grad(advection) * advection, v) * dx\n",
        "    \n",
        "    A = lambda u, p: (grad(u) * advection) + (grad(advection) * u) + grad(p)\n",
        "    A_SS = lambda u, p: (grad(u) * advection) + (grad(advection) * u) + grad(p)\n",
        "  \n",
        "    h = CellDiameter(X.mesh())\n",
        "    anorm = sqrt(dot(advection, advection))\n",
        "    tau_K = 0.5 * h / conditional(anorm * h * Re > 12, anorm, 12 / h / Re)\n",
        "  \n",
        "    a += tau_K * (dot(A(u, p), A_SS(v, q)) + div(u) * div(v)) * dx\n",
        "    L += tau_K * (dot(f + grad(advection) * advection, A_SS(v, q))) * dx\n",
        "  \n",
        "  x = Function(X)\n",
        "  solve(a == L, x, bc)\n",
        "  \n",
        "  return x.split()\n",
        "\n",
        "# 4. solution\n",
        "Re = Constant(100)\n",
        "f = Constant((0, 0))\n",
        "\n",
        "uh, ph = solve_linear(X, None, Re, f, bc)\n",
        "\n",
        "niter = 50\n",
        "tolerance = 1e-6\n",
        "for i in range(niter):\n",
        "  uh_old, ph_old = uh, ph\n",
        "  uh, ph = solve_linear(X, uh_old, Re, f, bc)\n",
        "  \n",
        "  error = (errornorm(uh, uh_old, 'H1') / norm(uh_old, 'H1') +\n",
        "           errornorm(ph, ph_old, 'L2') / norm(ph_old, 'L2'))\n",
        "  print('step {}: {:.3e}'.format(i, error))\n",
        "  if error < tolerance:\n",
        "    break\n",
        "\n",
        "uh.rename('velocity_stab', 'velocity')\n",
        "ph.rename('pressure_stab', 'pressure')\n",
        "\n",
        "File('velocity_stab.pvd') << uh\n",
        "File('pressure_stab.pvd') << ph"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "step 0: 1.266e+00\n",
            "step 1: 1.047e-01\n",
            "step 2: 9.925e-04\n",
            "step 3: 9.449e-06\n",
            "step 4: 1.828e-07\n"
          ],
          "name": "stdout"
        }
      ]
    }
  ]
}