NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2fd77000",
   "metadata": {},
   "source": [
    "# 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": "8b11665a",
   "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",
    "disc = Cylinder( Axes((0.5,0,0), X, Y), r=0.4, h=0.2)\n",
    "shape = box-disc\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": "3dbda538",
   "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": "1c19b6ce",
   "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": "835e8a51",
   "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": "7903ef8f",
   "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": "3036f93d",
   "metadata": {},
   "outputs": [],
   "source": [
    "invmassp = fes_p.Mass(1).Inverse()\n",
    "invmassu = fes_u.Mass(1).Inverse()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d49f6cff",
   "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.3 * h / (order+1)**2 \n",
    "tend = 1.2\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": "719e9c5e",
   "metadata": {},
   "source": [
    "## Time-stepping on the device"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca806f69",
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import ngsolve.ngscuda\n",
    "except:\n",
    "    print (\"Sorry, no cuda device\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33406cbb",
   "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": "77281b6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print (op1.GetOperatorInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "575aa7b0",
   "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": "93f8f355",
   "metadata": {},
   "source": [
    "On the A-100:\n",
    "\n",
    "    ndof =  651035 + 1953105 , time per step = 0.0023579835891723634"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c176e930",
   "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.\n",
    "\n",
    "Formulation of B. Kapidani, J. Schöberl, [arxiv](https://arxiv.org/abs/2002.08733)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52a88f41",
   "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": "362613ed-ec18-40f6-87c9-3c12e437e017",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.geom2d import SplineGeometry\n",
    "\n",
    "geo = SplineGeometry()\n",
    "\n",
    "xneg  =-0.43\n",
    "xpos  = 0.43\n",
    "yneg  =-0.48\n",
    "ypos  = 0.48\n",
    "wslab = 0.04\n",
    "cringx = 0.0\n",
    "cringy = 0.0\n",
    "rring = 0.4\n",
    "gap   = 0.005\n",
    "\n",
    "pntx = [xneg,xpos]\n",
    "pnty = [yneg,-rring-gap-wslab,-rring-gap,rring+gap,rring+gap+wslab,ypos]\n",
    "\n",
    "pts = [geo.AddPoint(xi,yi) for yi in pnty for xi in pntx]\n",
    "\n",
    "### geometry #######################################################\n",
    "geo.Append ([\"line\", pts[0], pts[1]], leftdomain=1, rightdomain=0)\n",
    "geo.Append ([\"line\", pts[1], pts[3]], leftdomain=1, rightdomain=0)\n",
    "geo.Append ([\"line\", pts[3], pts[2]], leftdomain=1, rightdomain=2)\n",
    "geo.Append ([\"line\", pts[2], pts[0]], leftdomain=1, rightdomain=0)\n",
    "\n",
    "geo.Append ([\"line\", pts[3], pts[5]], leftdomain=2, rightdomain=0,bc=\"normal_wg_rightbottom\")\n",
    "geo.Append ([\"line\", pts[5], pts[4]], leftdomain=2, rightdomain=3)\n",
    "geo.Append ([\"line\", pts[4], pts[2]], leftdomain=2, rightdomain=0,bc=\"normal_wg_leftbottom\")\n",
    "\n",
    "geo.Append ([\"line\", pts[5], pts[7]], leftdomain=3, rightdomain=0)\n",
    "geo.Append ([\"line\", pts[7], pts[6]], leftdomain=3, rightdomain=4)\n",
    "geo.Append ([\"line\", pts[6], pts[4]], leftdomain=3, rightdomain=0)\n",
    "\n",
    "geo.Append ([\"line\", pts[7], pts[9]], leftdomain=4, rightdomain=0,bc=\"normal_wg_righttop\")\n",
    "geo.Append ([\"line\", pts[9], pts[8]], leftdomain=4, rightdomain=5)\n",
    "geo.Append ([\"line\", pts[8], pts[6]], leftdomain=4, rightdomain=0,bc=\"normal_wg_lefttop\")\n",
    "\n",
    "geo.Append ([\"line\", pts[9], pts[11]], leftdomain=5, rightdomain=0)\n",
    "geo.Append ([\"line\", pts[11], pts[10]], leftdomain=5, rightdomain=0)\n",
    "geo.Append ([\"line\", pts[10], pts[8]], leftdomain=5, rightdomain=0)\n",
    "\n",
    "geo.AddCircle(c=(cringx,cringy), r=rring, leftdomain=6, rightdomain=3)\n",
    "geo.AddCircle(c=(cringx,cringy), r=rring-wslab, leftdomain=7, rightdomain=6)\n",
    "\n",
    "for i in (1,3,5,7):\n",
    "    geo.SetMaterial(i, \"air\")\n",
    "for i in (2,4,6):\n",
    "    geo.SetMaterial(i, \"eps_nine\")\n",
    "data = geo.CreatePML(0.05)\n",
    "normals = data[\"normals\"]\n",
    "\n",
    "\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))\n",
    "\n",
    "eps_r = {\"air\" : 1, \"eps_nine\" : 3**3}\n",
    "\n",
    "for mat in mesh.GetMaterials():\n",
    "    if mat.startswith(\"pml_normal_wg\"):\n",
    "        eps_r[mat] = eps_r[\"eps_nine\"]\n",
    "\n",
    "for mat in mesh.GetMaterials():\n",
    "    if mat not in eps_r:\n",
    "        eps_r[mat] = eps_r[\"air\"]\n",
    "\n",
    "### Parameters for Source field ######\n",
    "wavelength = 1.542\n",
    "fcen       = 5/wavelength\n",
    "df         = 0.1\n",
    "tpeak      = 1\n",
    "order = 3\n",
    "\n",
    "mesh.Curve(order)\n",
    "fes_facet = FacetFESpace(mesh, order=order+1)\n",
    "gfsource = GridFunction(fes_facet)\n",
    "\n",
    "source_cf =  sin( (pi/wslab)*(y-pnty[3]) ) \n",
    "gfsource.Set(source_cf,definedon=mesh.Boundaries(\"normal_wg_lefttop\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41a7b803",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes_u = VectorL2(mesh, order=order, piola=True, order_policy=ORDER_POLICY.CONSTANT)\n",
    "fes_p = L2(mesh, order=order+1, all_dofs_together=True, order_policy=ORDER_POLICY.CONSTANT)\n",
    "fes_tr = FacetFESpace(mesh, order=order+1)\n",
    "fes_hdiv = HDiv(mesh, order=order+1, orderinner=1)\n",
    "\n",
    "n = specialcf.normal(2) \n",
    "\n",
    "p,q = fes_p.TnT()\n",
    "u,v = fes_u.TnT()\n",
    "ptr,qtr = fes_tr.TnT()\n",
    "uhdiv,vhdiv = fes_hdiv.TnT()\n",
    "\n",
    "Bel = BilinearForm(grad(p)*v*dx - p*(v*n)*dx(element_boundary=True), geom_free=True).Assemble()\n",
    "Btr = BilinearForm(0.5*ptr*(n*v)*dx(element_boundary=True), geom_free=True).Assemble()\n",
    "Bstab = BilinearForm(p*(vhdiv*n)*dx(element_boundary=True),  geom_free=True).Assemble()\n",
    "\n",
    "\n",
    "nvec = { mat : ((normals[mat][0], normals[mat][1]) if mat in normals else (0,0)) for mat in mesh.GetMaterials() }\n",
    "\n",
    "cfn = CF( [CF(nvec[mat]) for mat in mesh.GetMaterials()])\n",
    "cft = CF( ( cfn[1], -cfn[0] ) )\n",
    "\n",
    "pml1d = mesh.Materials(\"pml_default.*|pml_normal.*\")\n",
    "eps = CF([eps_r[mat] for mat in mesh.GetMaterials()])\n",
    "\n",
    "\n",
    "fes = fes_p*fes_p*fes_u*fes_u*fes_hdiv\n",
    "emb_p, emb_phat, emb_u, emb_uhat, emb_hdiv = fes.embeddings\n",
    "\n",
    "# gradient operator\n",
    "traceop = fes_p.TraceOperator(fes_tr, False)\n",
    "fullB = emb_u @ (Bel.mat + Btr.mat @ traceop) @ emb_p.T + emb_hdiv@Bstab.mat@emb_p.T\n",
    "\n",
    "# mass matrices\n",
    "invmassp = fes_p.Mass(eps).Inverse()\n",
    "invmassu = fes_u.Mass(Id(mesh.dim)).Inverse()\n",
    "Mstab = BilinearForm(uhdiv*vhdiv*dx, diagonal=True).Assemble()\n",
    "Mstabinv = Mstab.mat.Inverse()\n",
    "\n",
    "\n",
    "invp = emb_p @ invmassp @ emb_p.T + emb_phat @ invmassp @ emb_phat.T\n",
    "invu = emb_u @ invmassu @ emb_u.T + emb_uhat @ invmassu @ emb_uhat.T + emb_hdiv@Mstabinv@emb_hdiv.T\n",
    "\n",
    "\n",
    "# damping matrices\n",
    "dampp1 = fes_p.Mass (eps, definedon=pml1d)\n",
    "dampp2 = fes_p.Mass (eps, definedon=mesh.Materials(\"pml_corner\"))\n",
    "dampu1 = fes_u.Mass (OuterProduct(cfn,cfn), definedon=pml1d)\n",
    "dampu2 = fes_u.Mass (OuterProduct(cft,cft), definedon=pml1d)\n",
    "\n",
    "dampingu = emb_u @ dampu1 @ emb_u.T + (-emb_u + emb_uhat) @ dampu2 @ (emb_u.T + emb_uhat.T)\n",
    "dampingp = emb_p @ dampp1 @ emb_p.T + emb_p @ dampp2 @ (2*emb_p.T-emb_phat.T) + emb_phat @ dampp2 @ emb_p.T\n",
    "\n",
    "\n",
    "# source term\n",
    "Lsrc  = LinearForm(gfsource*q*dx(element_boundary=True)).Assemble()\n",
    "srcvec = emb_p * (invmassp*Lsrc.vec).Evaluate()\n",
    "\n",
    "# time-envelope:\n",
    "def Envelope(t):\n",
    "    if abs((t-tpeak)/tpeak) < 1:\n",
    "        return (2*exp(1)/sqrt(pi))*sin(2*pi*fcen*t)*exp (-1/(1-((t-tpeak)/tpeak)**2))\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "369555b1",
   "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 = 12\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: scene.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3af4e308",
   "metadata": {},
   "source": [
    "The differential operators and coupling terms to the pml - variables are represented via linear operators;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff4cb3ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "print (fullB.GetOperatorInfo())\n",
    "print (dampingp.GetOperatorInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94737b9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "print (op1.GetOperatorInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0d9314a-79eb-45fc-91b1-506a41c33277",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fdeb9fc-a372-43ff-b071-accb954e8dd1",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}