multimode source injection

[oskooi]

Currently, the eigenmode source only permits launching a single mode specified either by the frequency (eig_match_freq=True) or wavevector (eig_kpoint). Since the time profile of the source is often a Gaussian pulse (which is necessary for calculations involving Fourier-transformed fields such as mode coefficients or S-parameters), It would be useful to be able to launch multiple modes simultaneously either at equal or unequally-spaced frequency points (similar to #1070). This would enable computing in a single broadband simulation what now requires multiple single-frequency calculations.

[stevengj]

Note that you can already just add multiple eigenmode sources at the same position, for multiple modes and/or frequencies.

For a single mode, you can just add a bunch of narrow-band sources at different frequencies. The disadvantage of this is that the runtime increases as you add more points (narrowing the bandwidth of each).

[stevengj]

Various strategies and tradeoffs are also discussed in section 4.2.2 of our paper: https://arxiv.org/abs/1301.5366

[smartalecH]

One possible implementation (as described in your paper) is to specify a filter from the dispersive eigenmode coefficients. There are several FIR filter design tools out there (e.g. Parks-Mcclellan) but I think this approach may work better:

Once we have a Chebyshev interpolator in place (see #1070) we can allow the user to specify however many frequencies they care about, along with the start/stop points of the band of interest. The algorithm will add a filter tap for every point specified by the user. Since it's using a Chebyshev basis, it should handle Runge effects along the edges much better than a simple polynomial FIR filter fit. We would have to transform the continuous-time coefficients to discrete-time (i.e. using the bilinear transform in #920). We also need to implement a filter routine for the sources, obviously.

I think this approach would work better than a standard FIR filter tool, would be easier to implement (since it leverages features from other PRs), and would always be stable (in contrast to a Pade approximant or IIR filter).

It would consume more taps and therefore more memory, but the user just has to understand that the memory consumption should increase linearly with the number of frequency points. Ideally, with the Chebyshev fit, the user would need far fewer frequency points than the brute force method also specified in the paper.

[stevengj]

My understanding is that Parks–McClellan is not susceptible to Runge phenomena — in fact, my understanding is that it is a variant of Remez exchange and gives the optimal filter in the L∞ sense (also called "Chebyshev" optimal, confusingly) of a given degree. In contrast, interpolation from the Chebyshev nodes, is only an approximation for the minimax polynomial (but can be computed much more quickly).

[oskooi]

One useful demonstration that may be worth adding to Tutorial/Eigenmode Source: when launching an eigenmode via eig_match_freq=True using a Gaussian pulse for a given single-mode ridge waveguide, compute the error in the mode profiles at the non-center frequencies. Depending on the bandwidth of the pulse and the band structure of the waveguide near the center frequency, perhaps using a single pulsed eigenmode source (rather than multiple individual eigenmode sources with different center frequencies as proposed by @stevengj above) is sufficient for launching eigenmodes over a given bandwidth spectrum?

[stevengj]

The simplest demonstration would be to launch a broad-band eigenmode source into a simple dielectric waveguide.

1. In the single-mode case, also compute the Poynting flux through the sides of the computational cell, and also compute the power going in the backwards direction, as well as the eigenmode coefficient in the forward direction. Normalize everything by the total power. You will see that, away from the center frequency, the power in the radiative and backward directions both increase.

2. In a multi-mode case, you could also compute power into other modes (of the same symmetry) in the transmitted directions. Again, you will see that they increase away from the center frequency.

[oskooi]

The following are results for the demonstration outlined above for the single- and multi-mode waveguide. As expected, the backward/scattered power at the off-center frequencies is significantly larger for the multi-mode waveguide compared to the single-mode waveguide.

import meep as mp
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt

resolution = 50 # pixels/μm                                                                                                                               

sxy = 10
cell_size = mp.Vector3(sxy,sxy)

dpml = 1
pml_layers = [mp.PML(thickness=dpml)]

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

w = 1.0 # width of waveguide                                                                                                                              

geometry = [mp.Block(center=mp.Vector3(),
                     size=mp.Vector3(mp.inf,w,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 source                                                                                                               
kx = 0.4    # initial guess for wavevector in x-direction of eigenmode                                                                                    
bnum = 1    # band number of eigenmode                                                                                                                    

df = 0.1    # frequency width                                                                                                                             
nfreq = 51  # number of frequencies                                                                                                                       

kpoint = mp.Vector3(kx).rotate(mp.Vector3(z=1), rot_angle)

symmetries = [mp.Mirror(mp.Y)]

sources = [mp.EigenModeSource(src=mp.GaussianSource(fsrc,fwidth=df),
                              center=mp.Vector3(0,0),
                              size=mp.Vector3(0,sxy),
                              direction=mp.NO_DIRECTION,
                              eig_parity=mp.EVEN_Y+mp.ODD_Z,
                              eig_kpoint=kpoint,
                              eig_band=bnum,
                              eig_match_freq=True)]

sim = mp.Simulation(cell_size=cell_size,
                    resolution=resolution,
                    boundary_layers=pml_layers,
                    sources=sources,
                    geometry=geometry,
                    symmetries=symmetries)

flux_mon_top = sim.add_flux(fsrc,df,nfreq,mp.FluxRegion(center=mp.Vector3(0,0.5*sxy-dpml),size=mp.Vector3(sxy-2*dpml,0),weight=+1))
flux_mon_bot = sim.add_flux(fsrc,df,nfreq,mp.FluxRegion(center=mp.Vector3(0,-0.5*sxy+dpml),size=mp.Vector3(sxy-2*dpml,0),weight=-1))
flux_mon_right = sim.add_flux(fsrc,df,nfreq,mp.FluxRegion(center=mp.Vector3(0.5*sxy-dpml,0),size=mp.Vector3(0,sxy-2*dpml),weight=+1))
flux_mon_left = sim.add_flux(fsrc,df,nfreq,mp.FluxRegion(center=mp.Vector3(-0.5*sxy+dpml,0),size=mp.Vector3(0,sxy-2*dpml),weight=-1))

sim.run(until_after_sources=50)

freqs = mp.get_flux_freqs(flux_mon_top)
flux_top = np.asarray(mp.get_fluxes(flux_mon_top))
flux_bot = np.asarray(mp.get_fluxes(flux_mon_bot))
flux_right = np.asarray(mp.get_fluxes(flux_mon_right))
flux_left = np.asarray(mp.get_fluxes(flux_mon_left))

res = sim.get_eigenmode_coefficients(flux_mon_right,[bnum],eig_parity=mp.EVEN_Y+mp.ODD_Z)
coeffs = res.alpha[0,:,0]

flux_tot = flux_top + flux_bot + flux_right + flux_left
refl = flux_left / flux_tot
scatt = (flux_tot - np.power(np.abs(coeffs),2)) / flux_tot

if mp.am_master():
    fig = plt.figure()
    plt.subplot(1,2,1)
    plt.plot(freqs,refl,'bo-')
    plt.xlabel('frequency')
    plt.ylabel('backward power (fraction of total power)')
    plt.subplot(1,2,2)
    plt.plot(freqs,scatt,'ro-')
    plt.xlabel('frequency')
    plt.ylabel('scattered power (fraction of total power)')
    fig.subplots_adjust(wspace=0.5, hspace=0)
    fig.suptitle("single mode waveguide with pulsed eigenmode source\n (center freq. = {})".format(fsrc))
    plt.savefig('single_mode_eigsource.png',dpi=150,bbox_inches='tight')

Results for the single-mode waveguide were obtained with a small modification to the script involving the eigenmode properties of the source:

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

[stevengj]

Looks good! As expected, it has a minimum of ≈ 0 at the center frequency.

(But the errors are pretty small even away from the center.)

[stevengj]

If you want an example with a big error away from the center frequency, try a periodic waveguide (as in our holey-waveguide tutorial) and put the center frequency 10% below the band edge of the first band. You should see that the error diverges (because the dispersion diverges) as the band edge is approached.

[smartalecH]

It now makes sense to leverage our fitting routines to enable an eigenmode source for multiple frequencies, similar to what's being done in #1329. This would be especially useful for broadband adjoint simulations that occur near band edges.

The machinery for sources seems slightly different than near2far transforms, however.

What's the recommended way to store all of the eigenmode fields on the yee grid, pass those fields to the fitting in routine in python, and then back to the add_volume_source constructor?

[stevengj]

You could call add_source as usual, but then strip off the sources most recently added to the fields chunk.

Probably the easiest way to do this is:

1. Save savesources[i] = f->chunks[i]->sources to a temporary array savesources[i].
2. Set f->chunks[i]->sources = NULL and call add_sources with the eigenmode source at the current frequency.
3. Grab the amplitudes and indices arrays from f->chunks[i]->sources.
4. Restore f->chunks[i]->sources = savesources[i].