Polar Solver Inner PML test

Checking the veracity and usefullness of having inner PML in the polar coordinate Helmholtz solver.

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

Computational domain starting at origin

The following results suggest that the current inner pml implementation can be problematic when the computational domain starts at the origin.

setting up geometry

# domain parameters
wvlgth0 = 1.0 # center wavelength of interest
omega0 = 2*np.pi / wvlgth0 # center frequency of interest

r_tot = 4.0 #  computational domain outer radius (includes all pml)
w_pml_outer = 0.5 # thickness of outer layer of pml
gpr = 80
dr = 1.0 / gpr
Nr = int(np.round(r_tot / dr))
Npml = int(np.round(w_pml_outer/dr))

Nphi_tot = 720 # total azimuthal gridlines around full circle
n_sectors = 6
Nphi_sector = Nphi_tot // n_sectors

# no inner PML solver, for reference
geo_ref = PolarFDFDGeometry(Nphi_sector, Nr, Npml, dr, n_sectors=n_sectors)
FDFD_ref = TM_Polar_FDFD(omega0, geo_ref)

# thin inner PML
w_pml_inner_thin = 0.5
Npml_inner_thin = int(np.round(w_pml_inner_thin / dr))
geo_IPML_thin = PolarFDFDGeometry(Nphi_sector, Nr, Npml, dr,
                                  n_sectors=n_sectors,
                                  Npml_inner=Npml_inner_thin)
FDFD_IPML_thin = TM_Polar_FDFD(omega0, geo_IPML_thin)

# thick inner PML
w_pml_inner_thick = 1.0
Npml_inner_thick = int(np.round(w_pml_inner_thick / dr))
geo_IPML_thick = PolarFDFDGeometry(Nphi_sector, Nr, Npml, dr,
                                  n_sectors=n_sectors,
                                  Npml_inner=Npml_inner_thick)
FDFD_IPML_thick = TM_Polar_FDFD(omega0, geo_IPML_thick)

# get coordinate grids
phi_grid_sector, r_grid, phi_grid_full = FDFD_ref.get_symmetric_grids()

uniform ring source

excites only m=0 waves.

# current source: a uniform delta(r - r0)
r0 = r_tot / 2
J_rgrid = np.zeros(Nr, dtype=complex)
J_rgrid[int(r0 / dr)] = 1.0 / (2*np.pi*r0*dr)
J_grid = np.kron(np.ones(Nphi_sector), J_rgrid)

print('current ring')
plot_real_polar_field(np.real(J_grid), phi_grid_sector, r_grid)

E_ref = FDFD_ref.get_TM_field(J_grid)
E_IPML_thin = FDFD_IPML_thin.get_TM_field(J_grid)
E_IPML_thick = FDFD_IPML_thick.get_TM_field(J_grid)

print('compare fields')
plot_cplx_polar_field(E_ref, phi_grid_sector, r_grid)
plot_cplx_polar_field(E_IPML_thin, phi_grid_sector, r_grid)
plot_cplx_polar_field(E_IPML_thick, phi_grid_sector, r_grid)

print('plot fields with various IPML thickness along r')
plt.figure()
plt.plot(r_grid, E_IPML_thin[:Nr], label='IPML 0.5')
plt.plot(r_grid, E_IPML_thick[:Nr], label='IPML 1.0')
plt.show()

Inner pml works fine for these m=0 waves.

off center dipoles.

These excite all waves with m an integer factor of n_sectors.

# current source: 6 fold rotational point sources
J_rgrid = np.zeros(Nr, dtype=complex)
J_rgrid[int(r0 / dr)] = 1.0 / (2*np.pi*r0*dr) * Nphi_tot
J_phigrid = np.zeros(Nphi_sector)
J_phigrid[Nphi_sector // 2] = 1.0
J_grid = np.kron(J_phigrid, J_rgrid)

plot_real_polar_field(np.real(J_grid), phi_grid_sector, r_grid)

E_ref = FDFD_ref.get_TM_field(J_grid)
E_IPML_thin = FDFD_IPML_thin.get_TM_field(J_grid)
E_IPML_thick = FDFD_IPML_thick.get_TM_field(J_grid)

print('compare fields')
plot_cplx_polar_field(E_ref, phi_grid_sector, r_grid)
plot_cplx_polar_field(E_IPML_thin, phi_grid_sector, r_grid)
plot_cplx_polar_field(E_IPML_thick, phi_grid_sector, r_grid)

print('plot fields with various IPML thickness along r')
plt.figure()
plt.plot(r_grid, E_IPML_thin[:Nr], label='IPML 0.5')
plt.plot(r_grid, E_IPML_thick[:Nr], label='IPML 1.0')
plt.show()

THESE RESULTS SUGGEST INNER - PML SHOULD BE USED WITH CARE: INCLUDING ORIGIN MAY CAUSE ISSUES

A cosine ring current source

Excites only m=n_sector waves

# try a current source that only has a m=6 rotational order
J_rgrid = np.zeros(Nr, dtype=complex)
J_rgrid[int(r0 / dr)] = 1.0 / (2*np.pi*r0*dr)
J_phigrid = np.cos(n_sectors * phi_grid_sector)
J_grid = np.kron(J_phigrid, J_rgrid)

plot_cplx_polar_field(np.real(J_grid), phi_grid_sector, r_grid)

E_ref = FDFD_ref.get_TM_field(J_grid)
E_IPML_thin = FDFD_IPML_thin.get_TM_field(J_grid)
E_IPML_thick = FDFD_IPML_thick.get_TM_field(J_grid)

print('compare fields')
plot_cplx_polar_field(E_ref, phi_grid_sector, r_grid)
plot_cplx_polar_field(E_IPML_thin, phi_grid_sector, r_grid)
plot_cplx_polar_field(E_IPML_thick, phi_grid_sector, r_grid)

print('plot fields with various IPML thickness along r')
plt.figure()
plt.plot(r_grid, E_IPML_thin[:Nr], label='IPML 0.5')
plt.plot(r_grid, E_IPML_thick[:Nr], label='IPML 1.0')
plt.show()

These results suggest that some issues remain with inner PML for spherical waves with order m > 0.

Domain far from origin

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.

setting up geometry

#r_inner = 0.1
r_inner = 5.0
r_center = 3.0
w_pml_outer = 0.5
w_pml_inner_thin = 0.5
w_pml_inner_thick = 1.0

r_delta = 2*w_pml_inner_thin + r_center
r_outer = r_inner + r_delta

gpr = 40
dr = 1.0 / gpr
Nr = int(r_delta / dr)

n_sectors = 6
Nphi_sector = int(2*np.pi*r_outer / n_sectors / dr)

Npml_outer = int(w_pml_outer / dr)
Npml_inner_thin = int(w_pml_inner_thin / dr)
Npml_inner_thick = int(w_pml_inner_thick / dr)

# thin inner PML
w_pml_inner_thin = 0.5
Npml_inner_thin = int(np.round(w_pml_inner_thin / dr))
print('Npml_inner_thin', Npml_inner_thin)
geo_IPML_thin = PolarFDFDGeometry(Nphi_sector, Nr, Npml_outer, dr,
                                  n_sectors=n_sectors,
                                  r_inner=r_inner,
                                  Npml_inner=Npml_inner_thin)
FDFD_IPML_thin = TM_Polar_FDFD(omega0, geo_IPML_thin)

# thick inner PML
w_pml_inner_thick = 1.0
Npml_inner_thick = int(np.round(w_pml_inner_thick / dr))
geo_IPML_thick = PolarFDFDGeometry(Nphi_sector, Nr, Npml_outer, dr,
                                   n_sectors=n_sectors,
                                   r_inner=r_inner,
                                   Npml_inner=Npml_inner_thick)
FDFD_IPML_thick = TM_Polar_FDFD(omega0, geo_IPML_thick)

# get coordinate grids
phi_grid_sector, r_grid, phi_grid_full = FDFD_IPML_thin.get_symmetric_grids()
print(FDFD_IPML_thin.Npml_inner)
print(FDFD_IPML_thick.Npml_inner)
print('Npml', Npml)

J_r_ind = Npml_inner_thick + gpr // 4
J_m = n_sectors * 4 # adjust wave order here
J_rgrid = np.zeros(Nr, dtype=complex)
J_rgrid[J_r_ind] = 1.0 / (2*np.pi*r_grid[J_r_ind]*dr)
J_phigrid = np.cos(phi_grid_sector * J_m)
J_grid = np.kron(J_phigrid, J_rgrid)

plot_cplx_polar_field(J_grid, phi_grid_sector, r_grid)

E_IPML_thin = FDFD_IPML_thin.get_TM_field(J_grid)
E_IPML_thick = FDFD_IPML_thick.get_TM_field(J_grid)

plot_cplx_polar_field(E_IPML_thin, phi_grid_sector, r_grid)
plot_cplx_polar_field(E_IPML_thick, phi_grid_sector, r_grid)

plt.figure()
plt.plot(r_grid, np.real(E_IPML_thin[:Nr]), label='IPML 0.5')
plt.plot(r_grid, np.real(E_IPML_thick[:Nr]), label='IPML 1.0')
plt.show()

plt.figure()
plt.plot(r_grid, np.imag(E_IPML_thin[:Nr]), label='IPML 0.5')
plt.plot(r_grid, np.imag(E_IPML_thick[:Nr]), label='IPML 1.0')
plt.show()

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.