Adjoint optimization through Poynting Flux

[MalharJ]

Is it possible to perform adjoint inverse design through the Poynting Flux over a certain range of frequencies? The meep.adjoint.EigenmodeCoefficient object has a register_monitors function to track the power across specified frequencies, but I'm unsure as to how one could use this to inverse design a device.

I tried to use it in this tutorial, where I added the following lines in the adjoint optimization tutorial:

frequencies = np.linspace(1/1.55 - 0.1, 1/1.55 + 0.1, 100)
TE_top.register_monitors(frequencies)

https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/adjoint_optimization/01-Introduction.ipynb

I replaced the TE0 variable with TE0.monitor, but this ended up giving me the following error:

AttributeError: 'DftFlux' object has no attribute 'register_monitors'.

[smartalecH]

As of right now, you can't use anything but a ModeMonitor in your objective function.

Using the Poynting flux (or any other quantity) is a bit tricky since the adjoint source varies in both time and space. We can get around this with the mode monitor by just using a mode source as the adjoint source.

We do have plans to add support for Poynting Flux, though. We have the framework in place to fit the adjoint source to the correct data. It's just not super high on our priorities -- especially since you can still use mode monitors in place of Poynting Flux monitors in several instances.

[stevengj]

(Note that the time/space-varying source support will have to be added anyway by @mochen4 at some point in support of the near2far adjoint feature (#1300), and it's conceptually straightforward as we discussed — we fit the frequency spectrum at each point to a sum of gaussians as we do now, but we treat each gaussian as a spatially distributed source with a position-dependent amplitude given by the fit at each point.)

[MalharJ]

2 follow up questions:

• By ModeMonitor, do you mean the meep.adjoint.EigenmodeCoefficients class?
• Is there any way to specify the wavelength of the EigenmodeCoefficients in the objective function? Like if I wanted to design a 1x2 WDM?

[smartalecH]

Take a look at the splitter example.

You can modify it to do WDM filtering.

[MalharJ]

I'm aware of the splitter example, but it doesn't mention how we can specify the wavelengths in the objective function, which was my original question.

[stevengj]

(@smartalecH, strictly speaking, the adjoint of the eigenmode coefficients should also be varying in both space and frequency. That's because the eigenmode coefficients recompute the mode profile at each frequency, whereas the eigenmode source computes just the mode profile at the center frequency and neglects modal dispersion.)

[MalharJ]

@smartalecH @stevengj could you please clarify on how we can modify the example to get a WDM? I thought about it, and there seem to be 2 approaches:

• We can run two separate simulations for each frequency band with the same objective function defined in terms of the frequencies
• We can use the frequencies argument of register_monitors.

[smartalecH]

There are lots of ways to do this. Your first suggestion would work. Your second suggestion won't give you what you want -- you don't need to modify the internals of the adjoint solver to get the gradient.

[smartalecH]

Keep in mind that the adjoint solver is under constant development right now and is rather "bleeding edge". Consequently, there's not much support to offer right now...especially when every major PR dramatically changes the adjoint codebase.

Hopefully, we can get some docs together soon.

[MalharJ]

@smartalecH I don't want to modify the internals of the adjoint solver to get the gradient though. I just want to select the eigenmodes based on their wavelength. Can I currently do this with the adjoint solver without 2 separate simulations?

[yunxiangw]

I am also interested in doing inverse design through the Poynting flux.
From my understanding, the poynting flux is defined in the frequency domain:

So the adjoint source should be:

where n is the normal vector

For example, if we consider a 2D simulation. The objective function is the poynting flux through a line that is normal to the y axis.
The adjoint source should be:

If we consider the flux over a range of frequencies, we have to fit the adjoint source profile in the frequency domain.
But if we only consider single frequency, the adjoint source should be a lot of gaussian sources with complex amplitude H*. And we can use mp.Source to define the adjoint source.
Am I correct?

Since the Poynting flux need both E and H, we have to interpolate the field to a common point (center of each yee cell). The adjoint source is also defined at the center of each yee cell. The linear interpretation will introduce some error. Is there any way to reduce the error introduced by the interpolation?

[stevengj]

Since the Poynting flux need both E and H, we have to interpolate the field to a common point (center of each yee cell). The adjoint source is also defined at the center of each yee cell. The linear interpretation will introduce some error. Is there any way to reduce the error introduced by the interpolation?

There doesn't have to be any error — you just need to multiply by the transpose of this interpolation matrix when you construct the adjoint source...

[smartalecH]

If we consider the flux over a range of frequencies, we have to fit the adjoint source profile in the frequency domain.
But if we only consider single frequency, the adjoint source should be a lot of gaussian sources with complex amplitude H*. And we can use mp.Source to define the adjoint source.
Am I correct?

The key here is if we only consider a single frequency use case. The adjoint solver is set up to work over a broad range of frequencies. Sure, you could pull together what you described (and modify quite a bit of the existing framework along the way), but I feel it's better to do it right from the beginning.

As @stevengj mentioned earlier, take a look at #1300. Rather than place Gaussians at each point, we are using the FilteredSource time profile at each point. Unfortunately, this doesn't scale well for large problems. We need to port a lot of this to C.

It's not a ton of work -- I just don't have the bandwidth to do it right now. Plus, as I said, eigenmode coefficients work in place of flux in several use cases anyway (at least all the designs I'm working on).

That being said, feel free to contribute a PR. Luckily the groundwork for everything is already in place -- there's no new physics or math... just some software dev. And it seems like you understand what needs to happen rather well.

[stevengj]

Unfortunately, this doesn't scale well for large problems. We need to port a lot of this to C.

Not necessarily a lot. I'm working with @mochen4 on a new API that exposes functionality from the internal srcvol structure, which will allow us to specify sources with an array of amplitudes directly onto the Yee grid. Given this, you will still be able to do the filtering/fitting process in Numpy efficiently (it's a least-square fit with a matrix right-hand-side, which vectorizes well), and use that to create a small number of Gaussian sources with spatially varying amplitudes. (I think it's desirable to maintain the flexibility of designing the filter/time-dependence in Python for now.)

(This is necessary to efficiently support parallel near2far adjoint sources, so that each process only needs to compute the fits for the process-local adjoint sources. The near2far code in C++ has easy access to the exact indices of the Yee-grid locations where the sources are required, so it's just a matter of passing these back to Python and exposing a srcvol-like structure so that Python can add them back as currents. The same thing should make life simpler for optimizing the DFT flux code and other DFT-field functions — it's easy to access Yee-grid indices in C++, so we just need to pass this back to Python.)

[smartalecH]

@stevengj you're right, the fit itself isn't computationally expensive.

The expensive part is having a python callback function for each point that has its own sum of Nuttall windows (we opted to use the Nuttall window over the continuous time Gaussians).

[stevengj]

The expensive part is having a python callback function for each point that has its own sum of Nuttall windows (we opted to use the Nuttall window over the continuous time Gaussians).

With the new architecture, the plan is to have have one Python callback per temporal basis function (e.g. Nutall window), applied to many (all) spatial points at once (possibly per component and per chunk), so Python performance shouldn't be a big issue even if you have a sum of 100s of basis functions.

[yunxiangw]

I realized the Poynting flux objective function 'by hand'. The result seems reasonable.


Right now I want to use the new feature #1285 to make the code more clear. But I found that this merged PR has not been included in the current nightly-build version. When will this new feature be included?

[smartalecH]

Do you want to submit a PR using the in-place adjoint framework?

[yunxiangw]

Ok, I can do it using the current avaiable version.