oblique-waveguide support in mode decomposition

[stevengj]

Once #675 is merged, it would be good to support the same feature in the mode-decomposition code.

Hopefully it will be straightforward, but may need a little thought about how mode-orthogonality works in this kind of oblique setting, since that's not the usual setting where this property is derived.

[stevengj]

In particular, @fesc3555 already worked out that the Poynting flux is invariant under coordinate transformations (https://github.com/fesc3555/meep_variable_resolution/blob/master/notes.pdf). It would be worthwhile going through a similar procedure to see if the mode inner product ∫(E₁×H₂ + E₂×H₁)⋅dA has the same feature.

[stevengj]

Looking through Felix's notes, I think it works out — he shows that ∫E×H⋅dA is invariant under coordinate transformations, but does not require that E and H belong to the same mode (i.e. one doesn't have to be the curl of the other), only that they transform in the same way. This should mean that the modal inner product is invariant as well.

As a result, I think we can use the existing mode-decomposition feature as-is as long as we plug the oblique-waveguide modes into it.

[stevengj]

In particular, it mainly looks like we need a way to pass NO_DIRECTION to get_eigenmode inside the mode-coefficient calculation, in order to tell it to get its lattice vector from kguess as in #675. This may require adding another parameter to the API.

For now, you could just hack this line to pass NO_DIRECTION instead of d for testing purposes.

[oskooi]

With the suggested one-line hack in src/mpb.cpp:767, the results for the eigenmode coefficients of an oblique waveguide seem to be correct. The example is based on Tutorials/Eigenmode Source/Index-Guided Modes in a Ridge Waveguide and involves comparing the Poynting flux values computed using two methods — get_fluxes and get_eigenmode_coefficients (with the addition of a kpoint_func parameter to specify the wavevector) — in an oblique waveguide at four different rotation angles: 10°, 20°, 30°, and 40°.

import meep as mp
import numpy as np

resolution = 50 # pixels/μm                                                                                                                                                                                

cell_size = mp.Vector3(14,14)

pml_layers = [mp.PML(thickness=2)]

# rotation angle (in degrees) of waveguide, counter clockwise (CCW) around z-axis                                                                                                                          
rot_angle = np.radians(40)

geometry = [mp.Block(center=mp.Vector3(),
                     size=mp.Vector3(mp.inf,1,mp.inf),
                     e1=mp.Vector3(1).rotate(mp.Vector3(z=1), rot_angle),
                     e2=mp.Vector3(y=1).rotate(mp.Vector3(z=1), rot_angle),
                     material=mp.Medium(epsilon=12))]

fsrc = 0.15 # frequency of eigenmode or constant-amplitude source                                                                                                                                          
kx = 0.4    # initial guess for wavevector in x-direction of eigenmode                                                                                                                                     
bnum = 1    # band number of eigenmode                                                                                                                                                                     

sources = [mp.EigenModeSource(src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc),
                              center=mp.Vector3(),
                              size=mp.Vector3(y=14),
                              direction=mp.AUTOMATIC if rot_angle == 0 else mp.NO_DIRECTION,
                              eig_kpoint=mp.Vector3(kx).rotate(mp.Vector3(z=1), rot_angle),
                              eig_band=bnum,
                              eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z,
                              eig_match_freq=True)]

sim = mp.Simulation(cell_size=cell_size,
                    resolution=resolution,
                    boundary_layers=pml_layers,
                    sources=sources,
                    geometry=geometry,
                    symmetries=[mp.Mirror(mp.Y)] if rot_angle == 0 else [])

tran = sim.add_flux(fsrc, 0, 1, mp.FluxRegion(center=mp.Vector3(x=5), size=mp.Vector3(y=14)))
sim.run(until_after_sources=50)

res = sim.get_eigenmode_coefficients(tran,
                                     [1],
                                     eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z,
                                     kpoint_func=lambda f,n: mp.Vector3(kx).rotate(mp.Vector3(z=1), rot_angle))

mode_coeffs = res.alpha

print("flux:, {:.6f}, {:.6f}".format(mp.get_fluxes(tran)[0], abs(mode_coeffs[0,0,0])**2))

rotation angle: 10°

flux:, 1109.973264, 1110.955242

rotation angle: 20°

flux:, 1109.565028, 1111.636355

rotation angle: 30°

flux:, 1110.547906, 1111.619587

rotation angle: 40°

flux:, 1108.759159, 1109.236536

These values are all within less than 2% of one another. For comparison, the following are the results from master.
rotation angle: 10°

flux:, 1109.973264, 1093.363832

rotation angle: 20°

flux:, 1109.565028, 1042.771962

rotation angle: 30°

flux:, 1110.547906, 958.539646

rotation angle: 40°

flux:, 1108.759159, 825.569254

The difference in the pairs of values increases with the rotation angle (it is ~26% for a rotation angle of 40°).

Since EigenmodeSource contains a parameter direction which can be set to mp.AUTOMATIC or mp.NO_DIRECTION, it seems that this should also be added get_eigenmode_coefficients for consistency.