{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1.7.1 Perfectly Matched Layer (PML)\n",
    "====\n",
    "\n",
    "\n",
    "The use of (Perfectly Matched Layer) PML is a standard technique to numerically solve for outgoing waves in unbounded domains. Although scattering problems are posed in unbounded domains,  by bounding the scatterer and any inhomogeneities within a PML, one is able to truncate to a bounded computational domain. In this tutorial we see how to use PML for \n",
    "- source problems, and \n",
    "- eigenvalue problems (resonances)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Source problem\n",
    "\n",
    "Consider the problem of finding $u$ satisfying \n",
    "\n",
    "$$\n",
    "-\\Delta u - \\omega^2 u = f \\qquad \\text{ in } \\mathbb{R}^2\n",
    "$$\n",
    "\n",
    "together with the Sommerfeld (outgoing) radiation condition at infinity\n",
    "\n",
    "$$\n",
    "\\lim_{r \\to \\infty} r^{1/2}\n",
    "\\bigg( \n",
    "\\frac{\\partial u }{ \\partial r} - i \\omega u \n",
    "\\bigg) = 0\n",
    "$$\n",
    "\n",
    "where $r$ is the radial coordinate.  Below, we set an $f$ that represents a time-harmonic pulse that is  almost zero except for a small region."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a geometry which is divided into an `inner` disk and and an outer \n",
    "annulus called `pmlregion`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outer = Circle((0,0), 1.4).Face()\n",
    "outer.edges.name = 'outerbnd'\n",
    "inner = Circle((0,0), 1).Face()\n",
    "inner.edges.name = 'innerbnd'\n",
    "inner.faces.name ='inner'\n",
    "pmlregion = outer - inner\n",
    "pmlregion.faces.name = 'pmlregion'\n",
    "geo = OCCGeometry(Glue([inner, pmlregion]), dim=2)\n",
    "\n",
    "mesh = Mesh(geo.GenerateMesh (maxh=0.1))\n",
    "mesh.Curve(3)\n",
    "\n",
    "f = exp(-20**2*((x-0.3)*(x-0.3)+y*y))\n",
    "Draw(f, mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The PML facility in NGSolve implements a complex coordinate transformation on a given mesh region (which in  this example is  `pmlregion`). When this complex variable change is applied to the outgoing solution in the PML region,  it becomes a a function that decays exponentially as $r \\to \\infty$, allowing us to truncate the unbounded domain.\n",
    "\n",
    "With the following single line, we tell NGSolve to turn on this coordinate transformation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.SetPML(pml.Radial(rad=1,alpha=1j,origin=(0,0)), \"pmlregion\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then a  radial PML is set in the exterior of a disk \n",
    "centered at `origin`  of radius `rad`. In addition to \n",
    "`origin` and `rad`, the keyword argument `alpha` may be used to set\n",
    "the PML strength, which determines the rate of increase in the imaginary\n",
    "part of the coordinate map as radius increases.\n",
    "\n",
    "Having set the PML, the rest of the code now  looks very much like that  in [Unit 1.7](../unit-1.7-helmholtz/helmholtz.ipynb):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=4, complex=True)\n",
    "u = fes.TrialFunction()\n",
    "v = fes.TestFunction()\n",
    "\n",
    "omega = 10\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += grad(u)*grad(v)*dx - omega**2*u*v*dx\n",
    "a += -1j*omega*u*v*ds(\"outerbnd\")\n",
    "a.Assemble()\n",
    "\n",
    "b = LinearForm(f * v * dx).Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = a.mat.Inverse() * b.vec\n",
    "Draw(gfu);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that above we have kept the same Robin boundary conditions as in [Unit 1.7](../unit-1.7-helmholtz/helmholtz.ipynb).\n",
    "(Since the solution exponentially decays within PML, it is also common practice to put zero Dirichlet boundary conditions at the outer edge of the PML, an option that is also easily implementable in NGSolve using the `dirichlet` flag.)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigenvalue problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PML can also be used to compute resonances on unbounded domains. Resonances of a scatterer $S$ are outgoing waves $u$ satisfying the homogeneous Helmholtz equation\n",
    "\n",
    "$$\n",
    "\\Delta u + \\omega^2 u = 0\n",
    " \\qquad \\text{ in } \\mathbb{R}^2 \\setminus S,\n",
    "$$\n",
    "\n",
    "i.e., $u$ is an eigenfunction of $-\\Delta$ associated to eigenvalue $\\lambda :=\\omega^2.$\n",
    "\n",
    "We truncate this unbounded-domain eigenproblem to a bounded-domain eigenproblem using PML. The following example computes the resonances of a cavity opening on one side into free space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shroom = MoveTo(-2, 0).Line(1.8).Rotate(-90).Line(0.8).Rotate(90).Line(0.4). \\\n",
    "    Rotate(90).Line(0.8).Rotate(-90).Line(1.8).Rotate(90).Arc(2, 180).Face()\n",
    "circ = Circle((0, 0), 1).Face()\n",
    "pmlregion = shroom - circ\n",
    "pmlregion.faces.name = 'pml'\n",
    "air = shroom * circ\n",
    "air.faces.name = 'air'\n",
    "shape = Glue([air, pmlregion])\n",
    "geo = OCCGeometry(shape, dim=2)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.1))\n",
    "mesh.Curve(5)\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The top edge of the rectangular cavity opens into air, which is included in the computational domain as a semicircular region (neglecting backward propagating waves), while the remaining edges of the cavity are isolated from rest of space by perfect reflecting boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.SetPML(pml.Radial(rad=1, alpha=1j, origin=(0,0)), \"pml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now make the complex symmetric system with PML incorporated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=4, complex=True, dirichlet=\"dir\")\n",
    "\n",
    "u = fes.TrialFunction()\n",
    "v = fes.TestFunction()\n",
    "\n",
    "a = BilinearForm(grad(u)*grad(v)*dx).Assemble()\n",
    "m = BilinearForm(u*v*dx).Assemble();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For solving the eigenproblem, we use an Arnoldi eigensolver, to which a `GridFunction` with many vectors (allocated via keyword argument `multidim`) is input, together with `a` and `b`. A `shift` argument to `ArnoldiSolver` indicates the approximate location of expected eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = GridFunction(fes, multidim=50, name='resonances')\n",
    "\n",
    "with TaskManager():\n",
    "    lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(), \n",
    "                        list(u.vecs), shift=400)\n",
    "Draw(u);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (\"lam: \", lam)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The more \"confined\" modes have resonance values $\\omega$ that are closer to the real axis. Here are the computed $\\omega$-values plotted in the complex plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lamr = [sqrt(l).real for l in lam]\n",
    "lami = [sqrt(l).imag for l in lam]\n",
    "plt.plot(lamr, lami, \".\")\n",
    "plt.grid(True)\n",
    "plt.title('Computed $\\omega$ values')\n",
    "plt.show()"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}