Add support for DiffractedPlanewave to EigenmodeCoefficient objective function of adjoint solver
#2047
Comments
|
Can't you already pass arbitrary keyword arguments corresponding to |
|
Unfortunately, it's not just a simple matter of passing a |
|
Note that for the inverse design of diffractive structures, an alternative to using an objective function based on an Perhaps this means that it is not absolutely necessary to add this feature after all? |
|
Here is an attempt to verify the claim that it is possible to compute the diffracted orders of a 2D grating (a cylinder with axis along z positioned on a substrate with a unit cell that is periodic in the xy plane with normally incident planewave along z) using the near-to-far field transformation. The test involves using the Poynting theorem to verify that the sum of the transmittance in z of all the transmitted orders computed using mode decomposition is equivalent to the sum of the radial Poynting flux of the transmitted orders computed using the near-to-far field transformation. This demonstration is similar to Tutorial/Near to Far-field Spectra/Radiation Pattern of an Antenna which verified that the total flux in the far fields is independent of the actual shape of the bounding surface used to compute the flux. Unfortunately, there is currently a large difference in the results from the two approaches because of a likely bug. import meep as mp
import math
import numpy as np
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.4
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()
P_pol = True
# plane of incidence is XZ
# P/TM polarization: Ex, S/TE polarization: Ey
src_cmpt = mp.Ex if P_pol else mp.Ey
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)
tran_pt = mp.Vector3(0,0,0.5*sz-dpml)
tran_flux = sim.add_mode_monitor(fcen,
0,
1,
mp.ModeRegion(center=tran_pt,
size=mp.Vector3(sx,sy,0)))
stop_cond = mp.stop_when_fields_decayed(10,src_cmpt,tran_pt,1e-7)
sim.run(until_after_sources=stop_cond)
input_flux = mp.get_fluxes(tran_flux)
input_flux_data = sim.get_flux_data(tran_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)
nearfield = sim.add_near2far(fcen,
0,
1,
mp.Near2FarRegion(center=tran_pt,
size=mp.Vector3(sx,sy,0)),
nperiods=100)
tran_flux = sim.add_mode_monitor(fcen,
0,
1,
mp.ModeRegion(center=tran_pt,
size=mp.Vector3(sx,sy,0)))
sim.run(until_after_sources=stop_cond)
# length of far-field points
r = 1000*wvl
# sum the radial Poynting flux from the far fields of all transmitted orders
Pr_sum = 0
# 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
# determine the far-field point by scaling the wavevector of the diffracted mode
kmode = mp.Vector3(m_x/sx,m_y/sy,(fcen**2-(m_x/sx)**2-(m_y/sy)**2)**0.5)
rpt = kmode.unit().scale(r)
ff = sim.get_farfield(nearfield,rpt)
Efar = [np.conj(ff[j]) for j in range(3)]
Hfar = [ff[j+3] for j in range(3)]
Px = np.real(Efar[1]*Hfar[2]-Efar[2]*Hfar[1]) # Ey*Hz - Ez*Hy
Py = np.real(Efar[2]*Hfar[0]-Efar[0]*Hfar[2]) # Ez*Hx - Ex*Hz
Pz = np.real(Efar[0]*Hfar[1]-Efar[1]*Hfar[0]) # Ex*Hy - Ey*Hx
Pr = np.sqrt(np.square(Px)+np.square(Py)+np.square(Pz))
print("n2f:, {}, {}, {:.6f}".format(m_x,m_y,Pr))
if m_x == 0 and m_y == 0:
Pr_sum += Pr
elif (m_x != 0 and m_y == 0) or (m_x == 0 and m_y != 0):
Pr_sum += 2*Pr
else:
Pr_sum += 4*Pr
t_flux = mp.get_fluxes(tran_flux)
Tflux = t_flux[0]/input_flux[0]
# integrate the radial Poynting flux over a hemisphere
radial_flux = Pr_sum*2*np.pi*r*r/(nm_x*nm_y)
print("tran:, {:.6f}, {:.6f}, {:.6f}".format(Tsum,Tflux,radial_flux)) |
This won't work. The diffracted orders are planewaves, so they fill all space above the surface — you can't look at individual points in the far field and only see individual diffracted orders. |
Currently, the
EigenmodeCoefficientobjective function of the adjoint solver can only be specified using an eigenmode/band number (an integer) via itsmodemember parameter. It would be useful to extend this feature to also support aDiffractedPlanewaveobject similar to what already exists for theEigenmodeSource(via itseig_bandproperty) andget_eigenmode_coefficients(via itsmodeproeprty). This would make it easier to use the adjoint solver for the design of diffractive structures.The text was updated successfully, but these errors were encountered: