Polar TM Maxwell Solver Test / Tutorial with discrete rotational symmetry

Sanity checks for correctness of the polar TM Maxwell solver with discrete rotational symmetry, which also demonstrates its basic usage. High level API is essentially similar to original TM_FDFD. NOTE: because of the polar coordinates, the FD pixels are no longer of unit size.

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
import matplotlib.pyplot as plt
import matplotlib.colors as colors

import sys
sys.path.append("../../")

from dolphindes.geometry import PolarFDFDGeometry
from dolphindes.maxwell import TM_Polar_FDFD, plot_real_polar_field, plot_cplx_polar_field, expand_symmetric_field

We start by specifying a domain of interest shared between all basic tests below: a circular computational domain with radius \(2\lambda\), surrounded by pml. We will consider only structures and sources with 12-fold discrete rotational symmetry, which allows us to use the n_sector feature of the polar solver to restrict computation within the rotational unit cell.

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 = 2.0 # computational domain radius
w_pml = 0.5 # surrounding pml thickness
r_tot = r_nonpml + w_pml # total computational domain radius

gpr = 30
#gpr = 60 # high res option
dr = 1.0/gpr # radial grid size
Nr = int(np.round(r_tot / dr))
Npml = int(np.round(w_pml / dr))

# setting azimuthal grid size. Note the azimuthal pixel width gets larger with radius
Nphi_full = 180 # this gives ~ 0.043 pixel width at the edge of computational domain for r_tot=2.5
#Nphi_full = 360 # high res option

# Example: 12-fold rotational symmetry (30° sectors)
n_sectors = 12

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

Nphi_sector = Nphi_full // n_sectors  # azimuthal points in one sector

# geometry
# Important: n_sectors=1 by default so make sure to specify otherwise
geo_sector = PolarFDFDGeometry(Nphi_sector, Nr, Npml, dr,
                               n_sectors=n_sectors)

# FDFD solver
FDFD_sector = TM_Polar_FDFD(omega, geo_sector)

# Get coordinate grids for plotting later
phi_grid_sector, r_grid, phi_grid_full = FDFD_sector.get_symmetric_grids()
print(phi_grid_sector.shape)
print(phi_grid_sector)
print(FDFD_sector.phi_grid)

Vacuum Dipole

from scipy.special import hankel1 # for checking against analytical dipole field

vacuum dipole at the center

The simplest test: a dipole source at the origin. Because of the radial grid offset there is no gridpoint exactly at the origin, so we specify the dipole using all the points closest to the origin.

# setup dipole source at origin
J_r = np.zeros(Nr)
J_r[0] = 1.0 / (np.pi* dr**2)
# note the normalization: each r gridpoint r_i covers the range [r_i - dr/2, r_i + dr/2)

# Use Nphi_sector instead of Nphi_full for the symmetric case
J_center_dipole = np.kron(np.ones(Nphi_sector), J_r)
# note the ordering of the flattened fields: all r for one ray then move counter-clockwise

E_center_dipole = FDFD_sector.get_TM_field(J_center_dipole)
plot_cplx_polar_field(E_center_dipole, phi_grid_sector, r_grid) # use sector grid for plotting

# check against analytical solution along phi=0 ray (first ray in sector)
analytic_E_grid = (-omega/4) * hankel1(0, omega*r_grid)
fig, (ax1,ax2) = plt.subplots(ncols=2, figsize=(10,5))
ax1.plot(r_grid, np.real(E_center_dipole[:Nr]), label='FD real(E)')
ax1.plot(r_grid, np.real(analytic_E_grid), '--', label='analytic real(E)')
ax1.set_xlabel('r'); ax1.legend()
ax2.plot(r_grid, np.imag(E_center_dipole[:Nr]), label='FD imag(E)')
ax2.plot(r_grid, np.imag(analytic_E_grid), '--', label='analytic imag(E)')
ax2.set_xlabel('r'); ax2.legend()
plt.show()
print(f'deviation starting around r={r_nonpml} due to pml')

# expand to full circle for visualization
E_center_dipole_full = expand_symmetric_field(E_center_dipole, n_sectors, Nr)
print("Full circle visualization:")
plot_cplx_polar_field(E_center_dipole_full, phi_grid_full, r_grid)

Putting in structures

Now try a non-vacuum structure and dipole in center. rotational unit cell in \(\phi \in [0, \pi / 6)\): annular pie with \(\phi \in [\pi/24,\pi/8]\), \(r \in [1.0,1.5]\), \(\chi=3+0.1i\)

r_ring_inner = 1.0
r_ring_outer =1.5

# Create radial mask for the ring
chi_r_grid = np.zeros(Nr)
chi_r_grid[int(r_ring_inner/dr):int(r_ring_outer/dr)] = 1.0

## setup structure grid
chi_phi_grid = np.zeros(Nphi_sector, dtype=complex)
phi_start = int(Nphi_sector * 0.25)  # start at 25% of sector
phi_end = int(Nphi_sector * 0.75)    # end at 75% of sector
chi_phi_grid[phi_start:phi_end] = 1.0
chi_grid = np.kron(chi_phi_grid, chi_r_grid)
plot_real_polar_field(np.real(chi_grid), phi_grid_sector, r_grid, cmap='gist_yarg') # mask of structure

chi_grid_full = expand_symmetric_field(chi_grid, n_sectors, Nr)
print("Full circle visualization:")
plot_real_polar_field(np.real(chi_grid_full), phi_grid_full, r_grid, cmap='gist_yarg')

chi = 3.0 + 0.1j
chi_grid *= chi

E_struct = FDFD_sector.get_TM_field(J_center_dipole, chi_grid)
# plot_cplx_polar_field(E, phi_grid_sector, r_grid)

print("Full circle visualization:")
E_struct_full = expand_symmetric_field(E_struct, n_sectors, Nr)
plot_cplx_polar_field(E_struct_full, phi_grid_full, r_grid)

## Verification with polar TM FDFD solver without symmetry

# Create full-circle FDFD solver
geo_full = PolarFDFDGeometry(Nphi_full, Nr, Npml, dr)
FDFD_full = TM_Polar_FDFD(omega, geo_full)

# Setup dipole source at center for full circle
J_r_full = np.zeros(Nr)
J_r_full[0] = 1.0 / (np.pi* dr**2)
J_full = np.kron(np.ones(Nphi_full), J_r_full)

# Setup structure for full circle (expand the symmetric structure)
chi_full = expand_symmetric_field(chi_grid, n_sectors, Nr)

# Solve the full problem
E_full_reference = FDFD_full.get_TM_field(J_full, chi_full)

# Plot the full reference solution
phi_grid_reference = np.linspace(0, 2*np.pi, Nphi_full, endpoint=False)
plot_cplx_polar_field(E_full_reference, phi_grid_reference, r_grid)

# Check specific points
print("\nField comparison at origin:")
print(f"Full solver: {E_full_reference[0]:.6f}")
print(f"Symmetric solver: {E_struct_full[0]:.6f}")

# Check a point in the middle of the domain
mid_point = Nphi_full//2 * Nr + Nr//2
print(f"\nField comparison at middle point:")
print(f"Full solver: {E_full_reference[mid_point]:.6f}")
print(f"Symmetric solver: {E_struct_full[mid_point]:.6f}")

Off-center dipoles

One off-center dipole source in each rotational unit cell.

# Source parameters
source_r = 1.0  # radial position
source_phi_fraction = 0.5  # fraction through the sector (0.5 = middle)

# Find grid indices
ind_r = int(source_r / dr)
ind_phi = int(source_phi_fraction * Nphi_sector)

print(f"Source at r = {r_grid[ind_r]:.3f}, phi = {phi_grid_sector[ind_phi]*180/np.pi:.1f}°")

# Setup source vector for one sector
J_dipoles = np.zeros(Nphi_sector * Nr)

# Calculate pixel area at this radius
pixel_area = r_grid[ind_r] * dr * (2*np.pi / n_sectors / Nphi_sector)
source_strength = 1.0 / pixel_area

# Place source at the specified location in the sector
J_dipoles[ind_phi * Nr + ind_r] = source_strength

print(f"Source strength: {source_strength:.2e}")
print(f"Pixel area: {pixel_area:.6f}")

# Solve with vacuum Maxwell operator (no structure)
E_dipoles = FDFD_sector.get_TM_field(J_dipoles)

E_dipoles_full = expand_symmetric_field(E_dipoles, n_sectors, Nr)
print("Full circle - showing 12-fold symmetry:")
plot_cplx_polar_field(E_dipoles_full, phi_grid_full, r_grid)

# Verification: compare with full solver
J_dipoles_full = expand_symmetric_field(J_dipoles, n_sectors, Nr)

# Solve full problem
E_dipoles_full_ref = FDFD_full.get_TM_field(J_dipoles_full)

# Compare solutions
print("Full solver reference:")
plot_cplx_polar_field(E_dipoles_full_ref, phi_grid_full, r_grid)

# Quantitative comparison
max_diff = np.max(np.abs(E_dipoles_full_ref - E_dipoles_full))
rms_diff = np.sqrt(np.mean(np.abs(E_dipoles_full_ref - E_dipoles_full)**2))

print(f"\nMaximum difference: {max_diff:.2e}")
print(f"RMS difference: {rms_diff:.2e}")
print(f"Relative RMS error: {rms_diff/np.sqrt(np.mean(np.abs(E_dipoles_full_ref)**2)):.2e}")

# Check field at source location
source_idx_sector = ind_phi * Nr + ind_r
source_idx_full = ind_phi * Nr + ind_r  # first sector
print(f"\nField at source location:")
print(f"Symmetric solver: {E_dipoles[source_idx_sector]:.6f}")
print(f"Full solver: {E_dipoles_full_ref[source_idx_full]:.6f}")
print(f"Difference: {abs(E_dipoles[source_idx_sector] - E_dipoles_full_ref[source_idx_full]):.2e}")