Polar TM Maxwell Solver Test / Tutorial with mirror symmetry and discrete rotational symmetry

Sanity checks for correctness of the polar TM Maxwell solver with discrete rotational symmetry and mirror symmetry.

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 6-fold discrete rotational symmetry and mirrr symmetry, which allows us to use the n_sector and mirror feature of the polar solver to restrict computation within the irreducible domain. In this case the irreducible domain is a wedge of angle \(30^\circ\), which under reflection forms one 6-fold 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 # center circle 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: 6-fold rotational symmetry and mirror symmetry (60° sectors, 30° half sectors)
n_sectors = 6
mirror = True
assert Nphi_full % (2*n_sectors) == 0, "Nphi_full must be divisible by 2*n_sectors"

Nphi_sector = Nphi_full // n_sectors
Nphi_halfsector = Nphi_sector // 2  # azimuthal points in one 30° irreducible domain

geo_halfsector = PolarFDFDGeometry(Nphi_halfsector, Nr, Npml, dr,
                                   n_sectors=n_sectors, mirror=True)

FDFD_halfsector = TM_Polar_FDFD(omega, geo_halfsector)

# Get coordinate grids
phi_grid_halfsector, r_grid, phi_grid_full = FDFD_halfsector.get_symmetric_grids()

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_halfsector instead of Nphi_full for the symmetric case
J_center_dipole = np.kron(np.ones(Nphi_halfsector), J_r)
# note the ordering of the flattened fields: all r for one ray then move counter-clockwise

E_center_dipole = FDFD_halfsector.get_TM_field(J_center_dipole) # get the center dipole field
plot_cplx_polar_field(E_center_dipole, phi_grid_halfsector, 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')

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

Putting in structures

Do a test with a non-vacuum structure and dipole in center. Structure in irreducible domain: annular pie with \(\phi \in [0,\pi/12]\), \(r \in [1.0,1.3]\), \(\chi=3+0.1i\)

## setup structure: ring at radius r_ring with thickness r_t
r_ring = 1.0  # center radius of ring
r_t = 0.3     # thickness of ring
r_inner = r_ring
r_outer = r_ring + r_t

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

## setup structure grid
chi_phi_grid = np.zeros(Nphi_halfsector, dtype=complex)
phi_start = 0  # start at phi = 0
phi_end = int(np.round(Nphi_halfsector / 2))   # end at 1/4 of a sector, i.e. 1/2 of a half sector, or pi/12 in this case
chi_phi_grid[phi_start:phi_end] = 1.0

chi_grid = np.kron(chi_phi_grid, chi_r_grid)
print("structure in irreducible domain:")
plot_real_polar_field(np.real(chi_grid), phi_grid_halfsector, r_grid, cmap='gist_yarg') # mask of structure over irreducible halfsector

chi_grid_full = expand_symmetric_field(chi_grid, n_sectors, Nr, mirror=True)
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_halfsector.get_TM_field(J_center_dipole, chi_grid)
plot_cplx_polar_field(E_struct, phi_grid_halfsector, r_grid)

print("Full circle visualization:")
# important to specify mirror=True
E_struct_full = expand_symmetric_field(E_struct, n_sectors, Nr, mirror=True)
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,
                                  mirror=True)

# 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)

ind_r = int(source_r / dr)
ind_phi = int(source_phi_fraction * Nphi_halfsector)

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

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

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

# Calculate pixel area at this radius
pixel_area = r_grid[ind_r] * dr * (2*np.pi / n_sectors / 2 / Nphi_halfsector)
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_halfsector.get_TM_field(J_dipoles)

E_dipoles_full = expand_symmetric_field(E_dipoles, n_sectors, Nr, mirror=True)
print("Full circle - showing D6 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, mirror=True)
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}")