Skip to content

Commit

Permalink
unit test for mode decomposition in 3d
Browse files Browse the repository at this point in the history
  • Loading branch information
@oskooi
oskooi committed Apr 5, 2022
1 parent e3578e1 commit 06d418d
Showing 1 changed file with 190 additions and 14 deletions.
204 changes: 190 additions & 14 deletions python/tests/test_mode_decomposition.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

class TestModeDecomposition(unittest.TestCase):

def test_mode_decomposition(self):
def test_linear_taper_2d(self):
resolution = 10
w1 = 1
w2 = 2
Expand All @@ -21,15 +21,15 @@ def test_mode_decomposition(self):
symmetries = [mp.Mirror(mp.Y)]
Lt = 2
sx = dpml + Lw + Lt + Lw + dpml
cell_size = mp.Vector3(sx, sy, 0)
cell_size = mp.Vector3(sx,sy,0)
prism_x = sx + 1
half_Lt = 0.5 * Lt
src_pt = mp.Vector3(-0.5 * sx + dpml + 0.2 * Lw, 0, 0)
sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen, fwidth=0.2 * fcen),
src_pt = mp.Vector3(-0.5*sx+dpml+0.2*Lw,0,0)
sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen, fwidth=0.2*fcen),
center=src_pt,
size=mp.Vector3(0, sy - 2 * dpml, 0),
size=mp.Vector3(0,sy-2*dpml,0),
eig_match_freq=True,
eig_parity=mp.ODD_Z + mp.EVEN_Y)]
eig_parity=mp.ODD_Z+mp.EVEN_Y)]

vertices = [mp.Vector3(-prism_x, half_w1),
mp.Vector3(prism_x, half_w1),
Expand All @@ -43,12 +43,18 @@ def test_mode_decomposition(self):
sources=sources,
symmetries=symmetries)

mon_pt = mp.Vector3(-0.5 * sx + dpml + 0.5 * Lw, 0, 0)
flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=mon_pt, size=mp.Vector3(0, sy - 2 * dpml, 0)))
mon_pt = mp.Vector3(-0.5*sx+dpml+0.5*Lw,0,0)
flux = sim.add_flux(fcen,
0,
1,
mp.FluxRegion(center=mon_pt,
size=mp.Vector3(0,sy-2*dpml,0)))

sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, src_pt, 1e-9))

res = sim.get_eigenmode_coefficients(flux, [1], eig_parity=mp.ODD_Z + mp.EVEN_Y)
res = sim.get_eigenmode_coefficients(flux,
[1],
eig_parity=mp.ODD_Z+mp.EVEN_Y)
incident_coeffs = res.alpha
incident_flux = mp.get_fluxes(flux)
incident_flux_data = sim.get_flux_data(flux)
Expand All @@ -71,17 +77,24 @@ def test_mode_decomposition(self):
sources=sources,
symmetries=symmetries)

refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=mon_pt, size=mp.Vector3(0, sy - 2 * dpml, 0)))
refl_flux = sim.add_flux(fcen,
0,
1,
mp.FluxRegion(center=mon_pt,
size=mp.Vector3(0,sy-2*dpml,0)))
sim.load_minus_flux_data(refl_flux, incident_flux_data)

sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, src_pt, 1e-9))

res = sim.get_eigenmode_coefficients(refl_flux, [1], eig_parity=mp.ODD_Z + mp.EVEN_Y)
res = sim.get_eigenmode_coefficients(refl_flux,
[1],
eig_parity=mp.ODD_Z+mp.EVEN_Y)
coeffs = res.alpha
taper_flux = mp.get_fluxes(refl_flux)

self.assertAlmostEqual(abs(coeffs[0, 0, 1])**2 / abs(incident_coeffs[0, 0, 0])**2,
-taper_flux[0] / incident_flux[0], places=4)
self.assertAlmostEqual(abs(coeffs[0,0,1])**2 / abs(incident_coeffs[0,0,0])**2,
-taper_flux[0] / incident_flux[0],
places=4)

def test_oblique_waveguide_backward_mode(self):
sxy = 12.0
Expand Down Expand Up @@ -118,7 +131,9 @@ def test_oblique_waveguide_backward_mode(self):
mp.FluxRegion(center=mp.Vector3(-0.5*sxy+dpml,0,0),
size=mp.Vector3(0,sxy,0)),
decimation_factor=1)
mode_decimated = sim.add_mode_monitor(fcen, 0, 1,
mode_decimated = sim.add_mode_monitor(fcen,
0,
1,
mp.FluxRegion(center=mp.Vector3(-0.5*sxy+dpml,0,0),
size=mp.Vector3(0,sxy,0)),
decimation_factor=10)
Expand All @@ -143,5 +158,166 @@ def test_oblique_waveguide_backward_mode(self):
self.assertAlmostEqual(coeff,coeff_decimated,places=3)


def test_grating_3d(self):
"""Unit test for mode decomposition in 3d.
Verifies that the reflectance and transmittance in the z
direction at a single wavelength for a unit cell of a
3d grating using a normally incident planewave is equivalent
to the sum of the Poynting flux (normalized by the flux
of the input source) for all the individual reflected
and transmitted diffracted orders.
"""
resolution = 25 # pixels/μm

nSi = 3.45
Si = mp.Medium(index=nSi)
nSiO2 = 1.45
SiO2 = mp.Medium(index=nSiO2)

wvl = 0.5 # wavelength
fcen = 1/wvl

dpml = 2.0 # PML thickness
dsub = 3.0 # substrate thickness
dair = 3.0 # air padding
hcyl = 0.5 # cylinder height
rcyl = 0.2 # cylinder radius

sx = 1.1
sy = 0.8
sz = dpml+dsub+hcyl+dair+dpml

cell_size = mp.Vector3(sx,sy,sz)

boundary_layers = [mp.PML(thickness=dpml,direction=mp.Z)]

# periodic boundary conditions
k_point = mp.Vector3()

src_cmpt = mp.Ex
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
size=mp.Vector3(sx,sy,0),
center=mp.Vector3(0,0,-0.5*sz+dpml),
component=src_cmpt)]

sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
sources=sources,
default_material=SiO2,
boundary_layers=boundary_layers,
k_point=k_point)

flux = sim.add_mode_monitor(fcen,
0,
1,
mp.ModeRegion(center=mp.Vector3(0,0,-0.5*sz+dpml+0.5*dsub),
size=mp.Vector3(sx,sy,0)))

stop_cond = mp.stop_when_fields_decayed(20,
src_cmpt,
mp.Vector3(0,0,0.5*sz-dpml-0.5*dair),
1e-6)
sim.run(until_after_sources=stop_cond)

input_flux = mp.get_fluxes(flux)
input_flux_data = sim.get_flux_data(flux)

sim.reset_meep()

geometry = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,dpml+dsub),
center=mp.Vector3(0,0,-0.5*sz+0.5*(dpml+dsub)),
material=SiO2),
mp.Cylinder(height=hcyl,
radius=rcyl,
center=mp.Vector3(0,0,-0.5*sz+dpml+dsub+0.5*hcyl),
material=Si)]

sim = mp.Simulation(resolution=resolution,
cell_size=cell_size,
sources=sources,
geometry=geometry,
boundary_layers=boundary_layers,
k_point=k_point)

refl_flux = sim.add_mode_monitor(fcen,
0,
1,
mp.ModeRegion(center=mp.Vector3(0,0,-0.5*sz+dpml+0.5*dsub),
size=mp.Vector3(sx,sy,0)))
sim.load_minus_flux_data(refl_flux,input_flux_data)

tran_flux = sim.add_mode_monitor(fcen,
0,
1,
mp.ModeRegion(center=mp.Vector3(0,0,0.5*sz-dpml),
size=mp.Vector3(sx,sy,0)))

sim.run(until_after_sources=stop_cond)

# sum the Poynting flux in z direction for all reflected orders
Rsum = 0

# number of reflected modes/orders in SiO2 in x and y directions (upper bound)
nm_x = int(fcen*nSiO2*sx) + 1
nm_y = int(fcen*nSiO2*sy) + 1
for m_x in range(nm_x):
for m_y in range(nm_y):
for S_pol in [False,True]:
res = sim.get_eigenmode_coefficients(refl_flux,
mp.DiffractedPlanewave([m_x,m_y,0],
mp.Vector3(1,0,0),
1 if S_pol else 0,
0 if S_pol else 1))
r_coeffs = res.alpha
Rmode = abs(r_coeffs[0,0,1])**2/input_flux[0]
print("refl-order:, {}, {}, {}, {:.6f}".format("s" if S_pol else "p",m_x,m_y,Rmode))
if m_x == 0 and m_y == 0:
Rsum += Rmode
elif (m_x != 0 and m_y == 0) or (m_x == 0 and m_y != 0):
Rsum += 2*Rmode
else:
Rsum += 4*Rmode


# sum the Poynting flux in z direction for all transmitted orders
Tsum = 0

# number of transmitted modes/orders in air in x and y directions (upper bound)
nm_x = int(fcen*sx) + 1
nm_y = int(fcen*sy) + 1
for m_x in range(nm_x):
for m_y in range(nm_y):
for S_pol in [False,True]:
res = sim.get_eigenmode_coefficients(tran_flux,
mp.DiffractedPlanewave([m_x,m_y,0],
mp.Vector3(1,0,0),
1 if S_pol else 0,
0 if S_pol else 1))
t_coeffs = res.alpha
Tmode = abs(t_coeffs[0,0,0])**2/input_flux[0]
print("tran-order:, {}, {}, {}, {:.6f}".format("s" if S_pol else "p",m_x,m_y,Tmode))
if m_x == 0 and m_y == 0:
Tsum += Tmode
elif (m_x != 0 and m_y == 0) or (m_x == 0 and m_y != 0):
Tsum += 2*Tmode
else:
Tsum += 4*Tmode


r_flux = mp.get_fluxes(refl_flux)
t_flux = mp.get_fluxes(tran_flux)
Rflux = -r_flux[0]/input_flux[0]
Tflux = t_flux[0]/input_flux[0]

print("refl:, {}, {}".format(Rsum,Rflux))
print("tran:, {}, {}".format(Tsum,Tflux))
print("sum:, {}, {}".format(Rsum+Tsum,Rflux+Tflux))

self.assertAlmostEqual(Rsum,Rflux,places=2)
self.assertAlmostEqual(Tsum,Tflux,places=2)
self.assertAlmostEqual(Rsum+Tsum,1.00,places=2)


if __name__ == '__main__':
unittest.main()

0 comments on commit 06d418d

Please sign in to comment.