{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Example: Planewave Absorption Maximization\n",
    "In this example, we compute bounds on the maximum absorption of an incident planewave by a structure. \n",
    "\n",
    "The harmonic time-averaged absorption is\n",
    "\n",
    "$$\n",
    "f(\\mathbf{e}) = \\frac{\\omega}{2}  \\int \\Im(\\epsilon(\\mathbf{r})) |\\mathbf{e}(\\mathbf{r})|^2 d \\mathbf{r}\n",
    "$$\n",
    "\n",
    "Using the fact that $\\Im(\\epsilon(\\mathbf{r})) |\\mathbf{e}(\\mathbf{r})|^2$ is only non-zero over the structure, we can rewrite this using the polarization field as\n",
    "\n",
    "$$\n",
    "f(\\mathbf{p}) = \\frac{\\omega}{2} \\Im\\left(\\chi \\frac{\\mathbf{p}^\\dagger}{\\chi^*} \\frac{\\mathbf{p}}{\\chi} \\right) = \\frac{\\omega}{2}\\frac{\\Im(\\chi)}{|\\chi|^2} \\mathbf{p}^\\dagger \\mathbf{p}.\n",
    "$$\n",
    "Therefore, the constant and linear parts of the objective are $0$, and the quadratic part for the dense form is $A_0 = -\\frac{\\omega}{2}\\frac{\\Im(\\chi)}{|\\chi|^2} I = \\frac{\\omega}{2} \\Im\\left(\\frac{1}{\\chi} \\right) I$. Note the $-$ sign when defining $A_0$: this is because the QCQP class writes the QCQP objective as\n",
    "$$\n",
    "\\max_{\\mathbf{x}} \\quad -\\mathbf{x}^\\dagger A_0 \\mathbf{x} + 2\\Re(\\mathbf{x^\\dagger} s_0) + c_0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import sys, time, os\n",
    "from dolphindes import photonics, geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "## wavelength, geometry and materials of the planewave absorption problem ##\n",
    "wavelength = 1.0\n",
    "omega = 2*np.pi / wavelength\n",
    "\n",
    "chi = 3 + 1e-2j\n",
    "\n",
    "px_per_length = 20 # pixels per length unit. If wavelength = 1.0, then this is pixels per wavelength.\n",
    "dl = 1 / px_per_length\n",
    "\n",
    "des_x = 1.5\n",
    "des_y = 1.5 # size of the design region for the absorbing structure\n",
    "pmlsep = 1.0\n",
    "pmlthick = 0.5\n",
    "Mx = int(des_x / dl)\n",
    "My = int(des_y / dl)\n",
    "\n",
    "Npmlsepx = Npmlsepy = int(pmlsep / dl)\n",
    "Npmlx = Npmly = int(pmlthick / dl)\n",
    "Nx = Mx + 2*(Npmlsepx + Npmlx)\n",
    "Ny = My + 2*(Npmlsepy + Npmly)\n",
    "\n",
    "des_mask = np.zeros((Nx,Ny), dtype=bool)\n",
    "des_mask[Npmlx+Npmlsepx:-(Npmlx+Npmlsepx) , Npmly+Npmlsepy:-(Npmly+Npmlsepy)] = True\n",
    "Ndes = int(np.sum(des_mask))\n",
    "\n",
    "## planewave source\n",
    "ji = np.zeros((Nx,Ny), dtype=complex)\n",
    "ji[Npmlx,:] = 2.0 / dl # linesource for unit amplitude planewave traveeling in x direction\n",
    "\n",
    "plt.matshow(des_mask + np.real(ji)*dl) # visualize where the mask and the source are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "## setup geometry\n",
    "geo = geometry.CartesianFDFDGeometry(\n",
    "    Nx=Nx, Ny=Ny, Npmlx=Npmlx, Npmly=Npmly, dx=dl, dy=dl\n",
    ")\n",
    "\n",
    "## setup the photonic TM FDFD class\n",
    "abs_problem = photonics.Photonics_TM_FDFD(\n",
    "    omega=omega, geometry=geo, chi=chi,\n",
    "    des_mask=des_mask, ji=ji, sparseQCQP=True\n",
    ")\n",
    "\n",
    "print(abs_problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "## obtain incident planewave and plot it\n",
    "ei = abs_problem.get_ei(ji, update=True)\n",
    "plt.imshow(np.real(ei), cmap='bwr')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "## We start by calculating limits with just the global conservation of power constraints\n",
    "c0 = 0.0\n",
    "s0_p = np.zeros(Ndes, dtype=complex)\n",
    "A0_p = (omega/2) * np.imag(1.0/chi) * sp.eye_array(Ndes)\n",
    "abs_problem.set_objective(s0=s0_p, A0=A0_p, c0=c0, denseToSparse=True) # using sparse formulation\n",
    "abs_problem.setup_QCQP(Pdiags='global', verbose=1)\n",
    "result = abs_problem.bound_QCQP('newton')\n",
    "print(f'global constraint bound is {result[0]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "## if desired, we can run generalized constraint descent to tighten the bounds\n",
    "abs_problem.QCQP.run_gcd()\n",
    "abs_gcd_bound = abs_problem.QCQP.current_dual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Bound for maximum absorption power found is', abs_gcd_bound)\n",
    "Sx = 0.5 # unit amplitude planewave has Poynting vector magnitude of 0.5\n",
    "print('Ratio between absorption cross section bound and geometric cross section is', abs_gcd_bound / (Sx*des_y))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "## Compare with Inverse design\n",
    "\n",
    "We can perform topology optimization to compare with the bounds.\n",
    "For convenience, dolphindes provides built-in adjoint calculation once the design objective is specified, which can be combined with a non-linear optimization package (like [NLopt](https://nlopt.readthedocs.io/en/latest/) used here) for easy inverse design. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import nlopt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "ndof = int(np.sum(des_mask))\n",
    "dof = 0.5 * np.ones(ndof) # half slab initialization\n",
    "\n",
    "opt_data = {'count':0} # dictionary keeping track of optimization progress\n",
    "def abs_objective(dof, grad):\n",
    "    # wrapper function around the structure objective in abs_problem to provide printouts of optimization progress\n",
    "    obj = abs_problem.structure_objective(dof, grad)\n",
    "    print('at iteration #', opt_data['count'], 'the absorption found is', obj)\n",
    "    opt_data['count'] += 1\n",
    "    return obj\n",
    "    \n",
    "# nlopt convergence criteria\n",
    "maxeval = 1000\n",
    "ftol_rel = 1e-4\n",
    "\n",
    "opt = nlopt.opt(nlopt.LD_MMA, ndof)\n",
    "\n",
    "lb = np.zeros(ndof)\n",
    "ub = np.ones(ndof)\n",
    "opt.set_lower_bounds(lb)\n",
    "opt.set_upper_bounds(ub)\n",
    "opt.set_ftol_rel(ftol_rel)\n",
    "opt.set_maxeval(maxeval)\n",
    "\n",
    "opt.set_max_objective(abs_objective)\n",
    "opt_dof = opt.optimize(dof)\n",
    "opt_abs = opt.last_optimum_value()\n",
    "\n",
    "print('TopOpt found structure with highest absorption of', opt_abs)\n",
    "print('Ratio between absorption and geometric cross sections is', opt_abs / (Sx*des_y))\n",
    "print('Compare with the bound ratio:', abs_gcd_bound / (Sx*des_y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dolphindes",
   "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}