Meep topology optimization of waveguide mode converter

[oskooi]

As a direct comparison to the two designs for the waveguide mode converter from Schubert et al., 2022, here are two designs obtained using topology optimization via Meep. The multi-objective function involves minimizing the worst case of $R+(1-T)$ where $R$ is $|S_{11}|^2$ for mode 1 and $T$ is $|S_{21}|^2$ for mode 2 across six different wavelengths. The minimum linewidth criteria is 90 nm. The optimization uses the method of moving asymptotes (MMA) algorithm from NLopt.

There are two sets of results: (1) no mirror symmetry and (2) mirror symmetry along the $x$ axis. A plot of the reflectance and transmittance spectra for each case demonstrates good performance across all wavelengths although (2) is better. It is possible that the results can be improved by e.g. modifying the objective function to be $\min\max(R,1-T)$ as well as the usual adjustment of the hyperparameters.

This PR includes two CSV files of the optimized designs and the Meep script used to generate these results.

Requires NanoComp/meep#2271.

cc @ianwilliamson, @smartalecH.

[stevengj]

Looks good!

The fact that the symmetric structure is better is probably just random luck finding a better local optimum?

We'll want to eventually impose explicit lengthscale constraints (a filter radius is not sufficient) Note, for example, the sharp corner present at the edges of the design region.

[smartalecH]

Cool!

We'll want to eventually impose explicit lengthscale constraints (a filter radius is not sufficient) Note, for example, the sharp corner present at the edges of the design region.

Unfortunately, a lengthscale constraint won't fix the violations along the edge of the domain. Recall that we pad the domain (quite a bit) to perform the filtering. This was intentional, so that the optimizer wouldn't try to round corners meant to connect with waveguides (where it has no awareness).

Instead, the easiest way to deal with these "corner cases" (oof) is to force the boundary of the domain to either be solid or void, like we do in the tutorials. Then the filters can properly pad. We typically use a simple bitmask, along with the jax/autograd function where() to enforce the boundaries (you could do it with the optimizer itself too, but this seems cleaner).

[smartalecH]

cc @mfschubert

[smartalecH]

@oskooi FYI design constraints look something like:

def glc(result,x,gradient,beta,opt):
    t = x[0] # dummy parameter
    v = x[1:] # design parameters
    a1 = 1e-6 # hyper parameter!! (this one is rather important)
    b1 = 0 # hyper parameter!! (this one is less important... maybe keep it as is...)
    gradient[:,0] = -a1
    
    N = opt.design_regions[0].num_design_params
    filter_f = lambda a: mpa.conic_filter(a.reshape(Nx,Ny),filter_radius, 
        design_region_size.x,
        design_region_size.y,
        design_region_resolution)
    threshold_f = lambda a: mpa.tanh_projection(a,beta,0.5)
    eta_e = 0.75
    eta_d = 0.25
    c0 = (filter_radius*1/resolution)**4

    # you can plot these for debugging...
    I1 = mpa.indicator_solid(v,c0,filter_f,threshold_f,1)
    I2 = mpa.indicator_void(v,c0,filter_f,threshold_f,1)

    M1 = lambda a: mpa.constraint_solid(a,c0,eta_e,filter_f,threshold_f,1)
    M2 = lambda a: mpa.constraint_void(a,c0,eta_d,filter_f,threshold_f,1)

    g1 = grad(M1)(v)
    g2 = grad(M2)(v)
    
    result[0] = M1(v) - a1*t - b1
    result[1] = M2(v) - a1*t - b1
        
    gradient[0,1:] = g1.flatten()
    gradient[1,1:] = g2.flatten()

and added to the optimizer with

    solver.add_inequality_mconstraint(lambda r,x,g: glc(r,x,g,beta[iters],opt),[1e-8]*2)

(using simplified code to build results from our early paper... the above is untested but gives the gist)

[smartalecH]

(also make sure to only add your glc constraint during your last $\beta$...)

[oskooi]

@smartalecH: thanks for the code snippet to implement the linewidth constraint. I modified the script in mode_converter_meep_opt.py in this repo to include this constraint and also, following your earlier suggestion, imposed an additional constraint of a one-pixel thick padding around the edges of the design region to prevent sharp corners from arising there. The modified script is available in this gist.

Here are a couple of sample designs generated using these constraints. The performance, particularly the transmittance spectrum, is noticeably worse compared to the unconstrained design. I think this is because the optimization is overly constrained as the final designs are too smooth and do not contain any small features. I tried adjusting the a1 parameter in the glc function (1e-5, 1e-6) but this did not change the results much.

What can we do to relax the constraints a little and enable the optimizer to explore a larger design space?

(I also tried changing the objective function to $\min\max(R,1-T)$ but this made the final results significantly worse. Thus, for the time being, I am still using $\min\max(R+1-T)$.)

[smartalecH]

@oskooi can you include the convergence plots?

This is where the "art" kicks in (not as much "science", unfortunately).

Suggestion: typically, when optimizing with glc constraints, I typically only run for three $\beta$ 's: $\beta=8,16,32$. I'll run each $\beta$ for about 30-70 iterations (depending on the device). As discussed earlier, only apply the GLC constraints during the last epoch (i.e. the last $\beta$).

[oskooi]

can you include the convergence plots?

From the convergence plots for these two designhs, it appears that the performance got significantly worse towards the end of the optimization right when the linewidth constraint was activated. I should have probably run for a few more iterations to find a better design.

With regards to the hyperparameters, what about eta_e and eta_d in your script? Should the values 0.75 and 0.25, respectively, also be tuned?

[smartalecH]

In the future I would plot it on a log scale so you can better track convergence. But you're right, you certainly need to run that last Beta longer. Furthermore, you should run fewer betas, each of them longer, so that each epoch converges to some "topology".

With regards to the hyperparameters, what about eta_e and eta_d in your script?

You could play with those, but I try to keep them at those default values. They correspond to a radius that is exactly the lengthscale constraint (which is somewhat elegant). By changing them, you effectively allow for smaller radii. But that could produce topologies in the earlier epochs that struggle satisfying the constraint. Conversely, you could tune them such that the radius is larger, but then you'll struggle getting good performance (typically). These default values seem to be a happy middle ground.