Tutorial example for computing the adjoint gradient of a level set

[oskooi]

Closes #2235.

In order to obtain early feedback, this is the first commit which currently just includes two explanatory figures and an image showing simulation output along with the simulation script in what will eventually become a new tutorial for the adjoint solver. Text and additional results are forthcoming.

[stevengj]

If you are willing to do finite difference, in principle you should be able to skip the subpixel smoothing entirely. You can instead do:

1. Construct your level set function φ at a very high resolution for a given width h. Call this φ(h)
2. Tell numpy to downsample (with bilinear interpolation) to a lower resolution, comparable to the Meep resolution. Call this ρ(φ(h)). This will "anti-alias" the boundaries. This is the material grid you input to Meep.
3. Meep gives you the derivative with respect to ρ of your objective function. To get the derivative with respect to h, you just need a custom vJp that computes dρ/dh, e.g. by finite differences.

(You need high resolution for φ for finite differences to be accurate, since you can only change h by steps of one pixel in φ, unless you are willing to do your own subpixel averaging.)

[oskooi]

1. Construct your level set function φ at a very high resolution for a given width h. Call this φ(h)

Unfortunately, even using this new approach (downsampling from a grid for the level set $\phi(h)$ that that has 10X the resolution of the Yee grid) produces similar results to what I had previously: good agreement in the directional derivative for the adjoint gradient and finite difference for the $\mathcal{S}$ but not the $\mathcal{P}$ polarization.

The test uses a resolution of 400 pixels/um which is fairly high. It is not obvious why the results for the $\mathcal{P}$ polarization are inaccurate.

1. $\mathcal{S}$ polarization

directional-deriv:, -0.00096108 (finite difference), -0.00090948 (adjoint), 0.053690 (error)

2. $\mathcal{P}$ polarization

directional-deriv:, -0.00015576 (finite difference), -0.00112504 (adjoint), 6.222962 (error)

Script used to generate these results is in this gist.

[stevengj]

When you carry out the chain rule, you multiply the Meep d/dρ with the level-set dρ/dh derivative. The latter is only nonzero on the boundary pixels. So, if the Meep derivative is inaccurate right on the boundary, that would be exacerbated by this d/dh derivative.

My hypothesis at the moment is that the p-polarization Meep derivatives are less accurate at the boundary (due to the the differentiate–then–discretize errors), since the p polarization has a discontinuity at the boundary, compared to the s polarization (which is continuous).

To check this, try checking the Meep d/dρ derivatives against finite differences, but instead of using a random perturbation to ρ, do a perturbation right at the the boundary pixels.

[oskooi]

My hypothesis at the moment is that the p-polarization Meep derivatives are less accurate at the boundary (due to the the differentiate–then–discretize errors), since the p polarization has a discontinuity at the boundary, compared to the s polarization (which is continuous).

To check this, try checking the Meep d/dρ derivatives against finite differences, but instead of using a random perturbation to ρ, do a perturbation right at the the boundary pixels.

The P polarization Meep derivatives are indeed less accurate at the boundaries compared to the S polarization. To verify this, I used the example of the usual 1D binary grating and added a 2-pixel perturbation for the height (see gist). Results are similar for a 1-pixel height perturbation.

As shown in the results below, the error in the directional derivative is small and decreases with resolution for the S polarization but not the P polarization.

We can explain this feature in this tutorial. For now, the tutorial is restricted to just the S polarization.

$\mathcal{S}$ polarization

RESOLUTION_UM = 50
dir-deriv:, -2.331625e-07 (finite difference), -2.371818e-07 (adjoint), 0.01723843268351378 (error)

RESOLUTION_UM = 100
dir-deriv:, -1.212851e-07 (finite difference), -1.212852e-07 (adjoint), 1.295385302514402e-06 (error)

$\mathcal{P}$ polarization

RESOLUTION_UM = 50
dir-deriv:, -1.200119e-07 (finite difference), -1.417497e-07 (adjoint), 0.18113050352310484 (error)

RESOLUTION_UM = 100
dir-deriv:, 1.785888e-08 (finite difference), 1.255981e-08 (adjoint), 0.2967190280473706 (error)

[oskooi]

The test failures are due to Autograd's dependence on deprecated Numpy features as reported (and fixed but not yet merged) in HIPS/autograd#618.

[stevengj]

Would be good to check the in-plane polarization for the $\beta = \infty$ subpixel smoothing as well.

Maybe there is a bug in what @smartalecH calls the "recombination" step (the overlap between adjoint and forward fields) that is lowering the accuracy.

[oskooi]

@stevengj, this is ready to be merged pending final review.