Computing gradients at any location within the simulation domain

[ianwilliamson]

Currently, Meep gradients may only be computed within a DesignRegion, where the material distribution is defined by a bilinear interpolation of a density from a numpy array.

However, it would be useful if the user could compute sensitivities / gradients at any location within the simulation domain (where the material may not necessarily be defined by a material grid). In principle, such a calculation should be possible since it only requires the forward and adjoint field distributions at each frequency. It would be great if a user could provide measurement point(s) that are not aligned with the Yee grid and Meep would take care of appropriately interpolating the gradients to those location(s).

As a hack, I attempted to achieve the above functionality by using the existing DesignRegion, filling the interpolated array with zeros, and overlaying other geometric objects:

However, Meep (correctly) returns a sensitivity map where the values overlapping with the other objects are zero:

[smartalecH]

Older versions of the code actually did this... as you say you can pull the forward and adjoint fields directly from the yee grid and interpolate wherever you need the gradient.

The advantage of the current infrastructure is that it allows multiple grids to interact (e.g. symmetries, rotations, constraints, etc.). The ugly part is backpropagating through the material grid interpolation scheme (especially when you have that much flexibility).

If you don't need that flexibility you can tweak the adjoint codebase to do what was done with the initial release of adjoint solver 2.0 (#1167).

Or if you are feeling extra ambitious, you can find a way to include both methods :-)

[ianwilliamson]

Thanks for the response Alec.

To handle symmetries, rotations, and other constraints on a grid parameterization, one could perform those operations entirely outside of Meep in a suitable automatic differentiation framework (autograd in Meep's case). In that situation you get the backpropagation through those operations for "free" and the user can compose their density array using any graph of differentiable operations. Meep would only need to receive a single array to interpolate onto the Yee grid. Here, it sounds like you're saying that your goal was to support the different cases of the grid_type parameter, which is fair. Or is there some advantage to baking those operations into Meep?

All of the above is tied to a particular material type (MaterialGrid). What I was suggesting is an interface for calculating d <some_quantity> / d epsilon_r at various points in the Meep simulation domain without regard to how the epsilon_r at those points was defined. My (possibly incorrect) understanding was that mp._get_gradient used by mpa.DesignRegion is intrinsically tied to the MaterialGrid.

you can tweak the adjoint codebase to do what was done with the initial release of adjoint solver 2.0

Are you referring to this line?

https://github.com/smartalecH/meep/blob/188ee252b60389b783d890b79ed986c61b855c62/python/adjoint/optimization_problem.py#L229

Is there any advantage / difference when using that approach instead of the call to mp._get_gradient (aside from handling the interpolation and details of different grid_type cases)?

https://github.com/NanoComp/meep/blob/master/python/adjoint/optimization_problem.py#L34

[smartalecH]

Glad to help.

To handle symmetries, rotations, and other constraints on a grid parameterization, one could perform those operations entirely outside of Meep in a suitable automatic differentiation framework (autograd in Meep's case).

This is true for basic structures, but not generalizable. Say you have a photonic crystal-esque structure with 6 fold symmetry (or some other topological insulator). To represent that you have to wrap your cartesian design grid outside of your rectangular simulation domain and back in on itself (which meep does automatically). The current setup ensures maximum flexibility.

My (possibly incorrect) understanding was that mp._get_gradient used by mpa.DesignRegion is intrinsically tied to the MaterialGrid.

You are completely correct. The _get_gradient() routine traverses through a tree of material grids at each point. If there are no material grids present, it doesn't do anything with the forward and adjoint fields.

Are you referring to this line?

Yes, that line calculates the gradient w.r.t. epsilon at each point. The original framework was meant to be agnostic to the design grid (as you point out) so that you could do whatever vector jacobian product was necessary to retrieve the gradient w.r.t to your actual design parameters (grid_type).

Even though the points are staggered on the yee grid, they are "brought together" with a weighted inner product from the definition of the adjoint. Things get hairy when there are off-diagonal components though (as I'm sure you know). All of this is implemented here.

You could actually generalize this routine to not check for material grids and simply run over every point in the cell. You could use a simpler version of the _get_gradient() swig wrapper too. It's probably not too much work, especially, since you are just cutting out features. The most annoying part is building the python API.

[ianwilliamson]

Thanks for the explanations.

The current setup ensures maximum flexibility.

Okay, I can't picture exactly what optimization you have in mind, but fair enough :-)

Yes, that line calculates the gradient w.r.t. epsilon at each point.

I'm probably just missing it, but where does the gradient of the system matrix enter into the computation on that line (https://github.com/smartalecH/meep/blob/188ee252b60389b783d890b79ed986c61b855c62/python/adjoint/optimization_problem.py#L229)? I notice you use finite differences to get it in get_material_gradient(), but I don't see the corresponding calculation in the older Python version.

[smartalecH]

I notice you use finite differences to get it in get_material_gradient(), but I don't see the corresponding calculation in the older Python version.

Yes, I didn't do a true "gradient of the system matrix" with the older version. Instead, I combined each vector component independently. That was a major drawback to the old method. The "flexibility" (laziness) made it hard to take that weighted inner product accurately.

Now, however, it shouldn't be a problem. Here's what I would do:

• Tell the OptimizationProblem to record the forward and adjoint fields across the whole domain if no design_regions are specified (currently it throws an error, I think).
• Slightly modify the current gradient calculation code (which only loops over design regions) to do something special if no design regions are detected.
• To do so, you can create a modified version of _get_gradient() that doesn't require the geometry tree -- you just pass the entire simulation volume and appropriate fields.
• This function will call material_grids_addgradient() which multiplies by the transpose of the interpolation matrix.
• The most important function for you is material_grids_addgradient_point() The current implementation assumes an arbitrary design parameter density. In contrast, you want to "backinterpolate" to the center of each grid voxel. But only when it detects MaterialGrid objects! You can add a flag (entire_grid) that skips the tree traversal. You also need to make sure the points are being interpolated in the proper locations.

There are lots of ways to get this to work, but I think this is the path of least resistance.

[ianwilliamson]

Thanks, I really appreciate the detailed response Alec!

I didn't do a true "gradient of the system matrix" with the older version.

I see. I am more familiar with this calculation under FDFD, which has an expression like the one shown below. The current Meep implementation (in C/C++) looks almost identical to this expression:

I'm not sure what the form of the system matrix is here, but it sounds like you're saying that the older approach introduces an error in the gradient, which makes sense.

There are lots of ways to get this to work, but I think this is the path of least resistance.

Thanks for outlining these steps... I think this seems reasonable. I don't think that I will have the time to contribute something for this in the near future, but it's good to know what would be required and I think I better understand the current setup / rationale.