Skip to content

Commit

Permalink
unit test for mode decomposition in 3d (#2032)
Browse files Browse the repository at this point in the history 
* unit test for mode decomposition in 3d

* increase check tolerance for reflectance
  • Loading branch information
@oskooi
oskooi authored Apr 12, 2022
1 parent 9c5582e commit 75ced52
Showing 1 changed file with 195 additions and 14 deletions.
209 changes: 195 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,171 @@ 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 = 1.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)]

symmetries = [mp.Mirror(direction=mp.X,phase=-1),
mp.Mirror(direction=mp.Y,phase=+1)]

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

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

stop_cond = mp.stop_when_energy_decayed(20,1e-6)
sim.run(until_after_sources=stop_cond)

input_flux = mp.get_fluxes(refl_flux)
input_flux_data = sim.get_flux_data(refl_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,
symmetries=symmetries)

refl_flux = sim.add_mode_monitor(fcen,
0,
1,
mp.ModeRegion(center=refl_pt,
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))

## to obtain agreement for two decimal digits,
## the resolution must be increased to 200
self.assertAlmostEqual(Rsum,Rflux,places=1)
self.assertAlmostEqual(Tsum,Tflux,places=2)
self.assertAlmostEqual(Rsum+Tsum,1.00,places=1)


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

0 comments on commit 75ced52

Please sign in to comment.