{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 5.5.2 Discontinuous Galerkin for the Wave Equation\n",
    "\n",
    "We solve the first order wave equation by a matrix-free explicit DG method:\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\frac{\\partial p}{\\partial t} & = & \\operatorname{div} u \\\\\n",
    "\\frac{\\partial u}{\\partial t} & = & \\nabla p\n",
    "\\end{eqnarray*}\n",
    "\n",
    "\n",
    "Using DG discretization in space we obtain the ODE system\n",
    "\\begin{eqnarray*}\n",
    "M_p \\dot{p} & = & -B^T u \\\\\n",
    "M_u \\dot{u} & = & B p,\n",
    "\\end{eqnarray*}\n",
    "\n",
    "with mass-matrices $M_p$ and $M_u$, and the discretization matrix $B$ for the gradient, and $-B^T$ for the divergence. \n",
    "\n",
    "\n",
    "The DG bilinear-form for the gradient is:\n",
    "\n",
    "$$\n",
    "b(p,v) = \\sum_{T}\n",
    "\\Big\\{ \\int_T \\nabla p  \\, v + \\int_{\\partial T} (\\{ p \\} - p) \\, v_n \\, ds \\Big\\},\n",
    "$$\n",
    "where $\\{ p \\}$ is the average of $p$ on the facet.\n",
    "\n",
    "The computational advantages of a proper version of DG is:\n",
    "\n",
    "* universal element-matrices for the differntial operator $B$\n",
    "* cheaply invertible mass matrices (orthogonal polynomials, sum-factorization)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.occ import *\n",
    "from time import time\n",
    "\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.webgui import Draw as DrawGeo\n",
    "\n",
    "box = Box((-1,-1,-1), (1,1,0))\n",
    "sp = Sphere((0.5,0,0), 0.2)\n",
    "shape = box-sp\n",
    "geo = OCCGeometry(shape)\n",
    "\n",
    "h = 0.1\n",
    "        \n",
    "mesh = Mesh( geo.GenerateMesh(maxh=h))\n",
    "mesh.Curve(3)\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 4\n",
    "fes_p = L2(mesh, order=order, all_dofs_together=True)\n",
    "fes_u = VectorL2(mesh, order=order, piola=True)\n",
    "fes_tr = FacetFESpace(mesh, order=order)\n",
    "\n",
    "print (\"ndof_p = \", fes_p.ndof, \"+\", fes_tr.ndof, \", ndof_u =\", fes_u.ndof)\n",
    "\n",
    "traceop = fes_p.TraceOperator(fes_tr, average=True) \n",
    "\n",
    "gfu = GridFunction(fes_u)\n",
    "gfp = GridFunction(fes_p)\n",
    "gftr = GridFunction(fes_tr)\n",
    "\n",
    "gfp.Interpolate( exp(-400*(x**2+y**2+z**2)))\n",
    "gftr.vec.data = traceop * gfp.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = fes_p.TrialFunction()\n",
    "v = fes_u.TestFunction()\n",
    "phat = fes_tr.TrialFunction()\n",
    "\n",
    "n = specialcf.normal(mesh.dim)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "$$\n",
    "b(p,v) = \\sum_{T}\n",
    "\\Big\\{ \\int_T \\nabla p  \\, v + \\int_{\\partial T} (\\{ p \\} - p) \\, v_n \\, ds \\Big\\},\n",
    "$$\n",
    "\n",
    "where $\\{ p \\}$ is the average of $p$ on the facet.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "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",
    "Bel.Assemble()\n",
    "\n",
    "Btr = BilinearForm(trialspace=fes_tr, testspace=fes_u, geom_free = True)\n",
    "Btr += phat * (v*n) *dx(element_boundary=True)\n",
    "Btr.Assemble();\n",
    "\n",
    "B = Bel.mat + Btr.mat @ traceop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "invmassp = fes_p.Mass(1).Inverse()\n",
    "invmassu = fes_u.Mass(1).Inverse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfp.Interpolate( exp(-100*(x**2+y**2+z**2)))\n",
    "gfu.vec[:] = 0\n",
    "\n",
    "scene = Draw (gftr, draw_vol=False, order=3, min=-0.05, max=0.05);\n",
    "\n",
    "t = 0\n",
    "dt = 0.3 * h / (order+1)**2 \n",
    "tend = 1\n",
    "\n",
    "op1 = dt * invmassu @ B\n",
    "op2 = dt * invmassp @ B.T\n",
    "\n",
    "cnt = 0\n",
    "with TaskManager(): \n",
    "    while t < tend:\n",
    "        t = t+dt\n",
    "        gfu.vec.data += op1 * gfp.vec\n",
    "        gfp.vec.data -= op2 * gfu.vec\n",
    "        cnt = cnt+1\n",
    "        if cnt%10 == 0:\n",
    "            gftr.vec.data = traceop * gfp.vec\n",
    "            scene.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## Time-stepping on the device"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import ngsolve.ngscuda\n",
    "except:\n",
    "    print (\"Sorry, no cuda device\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfp.Interpolate( exp(-100*(x**2+y**2+z**2)))\n",
    "gfu.vec[:] = 0\n",
    "\n",
    "scene = Draw (gftr, draw_vol=False, order=3);\n",
    "\n",
    "t = 0\n",
    "dt = 0.5 * h / (order+1)**2 / 2\n",
    "tend = 0.1\n",
    "\n",
    "op1 = (dt * invmassu @ B).CreateDeviceMatrix()\n",
    "op2 = (dt * invmassp @ B.T).CreateDeviceMatrix()\n",
    "\n",
    "devu = gfu.vec.CreateDeviceVector(copy=True)\n",
    "devp = gfp.vec.CreateDeviceVector(copy=True)\n",
    "\n",
    "cnt = 0\n",
    "with TaskManager(): \n",
    "    while t < tend:\n",
    "        t = t+dt\n",
    "        devu += op1 * devp\n",
    "        devp -= op2 * devu\n",
    "        cnt = cnt+1\n",
    "        if cnt%10 == 0:\n",
    "            gfp.vec.data = devp\n",
    "            gftr.vec.data = traceop * gfp.vec\n",
    "            scene.Redraw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "print (op1.GetOperatorInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = time()\n",
    "steps = 10\n",
    "for i in range(steps):\n",
    "    devu += op1 * devp\n",
    "    devp -= op2 * devu\n",
    "te = time()\n",
    "print (\"ndof = \", gfp.space.ndof, \"+\", gfu.space.ndof, \", time per step =\", (te-ts)/steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "On the A-100:"
   ]
  },
  {
   "cell_type": "raw",
   "id": "14",
   "metadata": {},
   "source": [
    "ndof =  651035 + 1953105 , time per step = 0.0023579835891723634"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {},
   "source": [
    "## Time-domain PML\n",
    "\n",
    "PML (perfectly matched layers) is a method for approximating outgoing waves on a truncated domain. Its time domain version leads to additional field variables in the PML region, which are coupled via time-derivatives only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ring_resonator_import import *\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes)\n",
    "scene = Draw (gfu.components[0], order=3, min=-0.05, max=0.05, autoscale=False)\n",
    "\n",
    "t = 0\n",
    "tend = 15\n",
    "tau = 2e-4\n",
    "i = 0\n",
    "sigma = 10   # pml damping parameter\n",
    "\n",
    "op1 = invp@(-fullB.T-sigma*dampingp) \n",
    "op2 = invu@(fullB-sigma*dampingu)\n",
    "with TaskManager(): \n",
    "    while t < tend:\n",
    "\n",
    "        gfu.vec.data += tau*Envelope(t)*srcvec\n",
    "        gfu.vec.data += tau*op1*gfu.vec\n",
    "        gfu.vec.data += tau*op2*gfu.vec        \n",
    "\n",
    "        t += tau\n",
    "        i += 1\n",
    "        if i%20 == 0:\n",
    "            scene.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18",
   "metadata": {},
   "source": [
    "The differential operators and coupling terms to the pml - variables are represented via linear operators;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19",
   "metadata": {},
   "outputs": [],
   "source": [
    "print (fullB.GetOperatorInfo())\n",
    "print (dampingp.GetOperatorInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20",
   "metadata": {},
   "outputs": [],
   "source": [
    "print (op1.GetOperatorInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "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": 5
}