Waveguide crossing inverse design example

[smartalecH]

Here we add a python inverse-design example that maximizes the transmission of a waveguide crossing. In particular, it's meant to showcase the new (first-order accurate) smoothed projection functionality. Using this example, you can:

Perform shape optimization:

Perform topology optimization:

Analyze the norm of the gradient as $\beta \to \infty$ for both the smoothed case and non-smoothed case:

Analyze the convergence properties of a shape optimization problem at a finite beta ($\beta=64$ here) for both the smoothed case and non-smoothed case:

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 14.28571% with 24 lines in your changes are missing coverage. Please review.

Project coverage is 73.75%. Comparing base (295fa63) to head (68e7ae1). Report is 26 commits behind head on master.

:exclamation: Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #2800 +/- ## ========================================== - Coverage 74.06% 73.75% -0.32% ========================================== Files 18 18 Lines 5395 5422 +27 ========================================== + Hits 3996 3999 +3 - Misses 1399 1423 +24 ``` | [Files](https://app.codecov.io/gh/NanoComp/meep/pull/2800?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=None) | Coverage Δ | | |---|---|---| | [python/adjoint/filters.py](https://app.codecov.io/gh/NanoComp/meep/pull/2800?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=None#diff-cHl0aG9uL2Fkam9pbnQvZmlsdGVycy5weQ==) | `69.09% <14.28%> (-8.12%)` | :arrow_down: |

[smartalecH]

Here is an example that runs for 1000 iterations. Interesting how the non-smoothing case eventually breaks through the floor and keeps descending.

[mochen4]

As discussed, if it is a greyscale issue, you can try using the damping term in the MaterialGrid, and it really helps with binarization. From my experience, typically damping=0.05ω is sufficient.

[smartalecH]

If we look at the final designs for both cases, we see that the optimizer is not simply relying on making the non-smoothed case less binary. Rather, it continues to perturb the geometry.

[stevengj]

I wonder if it's getting stuck from the subpixel-smoothed structure being only piecewise differentiable?

In principle, we could fix that by convolving with something like a Gaussian rather than a sphere.

[smartalecH]

So here is a plot of the gradients during the "flat region" of the smoothing optimization:

What's interesting, is that the gradients for the smoothed case are antisymmetric across x and y... this is problematic though because our mapping function forces a symmetry across x and y... so the optimizer can't actually achieve the change it wants to make (i.e. there may be a bug somewhere breaking down the gradient accuracy?)

Edit: Actually the mapping function is forcing a C4 symmetry only, so this is consistent... but the optimizer isn't wanting to deviate from mirror symmetries?

[stevengj]

The alternating signs of the gradient would cause the changes to cancel out to first order, so if that's due to a sign error it would be a problem.

[smartalecH]

@stevengj if I monitor the output of CCSA during the optimization, rho indeed blows up.

During the stalled region, it's not necessarily stuck in the same outer iteration. But rho keeps climbing (eg 1e15) and it does take more inner iterations to converge than the rest of the optimization.

[stevengj]

Some things to try:

1. Check smoothness numerically. Define ρ to be e.g. a Gaussian centered at (x,0), which projects to a cylinder, and plot the smoothed/projected ρ̂ at some pixel as a function of x — it should go smoothly (or at least with a continuous 1st derivative) from 0 to 1 as the cylinder passes through the pixel location.

2. If that looks fine, in principle you could do the same experiment with the Meep solution, e.g. plot DFT electric field at some point as a function of the Gaussian center x and see if it is changing smoothly.

3. Replace the smoothing function

   

   with some function of our choosing, e.g. a Gaussian or a twice-differentiable "bump" function with compact support (on a sphere). Implicitly, this defines a smoothing kernel that you could reverse-engineer if you wanted (since the smoothing kernel integrated over the spherical cap gives F_c, you should be able to take some derivatives of the latter to get the former)