solve_eigfreq tutorial

[stevengj]

Now that #1158 is merged, it might be nice to use solve_eigfreq in a tutorial.

Typically, to find resonances we do multiple simulations: first we use harminv with a broadband pulse, to identify the resonant frequencies. Then we do additional calculations with narrow-band pulses around each resonant frequency in order to see individual modes and also to get more accurate harminv data.

Instead, you can now replace the second step by calling solve_eigfreq with a guessfreq given by the broadband harminv result. This might be faster than a very narrowband calculation, I'm not sure yet (will probably depend on the problem). But it can probably give a more accurate result, especially for the Q of the mode (since high-Q calculations correspond to finding the complex ω to many decimal places).

(For example, suppose you have a resonator with a Q=10⁵, then if you want 3 digits of accuracy in Q you need to compute the complex ω to about 8 digits. This might require a very long simulation for Harminv, if it is possible at all.)

[oskooi]

In the following example of a ring resonator for a mode with Q of ~10⁶ (computed using harminv via a narrowband source with an error of ~10^-10), the mode frequency/Q pair computed by solve_eigfreq (with a guessfreq from harminv) has significantly worse accuracy (factor of 10⁴!) and requires a substantially longer runtime than harminv (1483.8517s vs. 182.8247s using two cores). Increasing the solve_eigfreq convergence parameters maxiters and cwmaxiters improves accuracy at the expense of longer runtime.

From these results, it seems that the eigensolver is impractical for these applications.

harminv result

harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error
harminv0:, 0.25763966078694694, -5.938000849363818e-07, 216941.41456257098, 0.06745820629360579, -0.06655029897341563+0.011030290245486258i, 2.0652221851419098e-10+0.0i

eigensolver result

Eigensolver step 30: 80 CG iters, freq = 0.256138-0.000940994i (change = 0.00169493+0.00010582i).
 -- CONVERGENCE FAILURE (1) in solve_cw!
eigfreq:, 0.25613766961518836, 136.0995780261287

import meep as mp

n = 3.4                 # index of waveguide                                                                                          
w = 1                   # width of waveguide                                                                                          
r = 1                   # inner radius of ring                                                                                        
pad = 4                 # padding between waveguide and edge of PML                                                                   
dpml = 2                # thickness of PML                                                                                            
sxy = 2*(r+w+pad+dpml)  # cell size                                                                                                   

c1 = mp.Cylinder(radius=r+w, material=mp.Medium(index=n))
c2 = mp.Cylinder(radius=r)

fcen = 0.257             # source center frequency                                                                                    

src = mp.Source(mp.ContinuousSource(fcen),
                component=mp.Ez,
                center=mp.Vector3(r+0.1))

sim = mp.Simulation(cell_size=mp.Vector3(sxy, sxy),
                    geometry=[c1, c2],
                    sources=[src],
                    resolution=40,
                    force_complex_fields=True,
                    symmetries=[mp.Mirror(mp.Y)],
                    boundary_layers=[mp.PML(dpml)])

sim.init_sim()
eigfreq = sim.solve_eigfreq(tol=1e-7, maxiters=30, guessfreq=fcen, cwtol=1e-6, cwmaxiters=80, L=10)

if mp.am_master():
    print("eigfreq:, {}, {}".format(eigfreq.real,eigfreq.real/(-2*eigfreq.imag)))

By comparison, for a mode with smaller Q of ~10⁴ the eigensolver produces results with slightly better accuracy but still far from adequate.

harminv result

harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error
harminv0:, 0.20345453628737467, -1.2503433555924088e-05, 8135.946633274128, 0.040612299864003616, -0.02673202846745665-0.030573805034029296i, 6.3668396774220965e-12+0.0i

eigensolver result

Eigensolver step 30: 80 CG iters, freq = 0.178157-8.45776e-05i (change = 0.0159243+2.18228e-05i).
 -- CONVERGENCE FAILURE (1) in solve_cw!
eigfreq:, 0.1781570997298217, 1053.2162229720261