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.

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

Problems in designing Metagrating3D with Meep #31

Open
mawc2019 opened this issue Nov 21, 2022 · 25 comments
Open

Problems in designing Metagrating3D with Meep #31

mawc2019 opened this issue Nov 21, 2022 · 25 comments

Comments

@mawc2019
Copy link
Collaborator

mawc2019 commented Nov 21, 2022

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.
image

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
Copy link
Collaborator

oskooi commented Nov 22, 2022

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
Copy link
Collaborator Author

mawc2019 commented Nov 22, 2022

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.
image

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

@mawc2019
Copy link
Collaborator Author

mawc2019 commented Nov 23, 2022

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.
image

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

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

@mawc2019
Copy link
Collaborator Author

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
Copy link
Collaborator Author

mawc2019 commented Dec 6, 2022

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
Copy link
Collaborator Author

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
Copy link
Collaborator Author

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.
image
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.
image

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
Copy link
Collaborator

oskooi commented Dec 13, 2022

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
Copy link
Collaborator Author

mawc2019 commented Dec 13, 2022

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
Copy link
Collaborator

oskooi commented Dec 15, 2022

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
Copy link
Collaborator Author

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
Copy link
Collaborator Author

mawc2019 commented Dec 20, 2022

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.
image

@mawc2019
Copy link
Collaborator Author

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.

@mawc2019
Copy link
Collaborator Author

mawc2019 commented Dec 23, 2022

The evaluation history of the objective during optimization is as follows. Before the 280th iteration, there are do_averaging = False and damping = 0 in MaterialGrid. At and after the 280th iteration, there are do_averaging = True and damping = 0.05*2*np.pi*fcen. The change of beta happens at the 70th, 140th, 210th, 280th, and 380th iterations.
image
The optimization here attempts to maximize the objective function, but after the 318th iteration, the objective function decreases significantly and never recovers. The objective at the 318th iteration is marked by the red dot in the above figure. The design pattern and the adjoint gradient at this step are as follows.
image
The value of rho given by nlopt is about 0.4126 or 0.04126 around this step.

@oskooi
Copy link
Collaborator

oskooi commented Jan 3, 2023

As discussed, when using the damping feature you should probably turn off subpixel smoothing since there could be a bad interaction between the two when the design is greyscale (i.e., the weights are intermediate between 0 and 1) which is degrading the accuracy of the adjoint gradient and thus spoiling the optimization.

Note that the left figure above does indicate the presence of intermediate weights which is where the adjoint gradient is largest in the right figure. Ideally, the adjoint gradient should be non zero only along the discontinuous contours.

@mawc2019
Copy link
Collaborator Author

mawc2019 commented Jan 10, 2023

when using the damping feature you should probably turn off subpixel smoothing since there could be a bad interaction between the two when the design is greyscale

I used damping=0.05*2*np.pi*fcen and do_averaging=False in the first few epochs, and later damping=0 and do_averaging=True. Meanwhile, I made the increase of beta more gradually. The history of the objective during optimization and the optimized design pattern are as follows.
image
The optimized deflection efficiency is 96.1%, which is slightly higher than the RCWA-optimized value. The minimum length scale estimated by our code is 73.8 nm, which is comparable to the minimum length scales of the RCWA-optimized design patterns. However, the optimized design pattern here is not completely binarized, even with beta=128 in the last epoch.

To binarize the design pattern, I tried to impose length-scale constraints, but all attempts so far lead to degradation of the objective function.
Then I abandoned the length-scale constraints, but instead, added a penalty term to the deflection efficiency, which resulted in a new objective function. The penalty term, which is proportional to np.sum((mapping(x, eta_i, beta)-0.5)**2), penalizes non-binarized values. The coefficient of this term needs to be tested and selected in order to effectively penalize non-binary values without ruining the deflection efficiency. The design pattern became more binarized using this method, as shown below. The estimated minimum length scale is 98 nm. However, when two structures above are enforced to become binarized, the deflection efficiencies are both only 89.3% under resolution=200.
image
I also tried to impose length-scale constraints, and meanwhile, impose another constraint to prevent the objective from decreasing too much. However, the design pattern did not become more binarized after more than a hundred iterations.

@mawc2019
Copy link
Collaborator Author

In constraint_solid(x, c, eta_e, filter_f, threshold_f, resolution) and constraint_void(x, c, eta_e, filter_f, threshold_f, resolution), I have been using a parameter c proportional to (filter_radius/design_region_resolution)**4 (with a coefficient to be determined via attempts), but in @mochen4's code for LDOS cavity optimization shown here, the coefficient is proportional to (filter_radius/design_region_resolution). Does the power of this expression also need to be determined via attempts?

@mochen4
Copy link
Collaborator

mochen4 commented Jan 22, 2023

In constraint_solid(x, c, eta_e, filter_f, threshold_f, resolution) and constraint_void(x, c, eta_e, filter_f, threshold_f, resolution), I have been using a parameter c proportional to (filter_radius/design_region_resolution)**4 (with a coefficient to be determined via attempts), but in @mochen4's code for LDOS cavity optimization shown here, the coefficient is proportional to (filter_radius/design_region_resolution). Does the power of this expression also need to be determined via attempts?

In my experience, as a hyper-parameter, c only needs to be at the right order of magnitude. So in practice it is equivalent. For the problems I worked on, I found the order is closer to (filter_radius/design_region_resolution) rather than (filter_radius/design_region_resolution)**4.

@mawc2019
Copy link
Collaborator Author

mawc2019 commented Jan 24, 2023

I found the order is closer to (filter_radius/design_region_resolution) rather than (filter_radius/design_region_resolution)**4.

I also tried this and found that the value of constraint_solid could be minimized when the constraint was imposed, in contrast to the previous cases where this value almost remained the same.
However, the value of the objective function still deteriorates significantly even if I use the new values of c and the pattern before imposing the constraint is largely binarized, such as the pattern below. The final beta here is only 64.
image

@mawc2019
Copy link
Collaborator Author

When I set beta=float("inf") and eps_averaging=True, I got RuntimeError: meep: simulation fields are NaN or Inf. Did I miss something when using this level-set method, @mochen4 ?

@mochen4
Copy link
Collaborator

mochen4 commented Jan 27, 2023

When I set beta=float("inf") and eps_averaging=True, I got RuntimeError: meep: simulation fields are NaN or Inf. Did I miss something when using this level-set method, @mochen4 ?

You should use the beta (projection) from the material grid (design variable); you should change beta there, something like design_variables.beta = .... Meanwhile, the mapping function should only have the filter but not the projection.

@mawc2019
Copy link
Collaborator Author

mawc2019 commented Jan 29, 2023

You should use the beta (projection) from the material grid (design variable); you should change beta there, something like design_variables.beta = ...

Sorry for my misunderstanding. I thought this only applied to the last epoch where beta was changed from a finite number to float("inf").

Now using this format with betas = [4, 8, 16, 32, 64, 96, 128, float("inf")], a binarized design pattern can be obtained. The evaluation history and the optimized design pattern are as follows.
image
The optimized diffraction efficiency is 93.6% when computed by the forward run of the adjoint solver, and 95.7% when computed via routine simulation without the adjoint solver. In both cases, the fdtd resolution is 50. The difference in diffraction efficiency may be due to different decay thresholds.

When the resolution is 200 in fdtd simulation, the peak diffraction efficiency becomes 95.9% and the peak frequency shifts to 0.964, which is 1.2% higher than the target frequency, i.e., 1/1.05=0.952. The diffraction spectrum is as follows.
image

The minimum length scale estimated by ruler.py is 73.8 nm. However, the design pattern clearly shows several sharp corners, implying that the minimum length scale should be near 0. The reason is related to a defect in ruler.py. Briefly speaking, those sharp corners have spurious minimum length scales, which can be quite large and sometimes masked by other smaller minimum length scales.

@mawc2019
Copy link
Collaborator Author

mawc2019 commented Jan 31, 2023

The diffraction efficiencies reported in README.md of Metagrating3D are not the maximum efficiencies around the target frequency. For example, the spectrum of diffraction efficiency of Device2 is as follows.
image
The peak diffraction efficiency is 95.5%, which is slightly higher than the reported ones, i.e., 93.3% (RETICOLO/RCWA) and 93.6% (MEEP/FDTD, but my test gives 93.4%). The peak frequency is 0.929, which is 2.5% lower than the target frequency, i.e., 1/1.05=0.952.
By the way, the 12nm hole in the design pattern of Device2 is unimportant. Filling that hole with solid material does not change the spectrum of diffraction efficiency obviously.

The spectrum of diffraction efficiency of Device3 is as follows.
image
The peak diffraction efficiency is 96.9%, which is slightly higher than the reported ones, i.e., 96.6% (RETICOLO/RCWA) and 95.0% (MEEP/FDTD). The peak frequency is 0.948, which is 0.5% lower than the target frequency, i.e., 1/1.05=0.952.

@mawc2019
Copy link
Collaborator Author

mawc2019 commented Feb 2, 2023

The sharp rectangular corners in the design pattern may result from rectangular shapes of filter kernels in Meep, as mentioned here. I am not sure whether the rectangular shapes have adverse effects on topology optimization, especially the binarization process and length-scale constraint, but I should probably change the filter kernels back to circular shapes and redo the optimization.

@mawc2019
Copy link
Collaborator Author

mawc2019 commented Feb 7, 2023

Here I use kernels with circular shapes and betas = [4, 8, 16, 32, 64, 96, 128, 256, float("inf")]. Starting from the same initial guess, i.e., x = 0.5*np.ones(Nx*Ny*Nz), the evaluation history and the optimized design pattern are as follows. The corners are not as sharp as those shown above. The optimized diffraction efficiency at the target frequency is 93.4%, which is calculated in the forward run of the adjoint solver.
image

Starting from a random guess, the evaluation history and the optimized design pattern are as follows. The optimized diffraction efficiency at the target frequency is 96.3%, which is calculated in the forward run of the adjoint solver.
image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants