{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.1.4 Element-wise BDDC Preconditioner\n",
    "\n",
    "The element-wise BDDC (Balancing Domain Decomposition preconditioner with Constraints) preconditioner in NGSolve is a good general purpose preconditioner that works well both in the shared memory parallel mode as well as in distributed memory mode. In this tutorial, we discuss this preconditioner, related built-in options, and customization from python.   \n",
    "\n",
    "Let us start with a simple description of the element-wise BDDC in the context of **Lagrange** finite element space $V$. Here the BDDC preconditioner is constructed on an auxiliary space $\\widetilde{V}$ obtained by connecting only element vertices (leaving edge and face shape functions discontinuous). Although larger, the auxiliary space allows local elimination of edge and face variables. Hence an analogue of the  original matrix $A$ on this space, named $\\widetilde A$, is less expensive to invert.  This inverse is used to construct a preconditioner for $A$ as follows:\n",
    "\n",
    "$$\n",
    "C_{BDDC}^{-1} = R {\\,\\widetilde{A}\\,}^{-1}\\, R^t\n",
    "$$\n",
    "\n",
    "Here, $R$ is the averaging operator for the discontinuous edge and face variables. \n",
    "\n",
    "To construct a general purpose BDDC preconditioner, NGSolve **generalizes** this idea to all its finite element spaces by a classification of degrees of freedom. Dofs are classified into (condensable) `LOCAL_DOF`s that we saw in [1.4](../unit-1.4-staticcond/staticcond.ipynb) and a remainder that includes these types: \n",
    "\n",
    "`WIREBASKET_DOF`  \n",
    "`INTERFACE_DOF`\n",
    "\n",
    "The original finite element space $V$ is obtained by requiring conformity of both the above types of dofs, while the auxiliary space $\\widetilde{V}$ is obtained by requiring conformity of `WIREBASKET_DOF`s only."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a figure of a typical function in the default $\\widetilde{V}$ (and the code to generate this is at the end of this tutorial) when $V$ is the Lagrange space:\n",
    "\n",
    "![title](figs/auxspace.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from ngsolve.la import EigenValues_Preconditioner\n",
    "SetHeapSize(100*1000*1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = Mesh(unit_cube.GenerateMesh(maxh=0.5))\n",
    "# mesh = Mesh(unit_square.GenerateMesh(maxh=0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Built-in options\n",
    "\n",
    "Let us define a simple function to study how the spectrum of the preconditioned matrix changes with various options.\n",
    "\n",
    "#### Effect of condensation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TestPreconditioner (p, condense=False, **args):\n",
    "    fes = H1(mesh, order=p, **args)\n",
    "    u,v = fes.TnT()\n",
    "    a = BilinearForm(fes, eliminate_internal=condense)\n",
    "    a += grad(u)*grad(v)*dx + u*v*dx\n",
    "    c = Preconditioner(a, \"bddc\")\n",
    "    a.Assemble()   \n",
    "    return EigenValues_Preconditioner(a.mat, c.mat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lams = TestPreconditioner(5)\n",
    "print (lams[0:3], \"...\\n\", lams[-3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the effect of static condensation on the BDDC preconditioner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lams = TestPreconditioner(5, condense=True)\n",
    "print (lams[0:3], \"...\\n\", lams[-3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tuning the auxiliary space\n",
    "\n",
    "Next, let us study the effect of a few built-in flags for *finite element spaces* that are useful for tweaking the behavior of the BDDC preconditioner. The effect of these flags varies in two (2D) and three dimensions (3D), e.g., \n",
    "\n",
    "- `wb_fulledges=True`: This option keeps **all** edge-dofs conforming (i.e. they are marked `WIREBASKET_DOF`s). This option is only meaningful in 3D. If used in 2D, the preconditioner becomes a direct solver.\n",
    "\n",
    "- `wb_withedges=True`: This option keeps only the **first** edge-dof conforming (i.e., the first edge-dof is marked `WIREBASKET_DOF` and the remaining edge-dofs are marked `INTERFACE_DOF`s). \n",
    "\n",
    "The complete situation is a bit more complex due to the fact these options  can take the three values *True, False,* or *Undefined*, the two options can be combined, and the space dimension can be 2 or 3. The *default* value of these flags in NGSolve is  *Undefined*. Here is a table with the summary of the effect of these options:\n",
    "\n",
    "| wb_fulledges      |  wb_withedges |  2D    |  3D    |\n",
    "|-------------------|---------------|--------|--------|\n",
    "|   True            |   any value   |  all   |   all  |\n",
    "|   False/Undefined |   Undefined   |  none  |  first |\n",
    "|   False/Undefined |    False      |  none  |  none  |\n",
    "|   False/Undefined |    True       |  first |  first |\n",
    "\n",
    "An entry $X \\in$ {all, none, first} of the last two columns is to be read as follows: $X$ of the edge-dofs is(are) `WIREBASKET_DOF`(s)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example of the effect of one of these flag values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lams = TestPreconditioner(5, condense=True,\n",
    "                          wb_withedges=False)\n",
    "print (lams[0:3], \"...\\n\", lams[-3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, when moving from the default case (where the first of the edge dofs are wire basket dofs) to the case (where none of the edge dofs are wire basket dofs), the conditioning became less favorable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Customize \n",
    "\n",
    "From within python, we can change the types of degrees of freedom of finite element spaces, thus affecting the behavior of the BDDC preconditioner.\n",
    "\n",
    "To depart from the default and commit the **first two** edge dofs to wire basket, we perform the next steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=10)\n",
    "u,v = fes.TnT()\n",
    "\n",
    "for ed in mesh.edges:\n",
    "    dofs = fes.GetDofNrs(ed)\n",
    "    for d in dofs:\n",
    "        fes.SetCouplingType(d, COUPLING_TYPE.INTERFACE_DOF)\n",
    "\n",
    "    # Set the first two edge dofs to be conforming\n",
    "    fes.SetCouplingType(dofs[0], COUPLING_TYPE.WIREBASKET_DOF)\n",
    "    fes.SetCouplingType(dofs[1], COUPLING_TYPE.WIREBASKET_DOF)\n",
    "\n",
    "a = BilinearForm(fes, eliminate_internal=True)\n",
    "a += grad(u)*grad(v)*dx + u*v*dx\n",
    "c = Preconditioner(a, \"bddc\")\n",
    "a.Assemble()\n",
    "\n",
    "lams=EigenValues_Preconditioner(a.mat, c.mat)\n",
    "max(lams)/min(lams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a slight improvement from the default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lams = TestPreconditioner (10, condense=True)\n",
    "max(lams)/min(lams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combine BDDC with AMG for large problems "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coarse inverse ${\\,\\widetilde{A}\\,}^{-1}$ of BDDC is expensive on fine meshes. Using the option `coarsetype=h1amg` flag, we can ask BDDC to replace \n",
    "${\\,\\widetilde{A}\\,}^{-1}$ by an Algebraic MultiGrid (AMG) preconditioner. Since NGSolve's `h1amg` is well-suited  \n",
    "for the lowest order Lagrange space,  we use the option \n",
    "`wb_withedges=False` to ensure that $\\widetilde{A}$ is made solely with vertex unknowns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 5\n",
    "mesh = Mesh(unit_cube.GenerateMesh(maxh=0.05))\n",
    "fes = H1(mesh, order=p, dirichlet=\"left|bottom|back\",\n",
    "         wb_withedges=False)\n",
    "print(\"NDOF = \", fes.ndof)\n",
    "u,v = fes.TnT()\n",
    "a = BilinearForm(fes)\n",
    "a += grad(u)*grad(v)*dx\n",
    "f = LinearForm(fes)\n",
    "f += v*dx\n",
    "\n",
    "with TaskManager():\n",
    "    pre = Preconditioner(a, \"bddc\", coarsetype=\"h1amg\")        \n",
    "    a.Assemble()\n",
    "    f.Assemble()\n",
    "    \n",
    "    gfu = GridFunction(fes)\n",
    "    solvers.CG(mat=a.mat, rhs=f.vec, sol=gfu.vec,\n",
    "               pre=pre, maxsteps=500)\n",
    "Draw(gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Postscript\n",
    "\n",
    "By popular demand, here is the code to draw the figure found at the beginning of this tutorial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.geom2d import unit_square\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))\n",
    "fes_ho = Discontinuous(H1(mesh, order=10))\n",
    "fes_lo = H1(mesh, order=1, dirichlet=\".*\")\n",
    "fes_lam = Discontinuous(H1(mesh, order=1))\n",
    "fes = fes_ho*fes_lo*fes_lam\n",
    "uho, ulo, lam = fes.TrialFunction()\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += Variation(0.5 * grad(uho)*grad(uho)*dx \n",
    "               - 1*uho*dx \n",
    "               + (uho-ulo)*lam*dx(element_vb=BBND))\n",
    "gfu = GridFunction(fes)\n",
    "solvers.Newton(a=a, u=gfu)\n",
    "Draw(gfu.components[0],deformation=True, settings={\"camera\": {\"transformations\": [{\"type\": \"rotateX\", \"angle\": -45}]}})"
   ]
  }
 ],
 "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
}