Yee grid gradients

[smartalecH]

Closes #1270

Before:

After:

[review-notebook-app]

Check out this pull request on 

Review Jupyter notebook visual diffs & provide feedback on notebooks.

---

Powered by ReviewNB

[stevengj]

To automatically fix merge conflicts due to clang-format, see #662:

# update repo
git checkout master && git pull origin master # update master

# switch to your branch
git checkout MYBRANCH

# rebase onto PREVIOUS commit
git rebase ad8986b622631232324fc7a100fb726a297f22f3

# reformat and rebase to clang-format commit
find . -name '*.*pp' | xargs clang-format -i
git commit -a -m "reformat"
git rebase -Xtheirs c363d2e609b7a970279a503d5cb85cb6201ef138

# finish rebase
git rebase master

[stevengj]

@zlin-opt has also been playing around with this, and recently set up an adjoint calculation "by hand" where he directly populated the Yee grid for his example problem and was able to get very good accuracy. It would be good to compare your (more general) code to Zin's for his special case to see if you get equivalent accuracy or if one of you is missing a trick.

[smartalecH]

Sounds great. Do you have a link?

[zlin-opt]

Sounds great. Do you have a link?

I will post my example case soon.

[zlin-opt]

@stevengj @smartalecH
Here is my test case for adjoint gradients on a 3d structure. The DOFs here are the variable heights.

my approach is to directly populate the Yee-grid through epsilon_func and amp_func for any given epsilon and current sources and also use yee-grid dft. I have made a repo for my codes.

[smartalecH]

Looks great, Zin. Thanks for taking the time to pull this together.

I carefully went through your repo (linked for posterity). It seems very consistent with how I implemented meep's "general" adjoint solver. In particular, it's important to use the yee grid for the adjoint source and the forward/adjoint fields in the design region.

This PR enables evaluating the forward/adjoint fields directly on the yee grid, a major improvement that's long overdue. A future PR needs to implement the first half of Zin's efforts (also discussed in #1270): using the amp_func for the adjoint source.

Currently, since the objective arguments are eigenmode coefficients, I simply use another eigenmode source pointing in the opposite direction as the adjoint source (with the correct scaling of course). There's some obvious handwaving here, especially when it comes to the yee grid.

It would be better to "generalize" this by storing the fields and directly implementing the adjoint source on the yee grid. (Homer's original code did this to some extent). This is somewhat problematic, however, for eigenmode arguments, Poynting flux, or anything with a cross product because the corresponding adjoint source uses the "opposite" field quantities (i.e. J sources use H fields and M sources use E fields). This may not play nicely with the yee grid. Is it possible to put a point J source on an H yee grid point? Or are sources agnostic to the yee grid?

Edit: After reviewing the Meep paper, I now realize that meep should automatically project the electric current (J) source specified at the magnetic field (H) locations to the corresponding Ex/Ey/Ez locations using a restriction (i.e. just the transpose of the interpolation matrix). While not a perfect solution, it should still be second-order accurate, and probably the best we can hope for in these scenarios.

This would require some extra machinery, however... so it's probably best to submit that patch in a separate PR.

[stevengj]

@smartalecH, I think the situation with eigenmode sources may be easier than you think; the main thing is to make sure that the overlap integral (for the coefficient extraction) occurs on the Yee grid rather than interpolating the DFT fields to the voxel centers first. Then the subsequent current sources should do the right thing because of the way they are implemented, I think.

[smartalecH]

I've tried evaluating the eigenmode coefficient using the yee grid directly, but the results are actually worse:

Note I also use the DFT fields directly for the adjoin source, but got identical results as and eigenmode source pointed in the opposite direction; so it's definitely something going on with the eigenmode coefficient.

[stevengj]

I think we need to add a centered_grid flag (defaults to true) to add_dft_flux so that you can passes centered_grid=false to add_dft here when you are adding DFT regions for mode coefficient calculations.

(The issue is that the add_dft_flux regions are designed for computing ExH of Meep fields, so it defaults to interpolating to the centered grid. Here, we don't want to do that, because the eigenmode sources do no centered-grid interpolation. It looks like the subsequent get_mode_flux_overlap calculations will be correct—they will use whatever grid and weights were employed when add_dft_flux was called.)

[smartalecH]

I think we need to add a centered_grid flag

This seemed to have done the trick (and was much simpler than my original fix -- which I also learned works after a bugfix). As discussed before, the gradients improve as the Courant factor decreases:

Note that I already implemented the iω discrete-time transform (so it's not that). Rather, I think the issue is with the forward source normalization. The forward source is typically a gaussian and normalized such that the center frequency has unit power. However, the normalization constant is performed using a continuous approximation of the gaussian, rather than the DTFT of the truncated gaussian (remember -- all of the Fourier transformed fields in meep are calculated using the DTFT).

This complex amplitude factor is used to normalize out the Gaussian pulse in the adjoint source -- but it's wrong. It makes sense that the adjoints are more accurate with smaller timesteps (but equal resolution) because this factor more closely approximates the continuous transform.

It would be nice to analytically derive the dtft of the truncated gaussian. For now, I'll try a brute force calculation and use that in my normalization factor.

[smartalecH]

Okay I think I got it. It's not the most elegant approach but it works well (Courant=0.5):

Next steps are to clean up the failing tests and resolve a few bugs in 3D I found during an unrelated exercise.

[smartalecH]

Note that a grating coupler example still produces an incorrect gradient by a common scalefactor:

The scale factor decreases with resolution. I'll post more information once I have it.