Add support for DiffractedPlanewave to EigenmodeCoefficient objective function of adjoint solver

[oskooi]

Currently, the EigenmodeCoefficient objective function of the adjoint solver can only be specified using an eigenmode/band number (an integer) via its mode member parameter. It would be useful to extend this feature to also support a DiffractedPlanewave object similar to what already exists for the EigenmodeSource (via its eig_band property) and get_eigenmode_coefficients (via its mode proeprty). This would make it easier to use the adjoint solver for the design of diffractive structures.

[smartalecH]

Can't you already pass arbitrary keyword arguments corresponding to get_eigenmode_coefficients to the EigenmodeCoefficient object?

[oskooi]

Unfortunately, it's not just a simple matter of passing a DiffractedPlanewave object as the mode parameter of the EigenmodeCoefficient class constructor to enable this feature because the __call__ and place_adjoint_source functions have been set up specifically to use a band number for defining an eigenmode.

[oskooi]

Note that for the inverse design of diffractive structures, an alternative to using an objective function based on an EigenmodeCoefficient object with a DiffractedPlanewave argument is to use a Near2FarFields object (which obviously already exists and has been debugged). These are functionally equivalent (i.e., for a given diffractive structure, you know exactly where to compute its diffracted orders in the far field). Also, the computational cost of each approach is similar.

Perhaps this means that it is not absolutely necessary to add this feature after all?

[oskooi]

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))

[stevengj]

Note that for the inverse design of diffractive structures, an alternative to using an objective function based on an EigenmodeCoefficient object with a DiffractedPlanewave argument is to use a Near2FarFields object (which obviously already exists and has been debugged). These are functionally equivalent (i.e., for a given diffractive structure, you know exactly where to compute its diffracted orders in the far field). Also, the computational cost of each approach is similar.

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.