Problems in designing Metagrating3D with Meep

[mawc2019]

The adjoint and finite-difference gradients are inconsistent when the EigenmodeCoefficient adjoint solver in Meep is applied to a case similar to that in @oskooi's metagrating-meep.py. The tests were based on the latest version of Meep, namely, v1.25.0. The code ran with 32 processes. The comparisons between adjoint and finite-difference gradients are as follows.

The code is as follows. It defines the same geometry as that in metagrating-meep.py, but does not include DiffractedPlanewave and the input flux.

import meep as mp
import meep.adjoint as mpa
import numpy as np
import autograd.numpy as npa
import math
from matplotlib import pyplot as plt

resolution = 50  # pixels/μm

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

theta_d = math.radians(50.0)  # deflection angle

wvl = 1.05  # wavelength
fcen = 1/wvl

px = wvl/math.sin(theta_d)  # period in x
py = 0.5*wvl  # period in y

dpml = 1.0  # PML thickness
gh = 0.325  # grating height
dsub = 5.0  # substrate thickness
dair = 5.0  # air padding

sz = dpml+dsub+gh+dair+dpml

cell_size = mp.Vector3(px,py,sz)

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


# plane of incidence is XZ
src_pt = mp.Vector3(0,0,-0.5*sz+dpml+0.5*dsub)
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=0.1*fcen),
                     size=mp.Vector3(px,py,0),
                     center=src_pt,
                     component=mp.Ex)]


design_region_resolution = 2*resolution
Nx, Ny, Nz = int(round(design_region_resolution * px)), int(round(design_region_resolution * py)), 1

design_variables = mp.MaterialGrid(mp.Vector3(Nx, Ny, Nz), mp.air, Si, do_averaging=False, grid_type="U_MEAN")
design_region = mpa.DesignRegion(design_variables,
                                 volume=mp.Volume(center=mp.Vector3(0,0,-0.5*sz+dpml+dsub+0.5*gh),
                                                  size=mp.Vector3(px,py,gh)))

geometry = [mp.Block(center=mp.Vector3(0,0,-0.5*sz+0.5*(dpml+dsub)),
                     size=mp.Vector3(mp.inf,mp.inf,dpml+dsub),
                     material=SiO2),
            mp.Block(center=design_region.center, 
                     size=design_region.size, 
                     material=design_variables)]

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

mode = 1
emc = mpa.EigenmodeCoefficient(sim,
                               mp.Volume(center=mp.Vector3(0,0,0.5*sz-dpml),size=mp.Vector3(px,py,0)),
                               mode)
ob_list = [emc]

def J(emc):
    return npa.abs(emc[0])**2

opt = mpa.OptimizationProblem(simulation=sim,objective_functions=J,objective_arguments=ob_list,
                              design_regions=[design_region],fcen=fcen,df=0,nf=1,decay_by=1e-12)

np.random.seed(20221119)
x0 = 0.5*np.ones(Nx*Ny)

f0, dJ_dn = opt([x0])

db,choose = 1e-4,5

g_discrete, idx = opt.calculate_fd_gradient(num_gradients=choose,db=db)
g_discrete = np.array(g_discrete).flatten()
g_adjoint = dJ_dn.flatten()

print("Randomly selected indices: ",idx)
print("Finite-difference gradient: ",g_discrete)
print("Adjoint gradient: ",g_adjoint[idx])

plt.plot(g_discrete,g_adjoint[idx],"ro")
plt.plot(g_discrete,g_discrete,"g-",label="y=x")
plt.xlabel("Finite-difference gradient")
plt.ylabel("Adjoint gradient")
plt.legend();

[oskooi]

Since this is a single-frequency problem and the accuracy of the adjoint gradients seems to degrade with increasing source bandwidth (as demonstrated in NanoComp/meep#2314), can you try redoing your test using a narrower bandwidth source, e.g. fwidth=0.02*fcen?

[mawc2019]

After replacing fwidth=0.1*fcen with fwidth=0.02*fcen in mp.Source, the adjoint and finite-difference gradients are still inconsistent. The comparisons are as follows.

In addition to v1.25.0, I also tried v1.23.0. The adjoint and finite-difference gradients are also inconsistent.

[mawc2019]

Apart from the issue above, the results may change when the number of processes changes. The above tests used 32 processes, but the results with Ex source become different when only 1 process is used. The comparison of the results between 1 process and 32 processes is as follows. The relative difference is calculated as 32-process result / 1-process result − 1.

On the other hand, the change in the results with Ey source is negligible.

The tests here also use fwidth=0.02*fcen instead of fwidth=0.1*fcen.

[mawc2019]

With random design parameters x0=np.random.rand(Nx*Ny), the adjoint and finite-difference gradients are still inconsistent. I also tried different resolution and/or design_region_resolution with x0=0.5*np.ones(Nx*Ny), but the adjoint and finite-difference gradients are inconsistent, either. So we probably cannot avoid the problem like here and here. All the tests here involve fwidth=0.02*fcen.
Any suggestions on the next step, @oskooi ?

[mawc2019]

Unlike waveguide_mode_converter where the adjoint and finite-difference gradients are inconsistent only in some special cases (such as uniform design patterns), in all tests here for Metagrating3D, the adjoint and finite-difference gradients are inconsistent.

However, to make the setup smaller, I did not add DiffractedPlanewave to EigenmodeCoefficient in the above tests. I can add DiffractedPlanewave to EigenmodeCoefficient to see whether the issue remains although it seems probably so. Can DiffractedPlanewave be readily added to the argument list of EigenmodeCoefficient, like
mpa.EigenmodeCoefficient(sim,mp.Volume(center=mp.Vector3(0,0,0.5*sz-dpml),size=mp.Vector3(px,py,0)),mode,mp.DiffractedPlanewave((1,0,0),mp.Vector3(1,0,0),0,1)), @oskooi ? It does not produce errors, but in the some existing examples that involve DiffractedPlanewave and get_eigenmode_coefficients, the mode does not need to specified.

[mawc2019]

Here is a different issue. I ran @oskooi's meep script metagrating-meep.py using the latest version of Meep. The calculated TM-polarization deflection efficiency for device1.csv is 88.8%, which is noticeably lower than the reported results in the README file, i.e., 95.7% (RETICOLO/RCWA) and 95.5% (MEEP/FDTD). Which version of Meep was used to obtain 95.5%?

[mawc2019]

The optimization is ongoing. The latest design pattern is as follows, which is not completely binarized. The resolution of FDTD simulation is 50 while the resolution of the design region is 100. The minimum length scale constraint is not imposed, but I set a minimum length of 50 nm to obtain a filter radius for the conic filter, which is 111.8 nm.

The evaluation history during optimization is as follows. The latest deflection efficiency is 93.6%, which is lower than 95.7% (RETICOLO/RCWA) but higher than 88.8% mentioned above.

The Meep script is as follows, which resembles @oskooi's script.

import meep as mp
import meep.adjoint as mpa
import numpy as np
import autograd.numpy as npa
from autograd import tensor_jacobian_product
import math
import nlopt

resolution = 50  # pixels/μm
design_region_resolution = 100
nSi = 3.45
Si = mp.Medium(index=nSi)
nSiO2 = 1.45
SiO2 = mp.Medium(index=nSiO2)

theta_d = math.radians(50.0)  # deflection angle

wvl = 1.05  # wavelength
fcen = 1/wvl

px = wvl/math.sin(theta_d)  # period in x
py = 0.5*wvl  # period in y

dpml = 1.0  # PML thickness
gh = 0.325  # grating height
dsub = 5.0  # substrate thickness
dair = 5.0  # air padding

sz = dpml+dsub+gh+dair+dpml

cell_size = mp.Vector3(px,py,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
src_cmpt = mp.Ex if P_pol else mp.Ey
src_pt = mp.Vector3(0,0,-0.5*sz+dpml+0.5*dsub)
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=0.1*fcen),
                     size=mp.Vector3(px,py,0),
                     center=src_pt,
                     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),
                                          size=mp.Vector3(px,py,0)))

stop_cond = mp.stop_when_fields_decayed(10,src_cmpt,src_pt,1e-7)
sim.run(until_after_sources=stop_cond)

input_flux = mp.get_fluxes(flux)

sim.reset_meep()

Nx, Ny, Nz = int(round(design_region_resolution * px)), int(round(design_region_resolution * py)), 1
design_variables = mp.MaterialGrid(mp.Vector3(Nx, Ny, Nz), mp.air, Si, do_averaging=False, grid_type="U_MEAN")
design_region = mpa.DesignRegion(design_variables,
                                 volume=mp.Volume(center=mp.Vector3(0,0,-0.5*sz+dpml+dsub+0.5*gh),
                                                  size=mp.Vector3(px,py,gh)))

geometry = [mp.Block(center=mp.Vector3(0,0,-0.5*sz+0.5*(dpml+dsub)),
                     size=mp.Vector3(mp.inf,mp.inf,dpml+dsub),
                     material=SiO2),
            mp.Block(center=design_region.center, 
                     size=design_region.size, 
                     material=design_variables)]

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

diffraction_order = 1
kdiff = mp.Vector3(diffraction_order/px,0,np.sqrt(fcen**2-(diffraction_order/px)**2))
emc = mpa.EigenmodeCoefficient(sim,mp.Volume(center=mp.Vector3(0,0,0.5*sz-dpml),size=mp.Vector3(px,py,0)),
                               mode=1,eig_parity = mp.EVEN_Y if P_pol else mp.ODD_Y,
                               kpoint_func=lambda *not_used: kdiff)
ob_list = [emc]

def J(emc):
    return npa.abs(emc[0])**2 / input_flux[0]

opt = mpa.OptimizationProblem(simulation=sim,objective_functions=J,objective_arguments=ob_list,
                              design_regions=[design_region],fcen=fcen,df=0,nf=1,decay_by=1e-6)

minimum_length = 0.05  # minimum length scale (microns)
eta_i = 0.5  # blueprint (or intermediate) design field thresholding point (between 0 and 1)
eta_e = 0.55  # erosion design field thresholding point (between 0 and 1)
eta_d = 1 - eta_e  # dilation design field thresholding point (between 0 and 1)
filter_radius = mpa.get_conic_radius_from_eta_e(minimum_length, eta_e)

def mapping(x, eta, beta):
    # filter
    filtered_field = mpa.conic_filter(x,filter_radius,px,py,design_region_resolution)
    # projection
    projected_field = mpa.tanh_projection(filtered_field, beta, eta)
    return projected_field.flatten()

evaluation_history = []
def f(v, gradient, beta):
    eff,dJ_dn = opt([mapping(v, eta_i, beta)])
    if gradient.size > 0:
        gradient[:] = tensor_jacobian_product(mapping, 0)(v, eta_i, beta, dJ_dn)
    evaluation_history.append(eff)
    np.savetxt('history.dat',evaluation_history)
    np.savetxt('pattern.dat',mapping(v, eta_i, beta).reshape(Nx,Ny))
    return eff

np.random.seed(20221119)
x = 0.5*np.ones(Nx*Ny*Nz)

lb, ub = np.zeros(Nx*Ny*Nz), np.ones(Nx*Ny*Nz) # lower and upper bounds

update_factor = 10000
algorithm = nlopt.LD_MMA

cur_beta = 4
beta_scale = 2
num_betas = 6
update_factor = 60
tol = 1e-8
for iters in range(num_betas):
    solver = nlopt.opt(algorithm, Nx*Ny*Nz)
    solver.set_lower_bounds(lb)
    solver.set_upper_bounds(ub)
    solver.set_max_objective(lambda a,g: f(a,g,cur_beta))
    solver.set_maxeval(update_factor)
    solver.set_xtol_rel(tol)
    x[:] = solver.optimize(x)
    cur_beta = cur_beta*beta_scale

print("Optimized indices: ",x)
np.savetxt('optimized_pattern.dat', x)

[oskooi]

For reference, here is the minimum lengthscale constraint function (borrowed from @mochen4's LDOS cavity optimization) and the optimization setup I used when testing the single-frequency waveguide mode converter while debugging #32. geom_c is the hyperparameter that requires some tuning.

Also, to ensure a binarized final design, you may also want to include damping in the MaterialGrid and turn on subpixel smoothing in order to use large β values.

def glc(x, gradient, eta, beta):
    """Constraint function for the minimum linewidth.                                                                     
                                                                                                                          
       Args:                                                                                                              
         x: 1d array of size Nx*Ny containing design weights.                                                             
         gradient: the Jacobian matrix with dimensions (Nx*Ny,                                                            
                   num. wavelengths) modified in place.                                                                   
         beta: bias parameter for projection.                                                                             
    """
    filter_f = lambda v: mpa.conic_filter(
        v.reshape(Nx,Ny),
        filter_radius,
        design_region_size.x,
        design_region_size.y,
        design_region_resolution,
    )
    threshold_f = lambda v: mpa.tanh_projection(v,beta,eta)
    geom_c = (filter_radius/design_region_resolution)*5000

    g_s = mpa.constraint_solid(
        x,
        geom_c,
        eta_e,
        filter_f,
        threshold_f,
        resolution,
    )

    g_s_grad = grad(mpa.constraint_solid,0)(
        x,
        geom_c,
        eta_e,
        filter_f,
        threshold_f,
        resolution,
    )

    gradient[:] = np.asarray(g_s_grad)

    return float(g_s)

    betas = [8, 16, 32, 64, 128, 256]
    max_evals = [70, 70, 70, 100, 100, 100]
    for beta, max_eval in zip(betas,max_evals):
        solver = nlopt.opt(algorithm, n)
        solver.set_lower_bounds(lb)
        solver.set_upper_bounds(ub)
        solver.set_min_objective(lambda x, g: f(x, g, eta_i, beta))
        solver.set_maxeval(max_eval)
        solver.set_param("dual_ftol_rel",1e-8)
        # apply the minimum linewidth constraint                                                                          
        # only in the final "epoch"                                                                                       
        tol_lw = 1e-5
        if beta == betas[-1]:
            solver.add_inequality_constraint(
                lambda x, g: glc(x, g, eta_i, beta),
                tol_lw,
            )

        x[:] = solver.optimize(x)

[mawc2019]

Thanks! However, minimum-length-scale constraints were not imposed when 3d metagratings were designed by RETICOLO/RCWA. Instead, only some filtering and projection were involved there. To make a fair comparison, perhaps we should not impose minimum-length-scale constraints either when performing optimization with Meep? Alternatively, should minimum-length-scale constraints be imposed on both RETICOLO/RCWA design and MEEP/FDTD design?

[oskooi]

Alternatively, should minimum-length-scale constraints be imposed on both RETICOLO/RCWA design and MEEP/FDTD design?

In principle, the Reticolo/RCWA designs should have been designed using minimum lengthscale constraints. For the three devices designed using Reticolo/RCWA and checked into this repository (device1.csv, device2.csv, device3.csv), the minimum lengthscales measured using your ruler are: 25 nm, 12 nm, and 37 nm. To obtain comparable designs using Meep, it might be useful to impose a lengthscale constraint.

[mawc2019]

For the three devices designed using Reticolo/RCWA and checked into this repository (device1.csv, device2.csv, device3.csv), the minimum lengthscales measured using your ruler are: 25 nm, 12 nm, and 37 nm. To obtain comparable designs using Meep, it might be useful to impose a lengthscale constraint.

Imposing constraints based on these minimum length scales may lead to additional issues. The details are as follows.

These minimum length scales arise from device1.csv, device2.csv, and device3.csv, but according to You Zhou's email, these design patterns were not obtained directly from optimization, but were interpolated from the optimized design patterns. The pixel size of the design region in optimization is about 12 nm, while the pixel size in the interpolated design patterns (device1.csv, device2.csv, and device3.csv) is about 3 nm. The minimum length scales you mentioned above are measured from the interpolated design patterns, not from the design patterns obtained directly from optimization. To make a fair comparison, we may need to use the same design-region resolution as that in their optimization, which corresponds to the pixel size of 12 nm. Consequently, the measured minimum length scales, i.e., 25 nm, 12 nm, and 37 nm, amount to the sizes of only 2, 1, and 3 pixels. Imposing constraints based on such small length scales may be problematic. Even if those minimum length scales are successfully imposed, checking them using ruler.py (whatever version) may also be problematic because of the error at 1-pixel level.
In addition, as shown in #28, the minimum length scales in device3.csv varies noticeably among different versions of ruler.py.

[mawc2019]

The calculated TM-polarization deflection efficiency for device1.csv is 88.8%, which is noticeably lower than the reported results in the README file, i.e., 95.7% (RETICOLO/RCWA) and 95.5% (MEEP/FDTD). Which version of Meep was used to obtain 95.5%?

When the resolution is 200, the TM-polarization deflection efficiency calculated by Meep for device1.csv becomes the same as that reported in the README file, i.e., 95.5%.

The optimization is ongoing.

After more iterations, the structure became nearly binarized. The optimized deflection efficiency is 90% when the resolution is 50. However, when the resolution is 200, the deflection efficiency for the same structure is only 45%, which may imply a frequency shift, i.e., the working wavelength of the device is different from the target wavelength. The spectrum is as follows, where the dashed line indicates the target frequency, i.e., 1/1/05=0.952.

[mawc2019]

In principle, the Reticolo/RCWA designs should have been designed using minimum lengthscale constraints. For the three devices designed using Reticolo/RCWA and checked into this repository (device1.csv, device2.csv, device3.csv), the minimum lengthscales measured using your ruler are: 25 nm, 12 nm, and 37 nm. To obtain comparable designs using Meep, it might be useful to impose a lengthscale constraint.

These results are based on interpolated design patterns, which are not raw design patterns obtained directly from optimization. The minimum length scales of the raw design patterns are different. As listed here, the minimum length scales of the raw design patterns are 58 to 74 nm. These values are not at the single-pixel level because the size of a pixel is 11.6 to 11.7 nm.
Therefore, when we design the structure with Meep, maybe we could impose a few minimum length scales such as 55, 60, 65, 70, and 75 nm.