Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

unit test for mode decomposition in 3d #2032

Merged
merged 2 commits into from
Apr 12, 2022
Merged

unit test for mode decomposition in 3d #2032

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
fac8638
unit test for mode decomposition in 3d
oskooi Apr 5, 2022
de9c776
increase check tolerance for reflectance
oskooi Apr 12, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
209 changes: 195 additions & 14 deletions python/tests/test_mode_decomposition.py
View file
Open in desktop
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()