{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.1.1 Preconditioners in NGSolve\n",
    "\n",
    "Preconditioners are approximate inverses which are used within iterative methods to solve linear or non-linear equations.\n",
    "\n",
    "Here are some built-in preconditioners in NGSolve:\n",
    "\n",
    "* Jacobi (`local`) and block Jacobi \n",
    "* Direct solvers using sparse factorization (`direct`)\n",
    "* Geometric multigrid with different block-smoothers (`multigrid`)\n",
    "* Algebraic multigrid preconditioner (`h1amg`)\n",
    "* p-version element-level BDDC (`bddc`)\n",
    "\n",
    "This tutorial quickly shows how to use some of these within a solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A simple test problem \n",
    "\n",
    "In order to experiment with various preconditioners, \n",
    "let us define a simple Poisson solver with the name\n",
    "of a preconditioner as the keyword argument `precond`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveProblem(h=0.5, p=1, levels=1, \n",
    "                 condense=False,\n",
    "                 precond=\"local\"):\n",
    "    \"\"\"\n",
    "    Solve Poisson problem on l refinement levels.\n",
    "    PARAMETERS:\n",
    "        h: coarse mesh size\n",
    "        p: polynomial degree \n",
    "        l: number of refinement levels\n",
    "        condense: if true, perform static condensations\n",
    "        precond: name of a built-in preconditioner\n",
    "    Returns: (ndof, niterations)\n",
    "        List of tuples of number of degrees of freedom and iterations\n",
    "    \"\"\"\n",
    "    mesh = Mesh(unit_square.GenerateMesh(maxh=h))\n",
    "    # mesh = Mesh(unit_cube.GenerateMesh(maxh=h))\n",
    "    fes = H1(mesh, order=p, dirichlet=\"bottom|left\")\n",
    "    \n",
    "    u, v = fes.TnT() \n",
    "    a = BilinearForm(grad(u)*grad(v)*dx, condense=condense)\n",
    "    f = LinearForm(v*dx)\n",
    "    gfu = GridFunction(fes)\n",
    "    Draw(gfu)\n",
    "    c = Preconditioner(a, precond) # Register c to a BEFORE assembly\n",
    "\n",
    "    steps = []\n",
    "    \n",
    "    for l in range(levels):\n",
    "        if l > 0: \n",
    "            mesh.Refine()\n",
    "        fes.Update()\n",
    "        gfu.Update()\n",
    "\n",
    "        with TaskManager():\n",
    "            a.Assemble()\n",
    "            f.Assemble()\n",
    "\n",
    "            # Conjugate gradient solver\n",
    "            inv = CGSolver(a.mat, c.mat, maxsteps=1000)\n",
    "\n",
    "            # Solve steps depend on condense \n",
    "            if condense:\n",
    "                f.vec.data += a.harmonic_extension_trans * f.vec\n",
    "            \n",
    "            gfu.vec.data = inv * f.vec\n",
    "            \n",
    "            if condense:\n",
    "                gfu.vec.data += a.harmonic_extension * gfu.vec\n",
    "                gfu.vec.data += a.inner_solve * f.vec\n",
    "        steps.append ( (fes.ndof, inv.GetSteps()) )\n",
    "        if fes.ndof < 15000:\n",
    "            Redraw()\n",
    "    return steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Preconditioner` registers itself to the `BilinearForm`. Whenever the `BilinearForm` is re-assembled, the `Preconditioner` is updated as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `local` preconditioner \n",
    "\n",
    "The `local` preconditioner is a simple Jacobi preconditioner. \n",
    "The number of CG-iterations with the local preconditioner is proportional to $h^{-1} \\sim 2^l$:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SolveProblem(precond=\"local\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_local = SolveProblem(levels=9, precond=\"local\")\n",
    "res_local"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multigrid preconditioner \n",
    "\n",
    "A geometric multigrid `Preconditioner` uses the sequence of refined meshes to define a preconditioner yielding optimal iteration numbers (and complexity). It uses a direct solve on the coarsest level, and block Gauss-Seidel smoothers on the refined levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_mg = SolveProblem(levels=9, precond=\"multigrid\")\n",
    "res_mg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.plot(*zip(*res_local), \"-*\")\n",
    "plt.plot(*zip(*res_mg), \"-+\")\n",
    "plt.legend(['local', 'mg'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Multigrid implementation for higher order spaces \n",
    "\n",
    "For high order finite element spaces NGSolve uses hierarchical bases.\n",
    "Let the (small) sub-spaces $V_E$, $V_F$, and $V_C$ be those generated by basis functions associated with an edge $E$, a face $F$, and a cell $C$, respectively. Then the full\n",
    "space admits the decomposition\n",
    "\n",
    "$$\n",
    "V_{hp} = V_{p=1} + \\sum_{\\text{edges }E} V_E + \\sum_{\\text{faces }F} V_F + \\sum_{\\text{cells }C} V_C\n",
    "$$\n",
    "\n",
    "where $V_{p=1}$ refers to the lowest order finite element subspace.\n",
    "The system matrix then takes the block form\n",
    "\n",
    "$$\n",
    "A = \\left( \\begin{array}{cccc}\n",
    "A_{VV} & A_{VE} & A_{VF} & A_{VC} \\\\\n",
    "A_{EV} & A_{EE} & A_{EF} & A_{EC} \\\\\n",
    "A_{FV} & A_{FE} & A_{FF} & A_{FC} \\\\\n",
    "A_{CV} & A_{CE} & A_{CF} & A_{CC} \\\\\n",
    "\\end{array} \\right)\n",
    "$$\n",
    "\n",
    "where the $A_{VV}$-block is exactly the system matrix of a lowest order method, representing the part of operator acting on the $V_{p=1}$ subspace.\n",
    "\n",
    "NGSolve's *multigrid implementation for a high order method uses h-version multigrid for the lowest order block,* and  local block-smoothing for the high order bubble basis functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for p in range(1,10):\n",
    "    r = SolveProblem(h=0.5, p=p, levels=4, condense=False, \n",
    "                     precond=\"multigrid\")\n",
    "    print (\"p=\", p, \": ndof,nsteps=\", r)          "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that the number of iterations grows mildly with the degree $p$ while remaining  bounded with mesh refinement.\n",
    "\n",
    "Performing static condensation improves the situation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for p in range(1,10):\n",
    "    r = SolveProblem(h=0.5, p=p, levels=4, condense=True, \n",
    "                     precond=\"multigrid\")\n",
    "    print (\"p=\", p, \": ndof,nsteps=\", r)       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Element-wise BDDC preconditioner\n",
    "\n",
    "A built-in element-wise BDDC (Balancing Domain Decomposition preconditioner with Constraints) preconditioner is also available. In contrast to local or multigrid preconditioners, the BDDC preconditioner needs access to the element matrices. This is exactly why we need to register the preconditioner with the bilinear form `bfa` before calling `bfa.Assemble()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for p in range(1,10):\n",
    "    r = SolveProblem(h=0.5, p=p, levels=4, condense=True, \n",
    "                     precond=\"bddc\")\n",
    "    print (\"p=\", p, \": ndof,nsteps=\", r)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The BDDC preconditioner needs more iterations, but the work per iteration is less, so performance is similar to multigrid. **This element-wise BDDC preconditioner works well in shared memory parallel as well as in distributed memory mode.** See $\\S$[2.1.4](../unit-2.1.4-bddc/bddc.ipynb) for more about BDDC preconditioner and how to combine it with an algebraic multigrid coarse solver. "
   ]
  }
 ],
 "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}