{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.8 Discontinuous Galerkin Methods\n",
    "\n",
    "* Use discontinuous finite element spaces to solve PDEs. \n",
    "* Allows upwind-stabilization for convection-dominated problems\n",
    "* Requires additional jump terms for consistency \n",
    "\n",
    "Interior penalty DG form for $-\\Delta u$:\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "A(u,v) = \\sum_T \\int_T \\nabla u \\nabla v\n",
    "-  \\sum_F \\int_F \\{ n \\nabla u \\} [v] \n",
    "-  \\sum_F \\int_F \\{ n \\nabla v \\} [u] \n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F [u][v]\n",
    "$$\n",
    "\n",
    "with jump-term over facets:\n",
    "$$\n",
    "[u] = u_{left} - u_{right}\n",
    "$$\n",
    "\n",
    "and averaging operator\n",
    "$$\n",
    "\\{ n \\nabla u \\} = \\tfrac{1}{2} (n_{left} \\nabla u_{left} + n_{left} \\nabla u_{right})\n",
    "$$\n",
    "\n",
    "DG form for $\\Div (b u)$, where $b$ is the given wind:\n",
    "\n",
    "$$\n",
    "B(u,v) = -\\sum_T b u \\nabla v + \\sum_F \\int_F b\\cdot n   u^{upwind} v \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The space is responsible for allocating the matrix graph. Tell it that it should reserve entries for the coupling terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order=4\n",
    "fes = L2(mesh, order=order, dgjumps=True)\n",
    "u,v = fes.TnT()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every facet has a master element. The value from the other element is referred to via the\n",
    "`Other()` operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_u = u-u.Other()\n",
    "jump_v = v-v.Other()\n",
    "n = specialcf.normal(2)\n",
    "mean_dudn = 0.5*n * (grad(u)+grad(u.Other()))\n",
    "mean_dvdn = 0.5*n * (grad(v)+grad(v.Other()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Integrals on facets are computed by setting `skeleton=True`. \n",
    "* `dx(skeleton=True)` iterates over all internal faces\n",
    "* `ds(skeleton=True)` iterates over all boundary faces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 4\n",
    "h = specialcf.mesh_size\n",
    "diffusion = grad(u)*grad(v) * dx \\\n",
    "    +alpha*order**2/h*jump_u*jump_v * dx(skeleton=True) \\\n",
    "    +(-mean_dudn*jump_v-mean_dvdn*jump_u) * dx(skeleton=True) \\\n",
    "    +alpha*order**2/h*u*v * ds(skeleton=True) \\\n",
    "    +(-n*grad(u)*v-n*grad(v)*u)* ds(skeleton=True)\n",
    "\n",
    "a = BilinearForm(diffusion).Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = LinearForm(1*v*dx).Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes, name=\"uDG\")\n",
    "gfu.vec.data = a.mat.Inverse() * f.vec\n",
    "Draw (gfu);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DG requires a lot of additional matrix entries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes2 = L2(mesh, order=order)\n",
    "ul2,vl2 = fes2.TnT()\n",
    "a2 = BilinearForm(ul2*vl2*dx).Assemble()\n",
    "print (\"DG-matrix nze:\", a.mat.nze)\n",
    "print (\"L2-matrix nze:\", a2.mat.nze)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we are solving a convection-diffusion problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 4\n",
    "h = specialcf.mesh_size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `IfPos` checks whether the first argument is positive. Then it returns the second one, else the third one. This is used to define the upwind flux. The check is performed in every integration-point on the skeleton:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = CF( (20,5) )\n",
    "uup = IfPos(b*n, u, u.Other())\n",
    "\n",
    "convection = -b * u * grad(v)*dx + b*n*uup*jump_v * dx(skeleton=True)\n",
    "\n",
    "acd = BilinearForm(diffusion + convection).Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = acd.mat.Inverse(freedofs=fes.FreeDofs(),inverse=\"umfpack\") * f.vec\n",
    "Draw (gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hybrid Discontinuous Galerkin methods\n",
    "use additionally the *hybrid* facet variable on the skeleton:\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "A(u,\\widehat u; v, \\widehat v) = \n",
    "  \\sum_T \\int_T \\nabla u \\nabla v\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla u (v-\\widehat v)\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla v (u-\\widehat u)\n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F (u-\\widehat u)(v-\\widehat v)\n",
    "$$\n",
    "\n",
    "the jump-term is now replaced by the difference $u - \\widehat u$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No additional matrix entries across elements are produced. Dirichlet boundary conditions are set as usual to the facet variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order=4\n",
    "V = L2(mesh, order=order)\n",
    "F = FacetFESpace(mesh, order=order, dirichlet=\"bottom|left|right|top\")\n",
    "fes = V*F\n",
    "u,uhat = fes.TrialFunction()\n",
    "v,vhat = fes.TestFunction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the jump is the difference between element-term and facet-term:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jump_u = u-uhat\n",
    "jump_v = v-vhat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 4\n",
    "condense = True\n",
    "h = specialcf.mesh_size\n",
    "n = specialcf.normal(mesh.dim)\n",
    "\n",
    "a = BilinearForm(fes, condense=condense)\n",
    "dS = dx(element_boundary=True)\n",
    "a += grad(u)*grad(v)*dx + \\\n",
    "    alpha*order**2/h*jump_u*jump_v*dS + \\\n",
    "    (-grad(u)*n*jump_v - grad(v)*n*jump_u)*dS\n",
    "\n",
    "b = CF( (20,1) )\n",
    "uup = IfPos(b*n, u, uhat)\n",
    "a += -b * u * grad(v)*dx + b*n*uup*jump_v *dS\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(1*v*dx).Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "print (\"A non-zero elements:\", a.mat.nze)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not condense:\n",
    "    inv = a.mat.Inverse(fes.FreeDofs(), \"umfpack\")\n",
    "    gfu.vec.data = inv * f.vec\n",
    "else:\n",
    "    fmod = (f.vec + a.harmonic_extension_trans * f.vec).Evaluate()\n",
    "    \n",
    "    inv = a.mat.Inverse(fes.FreeDofs(True), \"umfpack\")\n",
    "    gfu.vec.data = inv * fmod\n",
    "    \n",
    "    gfu.vec.data += a.harmonic_extension * gfu.vec\n",
    "    gfu.vec.data += a.inner_solve * f.vec\n",
    "\n",
    "Draw (gfu.components[0], mesh, \"u-HDG\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Projected Jumps\n",
    "The polynomial order of the facet-space can be reduced by one by inserting an $L_2$-projection into the penalty term. This can be implemented by keeping the highest order basis functions of the facet space discontinuous, see [Lehrenfeld-Schöberl: High order exactly divergence-free Hybrid Discontinuous Galerkin Methods for unsteady incompressible flows](https://www.sciencedirect.com/science/article/pii/S004578251630264X?via%3Dihub):\n",
    "\n",
    "$$\n",
    "\\DeclareMathOperator{\\Div}{div}\n",
    "A(u,\\widehat u; v, \\widehat v) = \n",
    "  \\sum_T \\int_T \\nabla u \\nabla v\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla u (v-\\widehat v)\n",
    "- \\sum_T \\int_{\\partial T} n \\nabla v (u-\\widehat u)\n",
    "+ \\frac{\\alpha p^2}{h} \\sum_F \\int_F (\\Pi u-\\widehat u)(\\Pi v-\\widehat v)\n",
    "$$\n",
    "\n",
    "This is achieved by setting the flag `highest_order_dc=True` in the `FacetFESpace`. \n",
    "The number of matrix entries is reduced, while the order of convergence is preserved.\n",
    "This trick does not work in the case of significant convection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "order=4\n",
    "V = L2(mesh, order=order)\n",
    "F = FacetFESpace(mesh, order=order, dirichlet=\"bottom|left|right|top\", \\\n",
    "                          highest_order_dc=True)\n",
    "fes = V*F\n",
    "u,uhat = fes.TrialFunction()\n",
    "v,vhat = fes.TestFunction()\n",
    "\n",
    "jump_u = u-uhat\n",
    "jump_v = v-vhat\n",
    "\n",
    "alpha = 2\n",
    "h = specialcf.mesh_size\n",
    "n = specialcf.normal(mesh.dim)\n",
    "\n",
    "a = BilinearForm(fes, condense=True)\n",
    "dS = dx(element_boundary=True)\n",
    "a += grad(u)*grad(v)*dx + \\\n",
    "    alpha*(order+1)**2/h*jump_u*jump_v*dS + \\\n",
    "    (-grad(u)*n*jump_v - grad(v)*n*jump_u)*dS\n",
    "a.Assemble()\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f += 1*v*dx\n",
    "f.Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "\n",
    "f.vec.data += a.harmonic_extension_trans * f.vec \n",
    "    \n",
    "inv = a.mat.Inverse(fes.FreeDofs(True), \"sparsecholesky\")\n",
    "gfu.vec.data = inv * f.vec\n",
    "    \n",
    "gfu.vec.data += a.harmonic_extension * gfu.vec\n",
    "gfu.vec.data += a.inner_solve * f.vec\n",
    "\n",
    "Draw (gfu.components[0], mesh, \"u-HDG\")\n",
    "\n",
    "print (\"A non-zero elements:\", a.mat.nze)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remarks on sparsity pattern in NGSolve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Remark 1: The sparsity pattern is set up a-priorily\n",
    "* The sparsity pattern of a sparse matrix in NGSolve is independent of its entries (it's set up a-priorily). \n",
    "* We can have \"nonzero\" entries that have the value 0\n",
    "\n",
    "Below we show the reserved memory for the sparse matrix and the (numerically) non-zero entries in this sparse matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes2 = L2(mesh, order=order, dgjumps=True)\n",
    "u,v=fes2.TnT()\n",
    "a3 = BilinearForm(fes2)\n",
    "a3 += u*v*dx + (u+u.Other())*v*dx(skeleton=True)\n",
    "a3.Assemble();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse as sp\n",
    "import matplotlib.pylab as plt\n",
    "plt.rcParams['figure.figsize'] = (12, 12)\n",
    "A = sp.csr_matrix(a3.mat.CSR())\n",
    "fig = plt.figure(); ax1 = fig.add_subplot(121); ax2 = fig.add_subplot(122)\n",
    "ax1.set_xlabel(\"numerically non-zero\"); ax1.spy(A)\n",
    "ax2.set_xlabel(\"reserved entries (potentially non-zero)\"); ax2.spy(A,precision=-1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Remark 2: Sparsity pattern with and without `dgjumps=True` is different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a1 = BilinearForm(L2(mesh, order=order, dgjumps=False)); a1.Assemble()\n",
    "a2 = BilinearForm(L2(mesh, order=order, dgjumps=True)); a2.Assemble()\n",
    "A1 = sp.csr_matrix(a1.mat.CSR())\n",
    "A2 = sp.csr_matrix(a2.mat.CSR())\n",
    "fig = plt.figure(); ax1 = fig.add_subplot(121); ax2 = fig.add_subplot(122)\n",
    "ax1.set_xlabel(\"dgjumps=False\"); ax1.spy(A1,precision=-1)\n",
    "ax2.set_xlabel(\"dgjumps=True\"); ax2.spy(A2,precision=-1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Remark 3: Dof numbering of higher order FESpaces "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "* In `NGSolve` `FESpace`s typically have a numbering where the first block of dofs corresponds to a low order subspace (which is convenient for iterative solvers). \n",
    "* For L2 this means that the first dofs correspond to the constants on elements. \n",
    "\n",
    "* You can turn this behavior off for some spaces, e.g. for L2 by adding the flag `all_dofs_together`.\n",
    "\n",
    "We demonstrate this in the next comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (15, 15)\n",
    "fig = plt.figure()\n",
    "ax = [fig.add_subplot(131), fig.add_subplot(132), fig.add_subplot(133)]\n",
    "\n",
    "for i, (order, all_dofs_together, label) in enumerate([(0,False, \"non-zeros (p=0)\"),\n",
    "                                                    (1,False,\"non-zeros (p=1, low order + high order)\"),\n",
    "                                                    (1,True,\"non-zeros (p=1, all_dofs_together)\")]):\n",
    "    a = BilinearForm(L2(mesh,order=order,dgjumps=True,all_dofs_together=all_dofs_together))\n",
    "    a.Assemble()\n",
    "    ax[i].spy(sp.csr_matrix(a.mat.CSR()),markersize=3,precision=-1)\n",
    "    ax[i].set_xlabel(label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}