{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "2.1.3 Multigrid and Multilevel Methods\n",
    "===\n",
    "\n",
    "In the previous tutorial $\\S$[2.1.1](unit-2.1.1-preconditioners/preconditioner.ipynb), we saw how to construct and use a geometric multigrid preconditioner using a keyword argument to the `Preconditioner`. This tutorial delves deeper into the construction of such preconditioners.  You will see how additive and multiplicative preconditioners can be implemented efficiently by recursive algorithms using NGSolve's python interface.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "Additive Multilevel Preconditioner\n",
    "---\n",
    "\n",
    "Suppose we have a sequence of hierarchically refined meshes and a sequence of nested finite element spaces\n",
    "\n",
    "$$\n",
    "V_0 \\subset V_1 \\subset \\ldots V_L\n",
    "$$\n",
    "\n",
    "of dimension $N_l = \\operatorname{dim} V_l, l = 0 \\ldots L$.\n",
    "Suppose we also have the prolongation matrices  \n",
    "\n",
    "$$\n",
    "P_l \\in {\\mathbb R}^{N_l \\times N_{l-1}}\n",
    "$$\n",
    "\n",
    "with the property that \n",
    "$\\underline v_l = P_l \\underline v_{l-1}$ where $\\underline v_k$ is the vector of coefficients in the finite element basis expansion of a $v_k \\in V_k$.\n",
    "For the stiffness matrix $A_k$ of a Galerkin discretization, there holds \n",
    "$$\n",
    "A_{l-1} = P_l^T A_l P_l.\n",
    "$$\n",
    "Let $D_l = \\operatorname{diag} A_l$ be the Jacobi preconditioner (or some similar, cheap and local preconditioner).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "The construction of the preconditioner is motivated by the theory of 2-level preconditioners. \n",
    "The 2-level preconditioner involving levels $l-1$ and level $l$, namely \n",
    "$$\n",
    "C_{2L}^{-1} = D_l^{-1} + P_l A_{l-1}^{-1} P_l^T.\n",
    "$$\n",
    "has optimal condition number by the additive Schwarz lemma. \n",
    "However, a direct inversion of the matrix $A_{l-1}$ (is up to constant factor) as expensive as the inversion of $A_l$.\n",
    "\n",
    "The ideal of the multilevel preconditioner is to replace that inverse by a recursion:\n",
    "\\begin{eqnarray*}\n",
    "C_{ML,0}^{-1} & = & A_0^{-1} \\\\\n",
    "C_{ML,l}^{-1} & = & D_l^{-1} + P_l C_{ML,l-1}^{-1} P_l^T \\qquad \\text{for} \\; l = 1, \\ldots L\n",
    "\\end{eqnarray*}\n",
    "\n",
    "We can implement this using a python class with a recursive constructor. We show two implementation techniques, one in long form, and one in short form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.la import EigenValues_Preconditioner\n",
    "\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))\n",
    "\n",
    "fes = H1(mesh,order=1, dirichlet=\".*\", autoupdate=True)\n",
    "u,v = fes.TnT()\n",
    "a = BilinearForm(grad(u)*grad(v)*dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "This is the long form implementation, which is a useful starting point for improvising to more complex preconditioners you might need for other purposes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "class MLPreconditioner(BaseMatrix):\n",
    "    def __init__ (self, fes, level, mat, coarsepre):\n",
    "        super().__init__()\n",
    "        self.fes = fes\n",
    "        self.level = level\n",
    "        self.mat = mat\n",
    "        self.coarsepre = coarsepre\n",
    "        if level > 0:\n",
    "            self.localpre = mat.CreateSmoother(fes.FreeDofs())\n",
    "        else:\n",
    "            self.localpre = mat.Inverse(fes.FreeDofs())\n",
    "        \n",
    "    def Mult (self, x, y):\n",
    "        \"\"\"Return y = C[l] * x = (inv(D) + P * C[l-1] * P.T) * x \n",
    "           when l>0 and inv(A[0]) * x when l=0. \"\"\"\n",
    "        if self.level == 0:\n",
    "            y.data = self.localpre * x\n",
    "            return\n",
    "        hx = x.CreateVector(copy=True)        \n",
    "        self.fes.Prolongation().Restrict(self.level, hx)  # hx <- P.T * x\n",
    "        cdofs = self.fes.Prolongation().LevelDofs(self.level-1)\n",
    "        y[cdofs] = self.coarsepre * hx[cdofs]             # y = C[l-1] * hx  \n",
    "        self.fes.Prolongation().Prolongate(self.level, y) # y <- P * C[l-1] * P.T * x\n",
    "        y += self.localpre * x                            # y += inv(D) * x\n",
    "\n",
    "    def Shape (self):\n",
    "        return self.localpre.shape\n",
    "    def CreateVector (self, col):\n",
    "        return self.localpre.CreateVector(col)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "Here is the considerably shorter version implementing the same preconditioner using operator notation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def MLPreconditioner2(fes, level, mat, coarsepre):\n",
    "    prol = fes.Prolongation().Operator(level)     # get P \n",
    "    localpre = mat.CreateSmoother(fes.FreeDofs()) # get inv(D)\n",
    "    return localpre + prol @ coarsepre @ prol.T   # Return: inv(D) + P * C[l-1] * P.T * x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "a.Assemble()   # Make exact inverse at current/coarsest level\n",
    "pre = a.mat.Inverse(fes.FreeDofs())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "for l in range(9):\n",
    "    mesh.Refine()\n",
    "    a.Assemble()\n",
    "    pre = MLPreconditioner(fes,l+1, a.mat, pre)    \n",
    "    lam = EigenValues_Preconditioner(a.mat, pre)\n",
    "    print(\"ndof=%7d:  minew=%.4f  maxew=%1.4f  Cond# = %5.3f\" \n",
    "          %(fes.ndof, lam[0], lam[-1], lam[-1]/lam[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "Multiplicative Multigrid Preconditioning\n",
    "---\n",
    "\n",
    "The multigrid preconditioner combines the local preconditioner and the coarse grid step multiplicatively or sequentially. The multigrid preconditioning action is defined as follows:\n",
    "\n",
    "Multigrid $C_{MG,l}^{-1} : d \\mapsto w$\n",
    "\n",
    "* if l = 0, set $w = A_0^{-1} d$ and return\n",
    "* w = 0\n",
    "* presmoothing, $m_l$ steps:   \n",
    "$$ w \\leftarrow w + D_{pre}^{-1} (d - A w)$$ \n",
    "* coasre grid correction:\n",
    "$$\n",
    " w \\leftarrow w + P_l C_{MG,l-1}^{-1} P_l^{T} (d - A w)\n",
    "$$\n",
    "* postsmoothing, $m_l$ steps:   \n",
    "$$w \\leftarrow w + D_{post}^{-1} (d - A w)$$ \n",
    "\n",
    "\n",
    "If the preconditioners $D_{pre}$ and $D_{post}$ from the pre-smoothing and post-smoothing iterations are transposed to each other, the overall preconditioner is symmetric. If the pre- and post-smoothing iterations are non-expansive, the overall preconditioner is positive definite.\n",
    "These properties can be realized using forward Gauss-Seidel for pre-smoothing, and backward Gauss-Seidel for post-smoothing. Here is its implementation as a recursive python class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.la import EigenValues_Preconditioner\n",
    "\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))\n",
    "\n",
    "fes = H1(mesh,order=1, dirichlet=\".*\", autoupdate=True)\n",
    "u,v = fes.TnT()\n",
    "a = BilinearForm(grad(u)*grad(v)*dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "class MGPreconditioner(BaseMatrix):\n",
    "    def __init__ (self, fes, level, mat, coarsepre):\n",
    "        super().__init__()\n",
    "        self.fes = fes\n",
    "        self.level = level\n",
    "        self.mat = mat\n",
    "        self.coarsepre = coarsepre\n",
    "        if level > 0:\n",
    "            self.localpre = mat.CreateSmoother(fes.FreeDofs())\n",
    "        else:\n",
    "            self.localpre = mat.Inverse(fes.FreeDofs())\n",
    "        \n",
    "    def Mult (self, d, w):\n",
    "        if self.level == 0:\n",
    "            w.data = self.localpre * d\n",
    "            return\n",
    "        \n",
    "        prol = self.fes.Prolongation().Operator(self.level)\n",
    "\n",
    "        w[:] = 0\n",
    "        self.localpre.Smooth(w,d)\n",
    "        res  = d - self.mat * w\n",
    "        w += prol @ self.coarsepre @ prol.T * res\n",
    "        self.localpre.SmoothBack(w,d)\n",
    "\n",
    "\n",
    "    def Shape (self):\n",
    "        return self.localpre.shape\n",
    "    def CreateVector (self, col):\n",
    "        return self.localpre.CreateVector(col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "a.Assemble()\n",
    "pre = MGPreconditioner(fes, 0, a.mat, None)\n",
    "    \n",
    "for l in range(9):\n",
    "    mesh.Refine()\n",
    "    a.Assemble()\n",
    "    pre = MGPreconditioner(fes,l+1, a.mat, pre)    \n",
    "    lam = EigenValues_Preconditioner(a.mat, pre)\n",
    "    print(\"ndof=%7d:  minew=%.4f  maxew=%1.4f  Cond# = %5.3f\" \n",
    "          %(fes.ndof, lam[0], lam[-1], lam[-1]/lam[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "A quick conjugate gradient call on this finest grid gives you an idea of typical number of iterations and practial efficiency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(1*v*dx).Assemble()\n",
    "gfu = GridFunction(fes)\n",
    "from ngsolve.krylovspace import CGSolver\n",
    "inv = CGSolver(mat=a.mat, pre=pre, printrates=True)\n",
    "gfu.vec.data = inv * f.vec"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16",
   "metadata": {},
   "source": [
    "Projection matrices from the finest level\n",
    "---\n",
    "It is often not feasible to assemble matrices on the coarse level, \n",
    "for example when solving non-linear problems. Then it is useful to \n",
    "calculate coarse grid matrices from the matrix on the finest level using the Galerkin property\n",
    "\n",
    "$$\n",
    "A_{l-1} = P_{l}^T A_l P_l.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.geom2d import unit_square\n",
    "from ngsolve import *\n",
    "from ngsolve.la import EigenValues_Preconditioner\n",
    "\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))\n",
    "\n",
    "fes = H1(mesh,order=1, dirichlet=\".*\", autoupdate=True)\n",
    "u,v = fes.TnT()\n",
    "a = BilinearForm(grad(u)*grad(v)*dx)\n",
    "\n",
    "for l in range(8):\n",
    "     mesh.Refine()\n",
    "    \n",
    "a.Assemble();  # only matrix at the finest level provided to mg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "class ProjectedMG(BaseMatrix):\n",
    "    def __init__ (self, fes, mat, level):\n",
    "        super(ProjectedMG, self).__init__()\n",
    "        self.fes = fes\n",
    "        self.level = level\n",
    "        self.mat = mat\n",
    "        if level > 0:\n",
    "            self.prol = fes.Prolongation().CreateMatrix(level)\n",
    "            self.rest = self.prol.CreateTranspose()\n",
    "            coarsemat = self.rest @ mat @ self.prol # multiply matrices\n",
    "            self.localpre = mat.CreateSmoother(fes.FreeDofs())\n",
    "                \n",
    "            self.coarsepre = ProjectedMG(fes, coarsemat, level-1)\n",
    "        else:\n",
    "            self.localpre = mat.Inverse(fes.FreeDofs())\n",
    "        \n",
    "    def Mult (self, d, w):\n",
    "        if self.level == 0:\n",
    "            w.data = self.localpre * d\n",
    "            return\n",
    "        w[:] = 0\n",
    "        self.localpre.Smooth(w,d)\n",
    "        res = d - self.mat * w\n",
    "        w += self.prol @ self.coarsepre @ self.rest * res\n",
    "        self.localpre.SmoothBack(w,d)\n",
    "        \n",
    "    def Shape (self):\n",
    "        return self.localpre.shape\n",
    "    def CreateVector (self, col):\n",
    "        return self.localpre.CreateVector(col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "pre = ProjectedMG(fes, a.mat, fes.mesh.levels-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "lam = EigenValues_Preconditioner(a.mat, pre)\n",
    "print(\"ndof=%7d:  minew=%.4f  maxew=%1.4f  Cond# = %5.3f\" \n",
    "          %(fes.ndof, lam[0], lam[-1], lam[-1]/lam[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(1*v*dx).Assemble()\n",
    "gfu = GridFunction(fes)\n",
    "from ngsolve.krylovspace import CGSolver\n",
    "inv = CGSolver(mat=a.mat, pre=pre, printrates=True)\n",
    "gfu.vec.data = inv * f.vec"
   ]
  }
 ],
 "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": 5
}