Cavity LDOS Bounds: Cartesian vs Polar Coordinates

This notebook demonstrates how to compute limits using the Polar coordinate formulation in Dolphindes. To verify consistency, we also compare the results with the Cartesian coordinate formulation.

%load_ext autoreload
%autoreload 2
import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt
import sys, time, os
sys.path.append("../../")

from dolphindes import photonics, geometry
from dolphindes.cvxopt import gcd
from dolphindes.maxwell import plot_real_polar_field, plot_cplx_polar_field, expand_symmetric_field, TM_Polar_FDFD

Problem Setup

We consider a 2D TM problem. We want to maximize the LDOS at the center of a circular design region. The design region is an annulus \(R_{in} < r < R_{out}\). We place a point source at the center.

wavelength = 1.0
omega = 2 * np.pi / wavelength
chi = 4 + 1e-4j

# Geometry parameters
R_inner = 0.2
R_outer = 0.6
Npml = 10
gpr = 50
dl = 1/gpr

1. Cartesian Coordinates

This should look pretty much identical to prior tutorials. We are just getting the cartesian value.

# Grid setup
L_domain = 2 * (R_outer + 0.5)
Nx = Ny = int(L_domain / dl) + 2 * Npml

# Source at center
ji_cart = np.zeros((Nx, Ny), dtype=complex)
cx, cy = Nx // 2, Ny // 2
ji_cart[cx, cy] = 1.0 / dl**2

# Design mask (Annulus)
x = (np.arange(Nx) - cx) * dl
y = (np.arange(Ny) - cy) * dl
X, Y = np.meshgrid(x, y, indexing='ij')
R = np.sqrt(X**2 + Y**2)
des_mask_cart = (R >= R_inner) & (R <= R_outer)
ndof = np.sum(des_mask_cart)

chi_bg_cart = np.zeros((Nx, Ny), dtype=complex)
plt.imshow(des_mask_cart.T + np.real(ji_cart.T).astype(bool), origin='lower', cmap='Greys')

# Solver Setup
geo_cart = geometry.CartesianFDFDGeometry(Nx, Ny, Npml, Npml, dl, dl)
prob_cart = photonics.Photonics_TM_FDFD(
    omega, geo_cart, chi, des_mask_cart, ji_cart, chi_background=chi_bg_cart, sparseQCQP=True
)

ei_cart = prob_cart.get_ei(ji_cart, update=True)
vac_ldos_cart = -np.sum(0.5 * np.real(ji_cart.conj() * ei_cart) * dl**2)
plt.imshow(ei_cart.T.real, origin='lower', cmap='RdBu')
plt.colorbar()
print(f"Cartesian Vacuum LDOS: {vac_ldos_cart:.4f}")

# Objective: Maximize LDOS
ei_des_cart = ei_cart[des_mask_cart]
s0_cart = -0.25 * 1j * omega * ei_des_cart.conj()
ndof = np.sum(des_mask_cart)
A0_cart = sp.csc_array(np.zeros((ndof, ndof), dtype=complex))
prob_cart.set_objective(A0=A0_cart, s0=s0_cart, c0=vac_ldos_cart, denseToSparse=True)

# QCQP Setup
prob_cart.setup_QCQP(Pdiags='global', verbose=1)
res_cart = prob_cart.bound_QCQP(method='bfgs')
print(f"Cartesian Bound: {res_cart[0]:.4f}")
print(f"Cartesian Enhancement Bound: {res_cart[0]/vac_ldos_cart:.4f}")

2. Polar Coordinates

wvlgth = 1.0
Qabs = np.inf # supports complex frequency, test by setting finite Qabs
omega = 2*np.pi / wvlgth * (1 + 1j/2/Qabs)

R_nonpml = 3.0 # center circle radius
w_pml = 0.5 # surrounding pml thickness
R_tot = R_nonpml + w_pml # total computational domain radius
R_i = R_inner # inner radius of the cavity
R_o = R_outer # outer radius of the cavity
assert R_o < R_nonpml, "Outer radius must be within non-PML region"

gpr = 20
dr = 1.0/gpr # radial grid size
Nr = int(np.round(R_tot / dr))
Npml = int(np.round(w_pml / dr))
Nr_i = int(np.round(R_i / dr)) # inner radius grid point
Nr_o = int(np.round(R_o / dr)) # outer radius grid point
Nphi = 100
n_sectors = 1
Nphi_sector = int(Nphi/n_sectors)  # azimuthal points in one sector

assert Nphi % n_sectors == 0, "Nphi must be divisible by n_sectors"

# Normalize source so integrated current in the central pixel = 1
J_r = np.zeros(Nr)
J_r[0] = 1.0 / (np.pi * dr**2)  # central pixel current density = 1/area
J = np.kron(np.ones(Nphi_sector), J_r)
J_full = expand_symmetric_field(J, n_sectors, Nr)
phi_grid = np.linspace(0, 2*np.pi, Nphi, endpoint=False)
r_grid = (np.arange(Nr) + 0.5) * dr

# Design Mask
des_mask_polar = np.zeros((Nr, Nphi), dtype=bool)
for i, r in enumerate(r_grid):
    if R_inner <= r <= R_outer:
        des_mask_polar[i, :] = True
ndof = np.sum(des_mask_polar)

# Use radial spacing 'dr' for polar geometry
geo_polar = geometry.PolarFDFDGeometry(Nphi, Nr, Npml, dr)

# Solver Setup
prob_polar = photonics.Photonics_TM_FDFD(
    omega, geo_polar, chi, des_mask_polar, J, sparseQCQP=True
)

prob_polar.setup_EM_solver()
phi_grid_sector, r_grid, phi_grid_full = prob_polar.EM_solver.get_symmetric_grids()
areas = prob_polar.geometry.get_pixel_areas()

ei = prob_polar.get_ei()
plot_cplx_polar_field(np.real(ei), phi_grid, r_grid)
# Evaluate vacuum LDOS (per-sector J, multiply by n_sectors)
vac_ldos_polar = -0.5 * np.real(J.conj().dot(areas * ei)) * n_sectors
print(f'vac_ldos_polar: {vac_ldos_polar}')
# print(f'vac_ldos_cart: {vac_ldos_cart}')

plot_real_polar_field(prob_polar.des_mask, phi_grid_sector, r_grid, cmap='Greys')

ei_des_polar = ei[prob_polar.des_mask.flatten(order='F')]
s0_polar = -0.25 * 1j * omega * ei_des_polar.conj()
ndof = np.sum(prob_polar.des_mask)
A0_polar = sp.csc_array(np.zeros((ndof, ndof), dtype=complex))
prob_polar.set_objective(A0=A0_polar, s0=s0_polar, c0=vac_ldos_polar, denseToSparse=True)

prob_polar.setup_QCQP(Pdiags='global', verbose=1)
res_polar = prob_polar.bound_QCQP(method='bfgs')
print(f"Polar Bound: {res_polar[0]:.4f}")
print(f"Polar Enhancement Bound: {res_polar[0]/vac_ldos_polar:.4f}")

3. Comparison and GCD

We expect the enhancement bounds (Purcell factors) to be similar, though differences in discretization near the singularity (source) will cause some discrepancy in the absolute vacuum LDOS.

print(f"Cartesian Enhancement: {res_cart[0]/vac_ldos_cart:.4f}")
print(f"Polar Enhancement:     {res_polar[0]/vac_ldos_polar:.4f}")

Running GCD

Now we run General Constraint Descent to tighten the bounds for both cases.

gcd_params = gcd.GCDHyperparameters(
    max_proj_cstrt_num=16,
    max_gcd_iter_num=20,
    gcd_iter_period=5,
    gcd_tol=1e-3
)

print("Running GCD on Cartesian...")
t0 = time.time()
prob_cart.QCQP.run_gcd(gcd_params)
print(f"Cartesian GCD Time: {time.time()-t0:.2f}s")
print(f"Cartesian GCD Bound: {prob_cart.QCQP.current_dual:.4f}")
print(f"Cartesian GCD Enhancement: {prob_cart.QCQP.current_dual/vac_ldos_cart:.4f}")

print("\nRunning GCD on Polar...")
t0 = time.time()
prob_polar.QCQP.run_gcd(gcd_params)
print(f"Polar GCD Time: {time.time()-t0:.2f}s")
print(f"Polar GCD Bound: {prob_polar.QCQP.current_dual:.4f}")
print(f"Polar GCD Enhancement: {prob_polar.QCQP.current_dual/vac_ldos_polar:.4f}")