add feature to orthogonalize degenerate modes from mode solver

[dmarek-flex]

draft at the moment for discussion. I added a lot of detail to this discussion, since I am still trying to figure out some of the details for this procedure when applied to mode in transmission lines.

Objective

The purpose of this PR is to improve our handling of degenerate modes. In particular, ensuring that the dot method gives the expected result of close to 0 when the mode indices are different. Before this PR, modes that are degenerate might not be orthogonal wrt the overlap integral calculation in dot. It seems to be worse when using scipy, where in one test calculating the dot of a pair of polarization degenerate modes at 3 different frequencies gave:

Without degenerate mode processing (absolute value of dot)

array([[0.47368421],
       [0.13973138],
       [0.99919483]])

With the degenerate mode processing (absolute value of dot), we get the expected result for the dot of different modes

array([[4.60334385e-16],
       [3.60442433e-17],
       [1.52273941e-14]])

Issue

This all works nicely for optical-type waveguides, at least the square core waveguide I have been using as a test case. But, the results are not as nice when attempting the same operations for a stripline transmission line, where the degeneracy now comes from the fact that there are 2 TEM modes.

Without degenerate mode processing (absolute value of dot)

array([[0.18542736],
       [0.07074218],
       [0.34136264]])

With degenerate mode processing (absolute value of dot)

array([[0.01743628],
       [0.01219984],
       [0.0188109 ]])

It is still unclear to me why these values are so much larger for the stripline. I was able to get a much lower value for dot by introducing the wip dot_yee method, which avoids interpolation and truncation of the grid.

With degenerate mode processing (absolute value of dot_yee)

array([[7.64227775e-17],
       [9.35805564e-16],
       [5.62172366e-16]])

This works because I basically forced the dot_yee method to return the same values I used to form the orthogonal basis in the first place. The problem with this solution is that the flux calculated for each mode will not be the same as the self-terms of dot_yee method, so I also introduced the flux_yee method temporarily so that they match.

Solutions

I see two main ways forward:

1. Just accept that the orthogonality of the degenerate stripline modes will not be as ideal as hoped for
2. Basically modify our calculation for flux everywhere in the code, so that it matches the result of the dot_yee calculation. This modification would essentially avoid interpolation when possible as I did in flux_yee.

Option 2 seems to me more numerically sound, since interpolation of the E field around metallic objects is prone to error, but potentially a lot of work and I foresee a lot of accuracy issue kind of bugs in the near future.

Wondering what everyone thinks @momchil-flex, @weiliangjin2021, @dbochkov-flexcompute, or if maybe I am missing another potential solution? My feeling is I should just go ahead with option 1, since I don't have time at the moment to make these sweeping changes to the flux/dot calculations.

[momchil-flex]

So, this is pretty interesting. My initial reaction to hearing about dot_yee and flux_yee is that these couldn't be correct, because you need to get the E and H fields to be at the same location before computing them. But then I realized that indeed, in-plane, all the relevant terms are already colocated (e.g. Ex and Hy, and Ey and Hx). Your methods do make sense.

Unfortunately, you might be right that switching to this will be annoying. I spent quite a bit of effort making sure that things are consistent for photonics, like FieldData.flux returns the same as FluxMonitor.flux when the two monitors have the same geometry, and same for FieldData.dot(ModeSolverData) vs. ModeData.amps. But I am not sure, maybe it will just work out? It seems like once you switch the .dot definition, the only other thing that needs to change is the on-the-fly flux computation? It might be worth trying that and seeing if all tests pass, and we could take it from there.

Thinking about it a little it could go like:

• Set the colocate argument for flux and field projection monitors to False. Note that these use similar handling also so that we can ensure consistency e.g. between near-field and far-field flux.
• Then, the most nontrivial bit is modifying this function. For example, there will now be different area prefactors for the two different tangential terms.
• Some related modifications might be needed in the solver code, but shouldn't be much there, just incorporating whatever is needed based on the python function modification.
• I think no modification to this function will be needed.

[dmarek-flex]

But I am not sure, maybe it will just work out? It seems like once you switch the .dot definition, the only other thing that needs to change is the on-the-fly flux computation? It might be worth trying that and seeing if all tests pass, and we could take it from there.

Yea I might let this draft PR float here for a bit, but I am interested in testing out the non-interpolation approach when I have a bit more time. Seems like it will be more accurate

[momchil-flex]

But I am not sure, maybe it will just work out? It seems like once you switch the .dot definition, the only other thing that needs to change is the on-the-fly flux computation? It might be worth trying that and seeing if all tests pass, and we could take it from there.

Yea I might let this draft PR float here for a bit, but I am interested in testing out the non-interpolation approach when I have a bit more time. Seems like it will be more accurate

Indeed, I am now curious too, I wonder if it might increase the accuracy in some photonic cases too, e.g. in mode injection/decomposition with angle and/or bend.