Open smartalecH opened 3 years ago
Mandelshtam (the author of the harminv) actually published a related algorithm specifically for spectra extrapolation. But the theory of robust pade extrapolation seems better developed.
In general, there are a number of series-extrapolation techniques that we could explore.
Two starting points:
While our FDTD adjoint handles large problems rather well, there are certain cases where the forward and adjoint solves take too long to converge in the frequency domain. For example, quasi-resonant devices or devices that are really large, often require several thousand (or more) time units to converge. The problem obviously compounds with size.
The source of the issue (at least with our implementation) stems from the Fourier uncertainty principle. Slowly decaying signals (tighter bandwidths) must propagate longer to ensure we compute everything correctly.
As we've discussed before (#375), we could get around this by using Pade approximates. The basic idea is that we store the time fields at every point of interest, fit the profiles, "extrapolate", and then compute the DFTs (or some variant of those steps). Fitting to a Pade approximant is rather tricky since it's such an ill-conditioned problem, so there are lots of algorithms to do this (e.g. using the SVD). The advantage here is that we don't need to propagate the fields for nearly as long.
This approach expresses a natural tradeoff between time-stepping and memory requirements. Since every design region point needs to store the fields each timestep, memory usage can become unwieldy (similar to a true-time domain adjoint calculation). But we can cheat here, decimating our stored profiles (or only choosing to store at certain timesteps) such that the sample rate corresponds to frequencies just outside the bandwidth of our objective function (just so we don't alias). The simulation itself is, of course, stable, but we can cut down our memory requirements significantly.
Furthermore, rather than implementing a new fitting algorithm that is "unconditionally" stable, we could instead fit to a series of decaying sinusoids using Harminv. Since the interface is already in place, this should be much easier to debug. The caveat here is that we are assuming the fields everywhere can be expressed in this form.