NGSolve24

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(complementary:complex)=\n",
    "# Complex-valued waves\n",
    "\n",
    "\n",
    "In NGSolve finite element spaces can be built over the complex field and the resulting  complex matrix systems can be solved. This tutorial shows how to compute the solution of the Helmholtz equation with impedance boundary conditions in complex arithmetic. The boundary value problem is to find $u$ satisfying \n",
    "\n",
    "$$\n",
    "-\\Delta u - \\omega^2 u = f\\qquad \\text{ in } \\Omega\n",
    "$$\n",
    "\n",
    "together with the impedance (outgoing) boundary condition\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u }{ \\partial n} - i \\omega u = 0 \n",
    "\\quad \\text{ on } \\partial \\Omega\n",
    "$$\n",
    "\n",
    "where $i =$ `1j` is the imaginary unit.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.occ import *\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "air: Face = Circle((0.5, 0.5), 0.8).Face()\n",
    "air.edges.name = 'outer'\n",
    "scatterer = MoveTo(0.7, 0.3).Rectangle(0.05, 0.4).Face()\n",
    "scatterer.edges.name = 'scat'\n",
    "geo = OCCGeometry(air - scatterer, dim=2)\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))\n",
    "Draw(mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Declare  a complex-valued finite element space "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=5, complex=True)\n",
    "u, v = fes.TnT()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = 10\n",
    "pulse = 5e4*exp(-(40**2)*((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5)))\n",
    "Draw(10**-5*pulse, mesh, order=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forming the system\n",
    "\n",
    "The weak form for $u \\in H^1$:\n",
    "\n",
    "$$\n",
    "\\int_\\Omega\\big[ \\nabla u \\cdot \\nabla \\bar v - \\omega^2 u \\bar v \\big]\n",
    "\\, dx - \n",
    "i \\,\\omega\\, \\int_{\\partial \\Omega} u \\bar v \\, ds = \\int_{\\Omega} f \\bar v\n",
    "$$\n",
    "\n",
    "for all $v$ in $H^1$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Forms\n",
    "\n",
    "a = BilinearForm(fes)\n",
    "a += grad(u)*grad(v)*dx - omega**2*u*v*dx\n",
    "a += -omega*1j*u*v * ds(\"outer\")\n",
    "a.Assemble()\n",
    "\n",
    "\n",
    "f = LinearForm(pulse * v * dx).Assemble();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfu = GridFunction(fes, name=\"u\")\n",
    "gfu.vec.data = a.mat.Inverse() * f.vec\n",
    "Draw(gfu, mesh, min=-10, max=10, order=3, animate_complex=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open controls and explore webgui's further viewing options such as viewing real and imaginary parts, viewing absolute value, and viewing with deformation turned on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add a PML \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the same problem as before\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",
    "define now an outer region called for the PML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "outer = Circle((0.5, 0.5), 1).Face()\n",
    "outer.edges.name = 'outer'\n",
    "air: Face = Circle((0.5, 0.5), 0.8).Face()\n",
    "#air.edges.name = 'outer'\n",
    "scatterer = MoveTo(0.7, 0.3).Rectangle(0.05, 0.4).Face()\n",
    "scatterer.edges.name = 'scat'\n",
    "\n",
    "pmlregion = outer - air\n",
    "pmlregion.faces.name = 'pmlregion'\n",
    "geo = OCCGeometry(Glue([air-scatterer, pmlregion]), dim=2)\n",
    "\n",
    "mesh = Mesh(geo.GenerateMesh(maxh=0.05))\n",
    "Draw(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=0.8, alpha=1j,origin=(0.5,0.5)), \"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",
    "\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",
    "f = LinearForm(pulse * v * dx).Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "gfu.vec.data = a.mat.Inverse() * f.vec\n",
    "\n",
    "settings = {}\n",
    "Draw(gfu, animate_complex=True, deformation=False, scale = 0.01);"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}