[2/2] Fix MaterialGrid design-variable attribution in the adjoint gradient (adjacent grids; #1984 symmetries) - #3277
Conversation
a25b915 to
beef257
Compare
0d12f05 to
f63909f
Compare
|
I just tested this pr on a two-layer vertical grating coupler with partial etch layer calculation and it dramatically lowered the adjoint error along a directional derivative to 0.2% compared to finite differences. In #3274 the relative error was around 26% and in the latest release the relative error was around 300%. I am a bit surprised that the error dropped so much given that the discretization error for the forward calculation is on the order of 20% at the resolution I am using. Nonetheless, I am looking forward to trying out this pr in an optimization |
smartalecH
changed the title
|
This is really nice work, especially the symmetry cases. I ran into a similar, though distinct, issue with symmetries while validating #3274, and it was very difficult to diagnose. I really appreciate the rotated and mirrored regression coverage here. |
e272d6f to
7f6f73f
Compare
smartalecH
changed the title
|
Glad to hear it, thanks @lxvm!
Which "resolution" is this? The simulation resolution, or the "resolution" of your directional derivative (e.g. your Δx)? Note that meep's hybrid time-/frequency-domain adjoint should match finite difference gradients regardless of the discretization. That's because we are differentiating the discretized operator, not the continuous operator (e.g. differentiate then optimize vs optimize then differentiate). Of course, the step of your finite difference matters (as does whether it's one sided vs centered). But you should get pretty good agreement whether your resolution is 10 or 1000. We haven't always, and that's because we have bugs 🙃 but let's see if we can track those down now that we have some additional support from coding tools. |
|
Thanks @smartalecH ! The resolution I was referring to was the simulation resolution (i.e. 30 for my simulation for a coupler based on 500nm x 220 nm waveguide dimensions). Thanks for your comment because I had thought meep was using a differentiate then discretize approach, which is not the case. Glad to hear that I should expect good agreement of the adjoint gradients with finite differences (in my test I was using a finite difference step size of 10^-5). |
|
@lxvm did you try running an optimization after rebasing on this PR? |
|
Yes, I ran an optimization with this pr and it gave improved results compared to the prs leading up to this one. I've been able to do robust optimization of a single mode coupler with length scale and partial etch constraints that attained 50% efficiency which is the best result I've had so far |
|
This looks reasonable to me. |
A design region whose MaterialGrid touches another MaterialGrid gets a wrong
adjoint gradient, and, when the neighbour's grid is larger, an out-of-bounds
write.
material_grids_addgradient_point() resolves the material grid by a spatial
lookup at the point (geom_tree_search) and then interpolates into the supplied
v[] using *that* grid's grid_size and box coordinates. It never checks that the
grid it found is the one the design region being processed actually owns; the
comment above the loop says as much, and leaves it to the user to "only have one
unique design grid in this volume".
That is not something the user can arrange. add_dft_fields() snaps the design
region's monitor outward to enclosing grid nodes, so two abutting grids give one
region a monitor that covers a row of nodes lying inside the other grid's
object. Instrumented on a 2d case with two 21x11-node regions sharing a
boundary, the lower region's gradient pass touched 210 nodes of its own grid and
21 of the neighbour's -- exactly the shared row. Their contributions are
interpolated in the neighbour's box coordinates and summed into this region's
v[]: wrong values when the two grids are the same size, and a write past the end
of v[] when the neighbour's grid is larger (reproducibly aborts in malloc).
Fix by giving every Python MaterialGrid a stable id, carrying it into
material_data, and passing the design region's id down so a point that resolves
to a different grid is skipped. The id is per Python object, not per
material_data, so the documented case of one MaterialGrid backing several
geometric objects (symmetries) still accumulates from all of them. Nothing is
lost by skipping: the neighbouring design region visits those same nodes on its
own pass. owner_grid_id < 0 disables the filter.
Measured, adjoint / directional finite difference:
before after
2d, one design region 1.0000 1.0000
2d, two regions, 1 px apart 1.0000 1.0000
2d, two regions abutting 1.3508 1.0001
2d, abutting, interface off-node 0.3360 1.0000
3d, two 220 nm sublayers, lower 0.33-1.31 1.0000-1.0002
3d, two 220 nm sublayers, upper 2.35-5.46 1.0000
Not a chunk effect: num_chunks is 1 in the 2d cases above. Not subpixel
smoothing either -- eps_averaging=False makes it worse, since the mixing is in
the grid lookup rather than in the averaging. A design region abutting a static
(non-MaterialGrid) block is unaffected, which is what isolates this to two
adjacent grids.
The regression test covers the adjacent case twice: with equal grid sizes, where
a regression is a wrong value and fails as an assertion, and with unequal sizes,
where it is also the out-of-bounds write. Both run in serial in ~10 s.
test_adjoint_solver 13/13, test_adjoint_chunks + test_adjoint_utils +
test_subpixel_3d 12/12.
Note: material_grids_addgradient() gains a required owner_grid_id parameter
before du. The only in-tree caller is _get_gradient in python/meep.i.
Fixes NanoComp#1984. A design constrained to a symmetry by overlapping one MaterialGrid with transformed copies of itself got an adjoint gradient that was not the gradient -- wrong sign, not merely wrong scale. material_grids_addgradient_point() walked the overlapping copies and differentiated each in turn, but add_interpolate_weights() was handed the *first* copy's weights array, hoisted before the loop, while the stencil index came from the current copy's box coordinates. For a transformed copy those disagree: the finite difference in get_material_gradient() perturbs a sample that takes no part in that copy's interpolation. Untransformed overlaps survived by cancellation -- the indices coincide, so the second copy just repeats the first, and `scalegrad /= matgrid_val_count` divides the duplicate back out. That divide was itself a double count: get_material_gradient() finite-differences through matgrid_val(), so the 1/N of a U_MEAN average is already in the result. Restructured so the unit of differentiation is the design variable rather than the copy. Collect the copies of the design region's own grid that cover the point, take the union of their interpolation stencils, and for each variable in it perturb that entry in *every* copy at once. Copies are recognised by the grid_id added in the parent commit, which is what makes "the same design variable" expressible at this level. `scalegrad /= matgrid_val_count` is gone. This subsumes the parent commit's U_DEFAULT-scoped guard: filtering at collection time is correct for every grid_type, so the scoping is dropped. Measured, adjoint / directional finite difference: before after overlapped with a 90-degree rotation -0.1598 1.0003 overlapped with a mirror 0.2433 1.0000 overlapped, untransformed 1.0000 1.0000 overlapped, distinct grids 1.0000 1.0000 two adjacent design regions 1.3508 0.9986 single design region 1.0000 1.0000 U_MIN and U_PROD keep their existing value-dependent adjustments, applied to the single copy the caller allows for those kinds; overlapping grids of those kinds are still rejected upstream. Not fixed here, and pre-existing: an overlap whose second grid is left at all zeros reads 1.0077. It is step-independent, so not finite-difference error, and it moves neither with this change nor with the parent commit. test_adjoint_solver, test_adjoint_chunks, test_adjoint_utils, test_adjoint_cyl, test_adjoint_jax, test_adjoint_adjacent_grids 37/37; test_material_grid and test_subpixel_3d 8/8.
c5d7cfd to
15e72c0
Compare
|
Failing test? |
| do { | ||
| material_data *mg_cur = (material_data *)tpc->objects[oic].o->material; | ||
| if (owner_grid_id < 0 || mg_cur->grid_id == owner_grid_id) { | ||
| if (udatas.empty()) sz = mg_cur->grid_size; |
There was a problem hiding this comment.
cc @stevengj claude picked this up during the rebase (another corner case). It added a test (above) which had a tolerance slightly too tight for single precision.
4 participants
Fixes #1984.
Stacked on #3274 — based on
fix/adjoint-chunk-pairing-minimal, so the diff here is just these two commits. Merge #3274 first, then this retargets tomastercleanly.While working on a new 3D grating coupler TO example for meep (2 design regions) I noticed that the gradients were rather inaccurate. Claude helped me trace things down to the unique case where we have two adjacent design regions (e.g. a partial etched layer). There were different issues depending on the sizes of the design regions (e.g. if one design region was bigger than the other, then we got a segfault!)
The core issue was that the logic we had for handling multiple design grids during the recombination step wasn't properly bookkeeping the effects near the boundaries. This was because we weren't keeping track of "which grid was which." In addition, we had some hairy "copies" of grids we were trying to use to keep track of multiple grids (hence the crash). So the resulting gradient near the boundary was completely off. To fix this, we simply add a grid id (automatically populated in python) and propagated to the cpp backend via swig.
Practically speaking, this means that the unit of differentiation is now the design variable, not the grid copy.
Turns out this bug also affects symmetries (see #1984) and the fix also resolves it.