python/examples/perturbation_theory.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sensitivity Analysis via Perturbation Theory\n",
    "--------------------------------------------\n",
    "\n",
    "For a given mode of the ring resonator, it is often useful to know how sensitive the resonant frequency $\\omega$ is to small changes in the ring radius $r$ by computing its derivative $\\partial\\omega/\\partial r$. The gradient is also a useful quantity for shape optimization because it can be paired with fast iterative methods such as [BFGS](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm) to find local optima. The \"brute-force\" approach for computing the gradient is via a finite-difference approximation requiring *two* simulations of the (1) unperturbed [$\\omega(r)$] and (2) perturbed [$\\omega(r+\\Delta r)$] structures. Since each simulation is potentially costly, an alternative approach based on semi analytics is to use [perturbation theory](https://en.wikipedia.org/wiki/Perturbation_theory) to obtain the gradient from the fields of the unperturbed structure. This involves a single simulation and is often more computationally efficient than the brute-force approach although some care is required to set up the calculation properly.  (More generally, [adjoint methods](https://math.mit.edu/~stevenj/18.336/adjoint.pdf) can be used to compute any number of derivatives with a single additional simulation.)\n",
    "\n",
    "[Pertubation theory for Maxwell equations involving high index-contrast dielectric interfaces](http://math.mit.edu/~stevenj/papers/KottkeFa08.pdf) is reviewed in Chapter 2 of [Photonics Crystals: Molding the Flow of Light, 2nd Edition (2008)](http://ab-initio.mit.edu/book/). The formula (equation 30 on p.19) for the frequency shift $\\Delta \\omega$ resulting from the displacement of a block of $\\varepsilon_1$-material towards $\\varepsilon_2$-material by a distance $\\Delta h$ (perpendicular to the boundary) is:\n",
    "\n",
    "<center>\n",
    "\n",
    "$$ \\Delta\\omega = -\\frac{\\omega}{2} \\frac{ \\iint d^2 \\vec{r} \\big[ (\\varepsilon_1 - \\varepsilon_2) |\\vec{E}_{\\parallel}(\\vec{r})|^2 - \\big(\\frac{1}{\\varepsilon_1} - \\frac{1}{\\varepsilon_2}\\big)|\\varepsilon\\vec{E}_{\\perp}|^2\\big] \\Delta h}{\\int d^3\\vec{r} \\varepsilon(\\vec{r})|\\vec{E}(\\vec{r})|^2} + O(\\Delta h^2)$$\n",
    "\n",
    "</center>\n",
    "\n",
    "In this expression, $\\vec{E}_{\\parallel}(\\vec{r})$ is the component of $\\vec{E}$ that is parallel to the surface, and $\\varepsilon\\vec{E}_{\\perp}$ is the component of $\\varepsilon\\vec{E}$ that is perpendicular to the surface. These two components are guaranteed to be continuous across an interface between two isotropic dielectric materials. In this demonstration, $\\partial\\omega/\\partial r$ is computed using this formula and the results are validated by comparing with the finite-difference approximation: $[\\omega(r+\\Delta r)-\\omega(r)]/\\Delta r$.\n",
    "\n",
    "There are three parts to the calculation: (1) find the resonant frequency of the unperturbed geometry using a broadband Gaussian source, (2) find the resonant mode profile of the unperturbed geometry using a narrowband source and from these fields compute the gradient via the perturbation-theory formula, and (3) find the resonant frequency of the perturbed geometry and from this compute the gradient using a finite-difference approximation. The perturbation is applied only to the inner and outer ring radii. The ring width is constant. There are two types of modes which are computed in separate simulations using different source polarizations: parallel ($E_z$) and perpendicular ($H_z$) to the interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------\n",
      "Initializing structure...\n",
      "time for choose_chunkdivision = 0.000187874 s\n",
      "Working in Cylindrical dimensions.\n",
      "Computational cell is 8 x 0 x 0.01 with resolution 100\n",
      "     block, center = (1.5,0,0)\n",
      "          size (1,1e+20,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (11.56,11.56,11.56)\n",
      "time for set_epsilon = 0.00455093 s\n",
      "-----------\n",
      "Meep: using complex fields.\n",
      "Meep progress: 63.24/394.1176452636719 = 16.0% done in 4.0s, 20.9s to go\n",
      "on time step 12648 (time=63.24), 0.000316286 s/step\n",
      "Meep progress: 127.035/394.1176452636719 = 32.2% done in 8.0s, 16.8s to go\n",
      "on time step 25410 (time=127.05), 0.00031345 s/step\n",
      "Meep progress: 190.945/394.1176452636719 = 48.4% done in 12.0s, 12.8s to go\n",
      "on time step 38195 (time=190.975), 0.00031288 s/step\n",
      "Meep progress: 254.66/394.1176452636719 = 64.6% done in 16.0s, 8.8s to go\n",
      "on time step 50941 (time=254.705), 0.000313832 s/step\n",
      "Meep progress: 316.92/394.1176452636719 = 80.4% done in 20.0s, 4.9s to go\n",
      "on time step 63396 (time=316.98), 0.000321182 s/step\n",
      "Meep progress: 376.46/394.1176452636719 = 95.5% done in 24.0s, 1.1s to go\n",
      "on time step 75307 (time=376.535), 0.000335853 s/step\n",
      "harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error\n",
      "harminv0:, 0.17578499227611236, -4.816700722801396e-05, 1824.7448034707427, 0.24613200861807683, -0.14560581011759283-0.19844372937023874i, 1.1797259436409042e-10+0.0i\n",
      "run 0 finished at t = 394.12 (78824 timesteps)\n",
      "-----------\n",
      "Initializing structure...\n",
      "time for choose_chunkdivision = 0.000282049 s\n",
      "Working in Cylindrical dimensions.\n",
      "Computational cell is 8 x 0 x 0.01 with resolution 100\n",
      "     block, center = (1.5,0,0)\n",
      "          size (1,1e+20,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (11.56,11.56,11.56)\n",
      "time for set_epsilon = 0.00452399 s\n",
      "-----------\n",
      "Meep: using complex fields.\n",
      "Meep progress: 64.11/1237.7535400390625 = 5.2% done in 4.0s, 73.2s to go\n",
      "on time step 12822 (time=64.11), 0.00031198 s/step\n",
      "Meep progress: 129.395/1237.7535400390625 = 10.5% done in 8.0s, 68.5s to go\n",
      "on time step 25882 (time=129.41), 0.0003063 s/step\n",
      "Meep progress: 194.32/1237.7535400390625 = 15.7% done in 12.0s, 64.4s to go\n",
      "on time step 38871 (time=194.355), 0.000307977 s/step\n",
      "Meep progress: 260.045/1237.7535400390625 = 21.0% done in 16.0s, 60.2s to go\n",
      "on time step 52019 (time=260.095), 0.00030423 s/step\n",
      "Meep progress: 325.825/1237.7535400390625 = 26.3% done in 20.0s, 56.0s to go\n",
      "on time step 65178 (time=325.89), 0.000303988 s/step\n",
      "Meep progress: 390.665/1237.7535400390625 = 31.6% done in 24.0s, 52.0s to go\n",
      "on time step 78150 (time=390.75), 0.00030838 s/step\n",
      "Meep progress: 455.2/1237.7535400390625 = 36.8% done in 28.0s, 48.1s to go\n",
      "on time step 91060 (time=455.3), 0.000309849 s/step\n",
      "Meep progress: 519.66/1237.7535400390625 = 42.0% done in 32.0s, 44.2s to go\n",
      "on time step 103958 (time=519.79), 0.000310148 s/step\n",
      "Meep progress: 584.73/1237.7535400390625 = 47.2% done in 36.0s, 40.2s to go\n",
      "on time step 116975 (time=584.875), 0.000307307 s/step\n",
      "Meep progress: 649.64/1237.7535400390625 = 52.5% done in 40.0s, 36.2s to go\n",
      "on time step 129957 (time=649.785), 0.000308119 s/step\n",
      "Meep progress: 714.945/1237.7535400390625 = 57.8% done in 44.0s, 32.2s to go\n",
      "on time step 143025 (time=715.125), 0.000306096 s/step\n",
      "Meep progress: 780.855/1237.7535400390625 = 63.1% done in 48.0s, 28.1s to go\n",
      "on time step 156209 (time=781.045), 0.000303408 s/step\n",
      "Meep progress: 846.395/1237.7535400390625 = 68.4% done in 52.0s, 24.0s to go\n",
      "on time step 169322 (time=846.61), 0.000305051 s/step\n",
      "Meep progress: 911.86/1237.7535400390625 = 73.7% done in 56.0s, 20.0s to go\n",
      "on time step 182419 (time=912.095), 0.000305427 s/step\n",
      "Meep progress: 977.405/1237.7535400390625 = 79.0% done in 60.0s, 16.0s to go\n",
      "on time step 195532 (time=977.66), 0.000305057 s/step\n",
      "Meep progress: 1042.795/1237.7535400390625 = 84.2% done in 64.0s, 12.0s to go\n",
      "on time step 208616 (time=1043.08), 0.000305738 s/step\n",
      "Meep progress: 1107.3500000000001/1237.7535400390625 = 89.5% done in 68.0s, 8.0s to go\n",
      "on time step 221526 (time=1107.63), 0.000309854 s/step\n",
      "Meep progress: 1172.115/1237.7535400390625 = 94.7% done in 72.0s, 4.0s to go\n",
      "on time step 234482 (time=1172.41), 0.000308753 s/step\n",
      "Meep progress: 1237.6100000000001/1237.7535400390625 = 100.0% done in 76.0s, 0.0s to go\n",
      "run 0 finished at t = 1237.755 (247551 timesteps)\n",
      "-----------\n",
      "Initializing structure...\n",
      "time for choose_chunkdivision = 0.000117064 s\n",
      "Working in Cylindrical dimensions.\n",
      "Computational cell is 8 x 0 x 0.01 with resolution 100\n",
      "     block, center = (1.51,0,0)\n",
      "          size (1,1e+20,1e+20)\n",
      "          axes (1,0,0), (0,1,0), (0,0,1)\n",
      "          dielectric constant epsilon diagonal = (11.56,11.56,11.56)\n",
      "time for set_epsilon = 0.0026741 s\n",
      "-----------\n",
      "Meep: using complex fields.\n",
      "Meep progress: 64.245/1237.7535400390625 = 5.2% done in 4.0s, 73.1s to go\n",
      "on time step 12849 (time=64.245), 0.000311334 s/step\n",
      "Meep progress: 128.08/1237.7535400390625 = 10.3% done in 8.0s, 69.3s to go\n",
      "on time step 25621 (time=128.105), 0.000313206 s/step\n",
      "Meep progress: 192.4/1237.7535400390625 = 15.5% done in 12.0s, 65.2s to go\n",
      "on time step 38487 (time=192.435), 0.0003109 s/step\n",
      "Meep progress: 256.76/1237.7535400390625 = 20.7% done in 16.0s, 61.1s to go\n",
      "on time step 51365 (time=256.825), 0.00031061 s/step\n",
      "Meep progress: 321.26/1237.7535400390625 = 26.0% done in 20.0s, 57.1s to go\n",
      "on time step 64271 (time=321.355), 0.000309956 s/step\n",
      "Meep progress: 385.66/1237.7535400390625 = 31.2% done in 24.0s, 53.0s to go\n",
      "on time step 77154 (time=385.77), 0.000310499 s/step\n",
      "Meep progress: 449.995/1237.7535400390625 = 36.4% done in 28.0s, 49.0s to go\n",
      "on time step 90026 (time=450.13), 0.000310771 s/step\n",
      "Meep progress: 514.22/1237.7535400390625 = 41.5% done in 32.0s, 45.0s to go\n",
      "on time step 102876 (time=514.38), 0.000311289 s/step\n",
      "Meep progress: 578.95/1237.7535400390625 = 46.8% done in 36.0s, 41.0s to go\n",
      "on time step 115826 (time=579.13), 0.000308899 s/step\n",
      "Meep progress: 642.9250000000001/1237.7535400390625 = 51.9% done in 40.0s, 37.0s to go\n",
      "on time step 128622 (time=643.11), 0.000312613 s/step\n",
      "Meep progress: 706.915/1237.7535400390625 = 57.1% done in 44.0s, 33.0s to go\n",
      "on time step 141424 (time=707.12), 0.00031246 s/step\n",
      "Meep progress: 770.58/1237.7535400390625 = 62.3% done in 48.0s, 29.1s to go\n",
      "on time step 154164 (time=770.82), 0.000313974 s/step\n",
      "Meep progress: 834.445/1237.7535400390625 = 67.4% done in 52.0s, 25.1s to go\n",
      "on time step 166942 (time=834.71), 0.00031304 s/step\n",
      "Meep progress: 898.195/1237.7535400390625 = 72.6% done in 56.0s, 21.2s to go\n",
      "on time step 179696 (time=898.48), 0.000313635 s/step\n",
      "Meep progress: 961.405/1237.7535400390625 = 77.7% done in 60.0s, 17.2s to go\n",
      "on time step 192338 (time=961.69), 0.000316413 s/step\n",
      "Meep progress: 1025.53/1237.7535400390625 = 82.9% done in 64.0s, 13.2s to go\n",
      "on time step 205173 (time=1025.87), 0.000311651 s/step\n",
      "Meep progress: 1089.78/1237.7535400390625 = 88.0% done in 68.0s, 9.2s to go\n",
      "on time step 218020 (time=1090.1), 0.000311357 s/step\n",
      "Meep progress: 1151.805/1237.7535400390625 = 93.1% done in 72.0s, 5.4s to go\n",
      "on time step 230424 (time=1152.12), 0.000322493 s/step\n",
      "Meep progress: 1210.695/1237.7535400390625 = 97.8% done in 76.0s, 1.7s to go\n",
      "on time step 242207 (time=1211.04), 0.000339497 s/step\n",
      "harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error\n",
      "harminv0:, 0.17493286736703903, -4.814897914072506e-05, 1816.5791932551942, 1.216470973025523, 1.2069307717323075-0.15205176901081904i, 2.316660865172475e-11+0.0i\n",
      "run 0 finished at t = 1237.755 (247551 timesteps)\n",
      "dwdR:, -0.08544696397218933 (pert. theory), -0.08521249090733263 (finite diff.)\n"
     ]
    }
   ],
   "source": [
    "import meep as mp\n",
    "import numpy as np\n",
    "\n",
    "resolution = 100        # pixels/um\n",
    "\n",
    "perpendicular = False   # perpendicular (Hz) or parallel (Ez) source?\n",
    "\n",
    "if perpendicular:\n",
    "    src_cmpt = mp.Hz\n",
    "    fcen = 0.21         # pulse center frequency\n",
    "else:\n",
    "    src_cmpt = mp.Ez\n",
    "    fcen = 0.17         # pulse center frequency\n",
    "\n",
    "n = 3.4                 # index of waveguide\n",
    "w = 1                   # ring width\n",
    "r = 1                   # inner radius of ring\n",
    "pad = 4                 # padding between waveguide and edge of PML\n",
    "dpml = 2                # thickness of PML\n",
    "m = 5                   # angular dependence\n",
    "\n",
    "pml_layers = [mp.PML(dpml)]\n",
    "\n",
    "sr = r + w + pad + dpml        # radial size (cell is from 0 to sr)\n",
    "dimensions = mp.CYLINDRICAL    # coordinate system is (r,phi,z) instead of (x,y,z)\n",
    "cell = mp.Vector3(sr)\n",
    "\n",
    "geometry = [mp.Block(center=mp.Vector3(r + (w / 2)),\n",
    "                     size=mp.Vector3(w, mp.inf, mp.inf),\n",
    "                     material=mp.Medium(index=n))]\n",
    "\n",
    "# find resonant frequency of unperturbed geometry using broadband source\n",
    "\n",
    "df = 0.2*fcen           # pulse width (in frequency)\n",
    "\n",
    "sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),\n",
    "                     component=src_cmpt,\n",
    "                     center=mp.Vector3(r+0.1))]\n",
    "\n",
    "sim = mp.Simulation(cell_size=cell,\n",
    "                    geometry=geometry,\n",
    "                    boundary_layers=pml_layers,\n",
    "                    resolution=resolution,\n",
    "                    sources=sources,\n",
    "                    dimensions=dimensions,\n",
    "                    m=m)\n",
    "\n",
    "h = mp.Harminv(src_cmpt, mp.Vector3(r+0.1), fcen, df)\n",
    "sim.run(mp.after_sources(h),\n",
    "        until_after_sources=100)\n",
    "\n",
    "frq_unperturbed = h.modes[0].freq\n",
    "\n",
    "sim.reset_meep()\n",
    "\n",
    "# unperturbed geometry with narrowband source centered at resonant frequency\n",
    "\n",
    "fcen = frq_unperturbed\n",
    "df = 0.05*fcen\n",
    "\n",
    "sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),\n",
    "                     component=src_cmpt,\n",
    "                     center=mp.Vector3(r+0.1))]\n",
    "\n",
    "sim = mp.Simulation(cell_size=cell,\n",
    "                    geometry=geometry,\n",
    "                    boundary_layers=pml_layers,\n",
    "                    resolution=resolution,\n",
    "                    sources=sources,\n",
    "                    dimensions=dimensions,\n",
    "                    m=m)\n",
    "\n",
    "sim.run(until_after_sources=100)\n",
    "\n",
    "deps = 1 - n**2\n",
    "deps_inv = 1 - 1/n**2\n",
    "\n",
    "if perpendicular:\n",
    "    para_integral = deps*2*np.pi*(r*abs(sim.get_field_point(mp.Ep, mp.Vector3(r)))**2 - (r+w)*abs(sim.get_field_point(mp.Ep, mp.Vector3(r+w)))**2)\n",
    "    perp_integral = deps_inv*2*np.pi*(-r*abs(sim.get_field_point(mp.Dr, mp.Vector3(r)))**2 + (r+w)*abs(sim.get_field_point(mp.Dr, mp.Vector3(r+w)))**2)\n",
    "    numerator_integral = para_integral + perp_integral\n",
    "else:\n",
    "    numerator_integral = deps*2*np.pi*(r*abs(sim.get_field_point(mp.Ez, mp.Vector3(r)))**2 - (r+w)*abs(sim.get_field_point(mp.Ez, mp.Vector3(r+w)))**2)\n",
    "\n",
    "denominator_integral = sim.electric_energy_in_box(center=mp.Vector3(0.5*sr), size=mp.Vector3(sr))\n",
    "\n",
    "perturb_theory_dw_dR = -frq_unperturbed * numerator_integral / (4 * denominator_integral)\n",
    "\n",
    "sim.reset_meep()\n",
    "\n",
    "# perturbed geometry with narrowband source\n",
    "\n",
    "dr = 0.01\n",
    "\n",
    "sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),\n",
    "                     component=src_cmpt,\n",
    "                     center=mp.Vector3(r + dr + 0.1))]\n",
    "\n",
    "geometry = [mp.Block(center=mp.Vector3(r + dr + (w / 2)),\n",
    "                     size=mp.Vector3(w, mp.inf, mp.inf),\n",
    "                     material=mp.Medium(index=n))]\n",
    "\n",
    "sim = mp.Simulation(cell_size=cell,\n",
    "                    geometry=geometry,\n",
    "                    boundary_layers=pml_layers,\n",
    "                    resolution=resolution,\n",
    "                    sources=sources,\n",
    "                    dimensions=dimensions,\n",
    "                    m=m)\n",
    "\n",
    "h = mp.Harminv(src_cmpt, mp.Vector3(r+dr+0.1), fcen, df)\n",
    "sim.run(mp.after_sources(h),\n",
    "        until_after_sources=100)\n",
    "\n",
    "frq_perturbed = h.modes[0].freq\n",
    "\n",
    "finite_diff_dw_dR = (frq_perturbed - frq_unperturbed) / dr\n",
    "\n",
    "print(\"dwdR:, {} (pert. theory), {} (finite diff.)\".format(perturb_theory_dw_dR,finite_diff_dw_dR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are three things to note. First, there is a built-in function `electric_energy_in_box` which calculates the integral of $\\vec{E}\\cdot\\vec{D}/2 = \\varepsilon|E|^2/2$. This is exactly the integral in the denominator, divided by 2. In cylindrical coordinates $(r,\\phi,z)$, the integrand is [multiplied](https://en.wikipedia.org/wiki/Cylindrical_coordinate_system#Line_and_volume_elements) by the circumference $2\\pi r$, or equivalently the integral is over an annular volume. Second, for the case involving the $H_z$ source, both parallel ($E_{\\parallel}=E_{\\phi}$) and perpendicular ($\\varepsilon E_{\\perp}=D_r$) fields are present which must be included in the numerator as separate terms. Field values anywhere in the grid obtained with `get_field_point` are [automatically interpolated](https://meep.readthedocs.io/en/latest/Introduction/#the-illusion-of-continuity); i.e., no additional post-processing is necessary. Third, when comparing the perturbation-theory result to the finite-difference approximation, there are *two* convergence parameters: the resolution and $\\Delta r$. In principle, to demonstrate agreement with perturbation theory, the limit of the resolution should be taken to infinity and *then* the limit of $\\Delta r$ to 0. In practice, this can be obtained by doubling the resolution at a given $\\Delta r$ until it is sufficiently converged, then halving $\\Delta r$ and repeating.\n",
    "\n",
    "For an $E_z$ source (parallel to the interface), at a resolution of 100 the results are -0.0854469639 (perturbation theory) and -0.0852124909 (finite difference). Doubling the resolution to 200, the results are -0.0854460732 (perturbation theory) and -0.0852115350 (finite difference). Both results have converged to at least five digits. The relative error at resolution 200 is 0.3%. The mode has a $\\omega$ of 0.175 and $Q$ of 1800.\n",
    "\n",
    "For an $H_z$ source (perpendicular to the interface), at a resolution of 100 the results are -0.0805038572 (perturbation theory) and -0.0798087331 (finite difference). Doubling the resolution to 200, the results are -0.0802028346 (perturbation theory) and -0.0798088015 (finite difference). Both results have converged to at least three digits. The relative error at resolution 200 is 0.5%. The error is larger in this case due to the presence of the discontinuous fields at the dielectric interface which are handled by [subpixel smoothing](https://meep.readthedocs.io/en/latest/Subpixel_Smoothing/). The mode has a $\\omega$ of 0.208 and $Q$ of 1200.\n",
    "\n",
    "Finally, as reference, the same calculation can be set up in Cartesian coordinates as a 2d simulation. There is one major difference: the mode produced by a point source in 2d is actually the $\\cos(m\\phi)$ mode, *not* $\\exp(im\\phi)$, or equivalently it is the superposition of the $\\pm m$ modes. This means that computing the numerator integral does not involve just multiplying the field at a single point on the surface by $2\\pi r$ &mdash; rather, it is the integral of $\\cos^2(m\\phi)$ which gives a factor of 1/2. (For non-circular shapes in 2d, the surface integral must be computed numerically.) The simulation script is in [examples/perturbation_theory_2d.py](https://github.com/NanoComp/meep/blob/master/python/examples/perturbation_theory_2d.py). The results are comparable to the cylindrical coordinate case (a 1d calculation) but the 2d simulation is much slower and less accurate at the same grid resolution."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}