{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 3.5.2 DG for the acoustic wave equation\n",
    "We consider the acoustic problem\n",
    "\\begin{align} \n",
    "\\partial_{t} \\mathbf{u}  -  c \\nabla p &= 0 \\quad \\text{ in } \\Omega \\times I, \\\\\n",
    "\\partial_{t} p  -  c \\operatorname{div}(\\mathbf{u}) &= 0 \\quad \\text{ in } \\Omega \\times I, \\\\\n",
    "%      p (0,\\cdot) &= p(1,\\cdot), \\quad\n",
    "%      p (\\cdot,0) = p(\\cdot,1), \\label{eq:per2b}\\\\\n",
    "%      \\mathbf{q} (0,\\cdot) &= \\mathbf{q}(1,\\cdot), \\quad\n",
    "%      \\mathbf{q} (\\cdot,0) = \\mathbf{q}(\\cdot,1), \\label{eq:per2bb}\\\\\n",
    "p &= p_0 \\!\\!\\! \\quad \\text{ on } \\Omega \\times \\{0\\},\\\\\n",
    "u &= 0 \\quad \\text{ on } \\Omega \\times \\{0\\}.\n",
    "\\end{align}\n",
    "\n",
    "Here $p$ is the acoustic pressure (the local deviation from the ambient pressure) and $\\mathbf{u}$ is the local velocity.\n",
    "\n",
    "`+` suitable boundary conditions "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A simple grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#imports, geometry and mesh:\n",
    "from netgen.occ import *\n",
    "from netgen.webgui import Draw as DrawGeo\n",
    "geo = OCCGeometry(unit_square_shape.Scale((0,0,0),2).Move((-1,-1,0)), dim=2)\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.1))\n",
    "Draw(mesh)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Find $p: [0,T] \\to \\bigoplus_{T\\in\\mathcal{T}_h} \\mathcal{P}^{k+1}(T)$ and $\\mathbf{u}: [0,T] \\to \\bigoplus_{T\\in\\mathcal{T}_h} [\\mathcal{P}^k(T)]^d$ so that \n",
    "\n",
    "\n",
    "\\begin{align*}\n",
    "(\\partial_t \\mathbf{u}, v) &= b_h(p,v) & &&& \\forall v,\\\\\n",
    "(\\partial_t p, q) &= & - b_h(q,\\mathbf{u}) &&& \\forall q\\\\    \n",
    "\\end{align*}\n",
    "\n",
    "with the centered flux approximation:\n",
    "$$\n",
    "b_h(p,v) = \\sum_{T} \\int_T \\nabla p \\cdot v + \\int_{\\partial T} ( \\hat{p} - p) v_n \n",
    "$$\n",
    "\n",
    "Here $\\hat{p}$ is the centered approximation i.e. $\\hat{p} = \\{\\!\\!\\{p\\}\\!\\!\\}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "order = 6\n",
    "fes_p = L2(mesh, order=order+1, all_dofs_together=True)\n",
    "fes_u = VectorL2(mesh, order=order, piola=True, all_dofs_together=True)\n",
    "fes = fes_p*fes_u\n",
    "gfu = GridFunction(fes)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### What is the flag `piola` doing?\n",
    "\n",
    "The `VectorL2` space uses the following definition of basis functions on the mesh. \n",
    "\n",
    "Let $\\hat{\\varphi}(\\hat{x})$ be a (vectorial) basis function on the reference element $\\hat{T}$, $\\Phi: \\hat{T} \\to T$ with $x = \\Phi(\\hat{x})$ and $F= D\\Phi$, then\n",
    "\n",
    "$$\n",
    "  \\varphi(x) := \\left\\{ \\begin{array}{rlll} \n",
    "  \\hat{\\varphi}(\\hat{x}), & \\texttt{piola=False}, & \\texttt{covariant=False}, & \\text{(default)} \\\\\n",
    "  \\frac{1}{|\\operatorname{det} F|} \\cdot F \\cdot \\hat{\\varphi}(\\hat{x}), & \\texttt{piola=True},  & \\texttt{covariant=False}, \\\\\n",
    "  F^{-T} \\cdot \\hat{\\varphi}(\\hat{x}), & \\texttt{piola=False}, & \\texttt{covariant=True}. \\\\\n",
    "  \\end{array}\n",
    "  \\right.\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Inverse mass matrix operations:\n",
    "Combining `Embedding` and the inverse mass matrix operation for `fes_p` allows to realize the following block inverse operations acting on `fes`:"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "$$\n",
    "\\underbrace{\\left( \\begin{array}{cc} M_{p}^{-1} & 0 \\\\ 0 & 0 \\end{array} \\right)}_{\\texttt{invp}} = \n",
    "\\underbrace{\\left( \\begin{array}{c} I \\\\ 0 \\end{array} \\right)}_{\\texttt{emb_p}} \\cdot \\underbrace{M_{p}^{-1}}_{\\texttt{invmassp}} \\cdot \\underbrace{\\left( \\begin{array}{cc} I & 0 \\end{array} \\right)}_{\\texttt{emb_p.T}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "pdofs = fes.Range(0);\n",
    "emb_p = Embedding(fes.ndof, pdofs)\n",
    "invmassp = fes_p.Mass(1).Inverse()\n",
    "invp = emb_p @ invmassp @ emb_p.T"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Analogously, combining `Embedding` and the inverse mass matrix operation for `fes_u` allows to realize the following block inverse operations on `fes`:"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "$$\n",
    "\\underbrace{\\left( \\begin{array}{cc} 0 & 0 \\\\ 0 & M_{\\mathbf{u}}^{-1} \\end{array} \\right)}_{\\texttt{invu}} = \n",
    "\\underbrace{\\left( \\begin{array}{c} I \\\\ 0 \\end{array} \\right)}_{\\texttt{emb_u}} \\cdot \\underbrace{M_{\\mathbf{u}}^{-1}}_{\\texttt{invmassu}} \\cdot \\underbrace{\\left( \\begin{array}{cc} I & 0 \\end{array} \\right)}_{\\texttt{emb_u.T}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "udofs = fes.Range(1)\n",
    "emb_u = Embedding(fes.ndof, udofs)\n",
    "invmassu = fes_u.Mass(Id(mesh.dim)).Inverse()\n",
    "invu = emb_u @ invmassu @ emb_u.T"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Time loop \n",
    "Assuming the operators `B`,`BT`: `fes_p`$\\times$ `fes_u` $\\to$(`fes_p`$\\times$ `fes_u` )$'$\n",
    "corresponding to $B_h((u,p),(v,q)=b_h(p,v)$ are given, the symplectic Euler time stepping method takes the form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def Run(B, BT, t0 = 0, tend = 0.25, tau = 1e-3, \n",
    "        backward = False, scenes = []):\n",
    "    t = 0\n",
    "    with TaskManager():\n",
    "        while t < (tend-t0) - tau/2:\n",
    "            t += tau\n",
    "            if not backward:\n",
    "                gfu.vec.data += -tau * invp @ BT * gfu.vec\n",
    "                gfu.vec.data += tau * invu @ B * gfu.vec\n",
    "                print(\"\\r t = {:}\".format(t0 + t),end=\"\")\n",
    "            else:\n",
    "                gfu.vec.data += -tau * invu @ B * gfu.vec\n",
    "                gfu.vec.data += tau * invp @ BT * gfu.vec\n",
    "                print(\"\\r t = {:}\".format(tend - t),end=\"\")\n",
    "            for sc in scenes: sc.Redraw() # blocking=False)\n",
    "        print(\"\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Initial values (density ring):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "gfu.components[0].Set (exp(-50*(x**2+y**2))-exp(-100*(x**2+y**2)))\n",
    "gfu.components[1].vec[:] = 0\n",
    "Draw(gfu.components[0], mesh, \"p\", min=-0.02, max=0.08, autoscale=False, order=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Version 1:\n",
    "## The bilinear form for application on the full space `fes`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "n = specialcf.normal(mesh.dim) \n",
    "(p,u),(q,v) = fes.TnT()\n",
    "B = BilinearForm(fes, nonassemble=True)\n",
    "B += grad(p)*v * dx + 0.5*(p.Other()-p)*(v*n) * dx(element_boundary=True)\n",
    "BT = BilinearForm(fes, nonassemble=True)\n",
    "BT += grad(q)*u * dx + 0.5*(q.Other()-q)*(u*n) * dx(element_boundary=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "scenep=Draw(gfu.components[0], mesh, \"p\", min=-0.02, max=0.08, \n",
    "            autoscale=False, order=3, deformation=True) \n",
    "%time Run(B.mat, BT.mat, backward=False, scenes=[scenep])\n",
    "%time Run(B.mat, BT.mat, backward=True, scenes=[scenep])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Version 2: \n",
    "## Using the `TraceOperator` and assembling of $b_h$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We want to use the `TraceOperator` again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fes_tr = FacetFESpace(mesh, order=order+1)\n",
    "traceop = fes_p.TraceOperator(fes_tr, False) "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to assemble the sub-block matrix and need test/trial functions for single spaces (not the product spaces):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "p = fes_p.TrialFunction()\n",
    "v = fes_u.TestFunction()\n",
    "phat = fes_tr.TrialFunction()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We split the operator to $b_h$ into\n",
    "* volume contributions (local) \n",
    "* and couplings between the trace (obtained through the trace op) and the volume:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Bel = BilinearForm(trialspace=fes_p, testspace=fes_u)\n",
    "Bel += grad(p)*v * dx -p*(v*n) * dx(element_boundary=True)\n",
    "%time Bel.Assemble()\n",
    "\n",
    "Btr = BilinearForm(trialspace=fes_tr, testspace=fes_u)\n",
    "Btr += 0.5 * phat * (v*n) * dx(element_boundary=True) \n",
    "%time Btr.Assemble()\n",
    "\n",
    "B = emb_u @ (Bel.mat + Btr.mat @ traceop) @ emb_p.T"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Version 2 in action (assembled matrices, `TraceOperators` and `Embeddings` can do transposed multiply):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%time Run(B, B.T, backward=False)\n",
    "%time Run(B, B.T, backward=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Version 3: \n",
    "## Using the `TraceOperator` and `geom_free=True` with assembling"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With $p(x) = \\hat{p}(\\hat{x})$ and $\\mathbf{v}(x) = \\frac{1}{|\\operatorname{det}(F)|} F \\cdot \\hat{\\mathbf{v}}(\\hat{x})$ (Piola mapped) there holds\n",
    "\n",
    "\\begin{align}\n",
    "\\int_T \\mathbf{v} \\cdot \\nabla p  = \\int_{\\hat{T}} \\hat{\\mathbf{v}} \\cdot \\hat{\\nabla} \\hat{p}\n",
    "\\end{align}\n",
    "\n",
    "and similarly for the facet integrals, cf. [unit-2.11](../unit-2.11-matrixfree/matrixfree.ipynb) and \n",
    "\n",
    "Integrals are independent of the \"physical element\" (up to ordering) \n",
    "\n",
    "$\\leadsto$ same element matrices for a large class of elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "Bel = BilinearForm(trialspace=fes_p, testspace=fes_u, geom_free = True)\n",
    "Bel += grad(p)*v * dx -p*(v*n) * dx(element_boundary=True)\n",
    "%time Bel.Assemble()\n",
    "\n",
    "Btr = BilinearForm(trialspace=fes_tr, testspace=fes_u, geom_free = True)\n",
    "Btr += 0.5 * phat * (v*n) *dx(element_boundary=True)\n",
    "%time Btr.Assemble()\n",
    "\n",
    "B = emb_u @ (Bel.mat + Btr.mat @ traceop) @ emb_p.T"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Version 3 in action:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%time Run(B, B.T, backward=False)\n",
    "%time Run(B, B.T, backward=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "References:\n",
    "\n",
    "* [Hesthaven, Warburton. Nodal Discontinuous Galerkin Methods—Algorithms, Analysis and Applications, 2007]\n",
    "* [[Koutschan, Lehrenfeld, Schöberl. Computer algebra meets finite elements: An efficient implementation for Maxwells equations. , 2012](https://arxiv.org/pdf/1104.4208v2.pdf)]\n"
   ]
  }
 ],
 "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.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}