{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Polar Solver Inner PML test\n",
    "Checking the veracity and usefullness of having inner PML in the polar coordinate Helmholtz solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "import scipy.sparse.linalg as spla\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "import sys\n",
    "sys.path.append(\"../../\")\n",
    "\n",
    "from dolphindes.geometry import PolarFDFDGeometry\n",
    "from dolphindes.maxwell import TM_Polar_FDFD, plot_real_polar_field, plot_cplx_polar_field, expand_symmetric_field"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2",
   "metadata": {},
   "source": [
    "## Computational domain starting at origin\n",
    "\n",
    "The following results suggest that the current inner pml implementation can be problematic when the computational domain starts at the origin.\n",
    "\n",
    "### setting up geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# domain parameters\n",
    "wvlgth0 = 1.0 # center wavelength of interest\n",
    "omega0 = 2*np.pi / wvlgth0 # center frequency of interest\n",
    "\n",
    "r_tot = 4.0 #  computational domain outer radius (includes all pml)\n",
    "w_pml_outer = 0.5 # thickness of outer layer of pml\n",
    "gpr = 80\n",
    "dr = 1.0 / gpr\n",
    "Nr = int(np.round(r_tot / dr))\n",
    "Npml = int(np.round(w_pml_outer/dr))\n",
    "\n",
    "Nphi_tot = 720 # total azimuthal gridlines around full circle\n",
    "n_sectors = 6\n",
    "Nphi_sector = Nphi_tot // n_sectors\n",
    "\n",
    "# no inner PML solver, for reference\n",
    "geo_ref = PolarFDFDGeometry(Nphi_sector, Nr, Npml, dr, n_sectors=n_sectors)\n",
    "FDFD_ref = TM_Polar_FDFD(omega0, geo_ref)\n",
    "\n",
    "# thin inner PML\n",
    "w_pml_inner_thin = 0.5\n",
    "Npml_inner_thin = int(np.round(w_pml_inner_thin / dr))\n",
    "geo_IPML_thin = PolarFDFDGeometry(Nphi_sector, Nr, Npml, dr, \n",
    "                                  n_sectors=n_sectors,\n",
    "                                  Npml_inner=Npml_inner_thin)\n",
    "FDFD_IPML_thin = TM_Polar_FDFD(omega0, geo_IPML_thin)\n",
    "\n",
    "# thick inner PML\n",
    "w_pml_inner_thick = 1.0\n",
    "Npml_inner_thick = int(np.round(w_pml_inner_thick / dr))\n",
    "geo_IPML_thick = PolarFDFDGeometry(Nphi_sector, Nr, Npml, dr, \n",
    "                                  n_sectors=n_sectors,\n",
    "                                  Npml_inner=Npml_inner_thick)\n",
    "FDFD_IPML_thick = TM_Polar_FDFD(omega0, geo_IPML_thick)\n",
    "\n",
    "# get coordinate grids\n",
    "phi_grid_sector, r_grid, phi_grid_full = FDFD_ref.get_symmetric_grids()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "### uniform ring source\n",
    "excites only m=0 waves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# current source: a uniform delta(r - r0)\n",
    "r0 = r_tot / 2\n",
    "J_rgrid = np.zeros(Nr, dtype=complex)\n",
    "J_rgrid[int(r0 / dr)] = 1.0 / (2*np.pi*r0*dr)\n",
    "J_grid = np.kron(np.ones(Nphi_sector), J_rgrid)\n",
    "\n",
    "print('current ring')\n",
    "plot_real_polar_field(np.real(J_grid), phi_grid_sector, r_grid)\n",
    "\n",
    "E_ref = FDFD_ref.get_TM_field(J_grid)\n",
    "E_IPML_thin = FDFD_IPML_thin.get_TM_field(J_grid)\n",
    "E_IPML_thick = FDFD_IPML_thick.get_TM_field(J_grid)\n",
    "\n",
    "print('compare fields')\n",
    "plot_cplx_polar_field(E_ref, phi_grid_sector, r_grid)\n",
    "plot_cplx_polar_field(E_IPML_thin, phi_grid_sector, r_grid)\n",
    "plot_cplx_polar_field(E_IPML_thick, phi_grid_sector, r_grid)\n",
    "\n",
    "print('plot fields with various IPML thickness along r')\n",
    "plt.figure()\n",
    "plt.plot(r_grid, E_IPML_thin[:Nr], label='IPML 0.5')\n",
    "plt.plot(r_grid, E_IPML_thick[:Nr], label='IPML 1.0')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6",
   "metadata": {},
   "source": [
    "Inner pml works fine for these m=0 waves.\n",
    "\n",
    "### off center dipoles.\n",
    "These excite all waves with m an integer factor of n_sectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# current source: 6 fold rotational point sources\n",
    "J_rgrid = np.zeros(Nr, dtype=complex)\n",
    "J_rgrid[int(r0 / dr)] = 1.0 / (2*np.pi*r0*dr) * Nphi_tot\n",
    "J_phigrid = np.zeros(Nphi_sector)\n",
    "J_phigrid[Nphi_sector // 2] = 1.0\n",
    "J_grid = np.kron(J_phigrid, J_rgrid)\n",
    "\n",
    "plot_real_polar_field(np.real(J_grid), phi_grid_sector, r_grid)\n",
    "\n",
    "E_ref = FDFD_ref.get_TM_field(J_grid)\n",
    "E_IPML_thin = FDFD_IPML_thin.get_TM_field(J_grid)\n",
    "E_IPML_thick = FDFD_IPML_thick.get_TM_field(J_grid)\n",
    "\n",
    "print('compare fields')\n",
    "plot_cplx_polar_field(E_ref, phi_grid_sector, r_grid)\n",
    "plot_cplx_polar_field(E_IPML_thin, phi_grid_sector, r_grid)\n",
    "plot_cplx_polar_field(E_IPML_thick, phi_grid_sector, r_grid)\n",
    "\n",
    "print('plot fields with various IPML thickness along r')\n",
    "plt.figure()\n",
    "plt.plot(r_grid, E_IPML_thin[:Nr], label='IPML 0.5')\n",
    "plt.plot(r_grid, E_IPML_thick[:Nr], label='IPML 1.0')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "#### THESE RESULTS SUGGEST INNER - PML SHOULD BE USED WITH CARE: INCLUDING ORIGIN MAY CAUSE ISSUES\n",
    "\n",
    "### A cosine ring current source \n",
    "Excites only m=n_sector waves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# try a current source that only has a m=6 rotational order\n",
    "J_rgrid = np.zeros(Nr, dtype=complex)\n",
    "J_rgrid[int(r0 / dr)] = 1.0 / (2*np.pi*r0*dr)\n",
    "J_phigrid = np.cos(n_sectors * phi_grid_sector)\n",
    "J_grid = np.kron(J_phigrid, J_rgrid)\n",
    "\n",
    "plot_cplx_polar_field(np.real(J_grid), phi_grid_sector, r_grid)\n",
    "\n",
    "E_ref = FDFD_ref.get_TM_field(J_grid)\n",
    "E_IPML_thin = FDFD_IPML_thin.get_TM_field(J_grid)\n",
    "E_IPML_thick = FDFD_IPML_thick.get_TM_field(J_grid)\n",
    "\n",
    "print('compare fields')\n",
    "plot_cplx_polar_field(E_ref, phi_grid_sector, r_grid)\n",
    "plot_cplx_polar_field(E_IPML_thin, phi_grid_sector, r_grid)\n",
    "plot_cplx_polar_field(E_IPML_thick, phi_grid_sector, r_grid)\n",
    "\n",
    "print('plot fields with various IPML thickness along r')\n",
    "plt.figure()\n",
    "plt.plot(r_grid, E_IPML_thin[:Nr], label='IPML 0.5')\n",
    "plt.plot(r_grid, E_IPML_thick[:Nr], label='IPML 1.0')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "These results suggest that some issues remain with inner PML for spherical waves with order m > 0.\n",
    "\n",
    "## Domain far from origin\n",
    "\n",
    "The prior results were with the default r_inner = 0. The following tests are for non-zero r_inner, to see if starting at larger radii improves inner PML performance.\n",
    "\n",
    "### setting up geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11",
   "metadata": {},
   "outputs": [],
   "source": [
    "#r_inner = 0.1\n",
    "r_inner = 5.0\n",
    "r_center = 3.0\n",
    "w_pml_outer = 0.5\n",
    "w_pml_inner_thin = 0.5\n",
    "w_pml_inner_thick = 1.0\n",
    "\n",
    "r_delta = 2*w_pml_inner_thin + r_center\n",
    "r_outer = r_inner + r_delta\n",
    "\n",
    "gpr = 40\n",
    "dr = 1.0 / gpr\n",
    "Nr = int(r_delta / dr)\n",
    "\n",
    "n_sectors = 6\n",
    "Nphi_sector = int(2*np.pi*r_outer / n_sectors / dr)\n",
    "\n",
    "Npml_outer = int(w_pml_outer / dr)\n",
    "Npml_inner_thin = int(w_pml_inner_thin / dr)\n",
    "Npml_inner_thick = int(w_pml_inner_thick / dr)\n",
    "\n",
    "# thin inner PML\n",
    "w_pml_inner_thin = 0.5\n",
    "Npml_inner_thin = int(np.round(w_pml_inner_thin / dr))\n",
    "print('Npml_inner_thin', Npml_inner_thin)\n",
    "geo_IPML_thin = PolarFDFDGeometry(Nphi_sector, Nr, Npml_outer, dr, \n",
    "                                  n_sectors=n_sectors,\n",
    "                                  r_inner=r_inner,\n",
    "                                  Npml_inner=Npml_inner_thin)\n",
    "FDFD_IPML_thin = TM_Polar_FDFD(omega0, geo_IPML_thin)\n",
    "\n",
    "# thick inner PML\n",
    "w_pml_inner_thick = 1.0\n",
    "Npml_inner_thick = int(np.round(w_pml_inner_thick / dr))\n",
    "geo_IPML_thick = PolarFDFDGeometry(Nphi_sector, Nr, Npml_outer, dr, \n",
    "                                   n_sectors=n_sectors,\n",
    "                                   r_inner=r_inner,\n",
    "                                   Npml_inner=Npml_inner_thick)\n",
    "FDFD_IPML_thick = TM_Polar_FDFD(omega0, geo_IPML_thick)\n",
    "\n",
    "# get coordinate grids\n",
    "phi_grid_sector, r_grid, phi_grid_full = FDFD_IPML_thin.get_symmetric_grids()\n",
    "print(FDFD_IPML_thin.Npml_inner)\n",
    "print(FDFD_IPML_thick.Npml_inner)\n",
    "print('Npml', Npml)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "J_r_ind = Npml_inner_thick + gpr // 4\n",
    "J_m = n_sectors * 4 # adjust wave order here\n",
    "J_rgrid = np.zeros(Nr, dtype=complex)\n",
    "J_rgrid[J_r_ind] = 1.0 / (2*np.pi*r_grid[J_r_ind]*dr)\n",
    "J_phigrid = np.cos(phi_grid_sector * J_m)\n",
    "J_grid = np.kron(J_phigrid, J_rgrid)\n",
    "\n",
    "plot_cplx_polar_field(J_grid, phi_grid_sector, r_grid)\n",
    "\n",
    "E_IPML_thin = FDFD_IPML_thin.get_TM_field(J_grid)\n",
    "E_IPML_thick = FDFD_IPML_thick.get_TM_field(J_grid)\n",
    "\n",
    "plot_cplx_polar_field(E_IPML_thin, phi_grid_sector, r_grid)\n",
    "plot_cplx_polar_field(E_IPML_thick, phi_grid_sector, r_grid)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(r_grid, np.real(E_IPML_thin[:Nr]), label='IPML 0.5')\n",
    "plt.plot(r_grid, np.real(E_IPML_thick[:Nr]), label='IPML 1.0')\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(r_grid, np.imag(E_IPML_thin[:Nr]), label='IPML 0.5')\n",
    "plt.plot(r_grid, np.imag(E_IPML_thick[:Nr]), label='IPML 1.0')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "To summarize: inner pml seems to be working well for r_inner large ($>0.5\\lambda$). If the origin is included some wonky behavior happens for larger angular momentum order waves. "
   ]
  }
 ],
 "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}