shift-and-invert eigensolver based on solve_cw

[stevengj]

Closes #1153. Adds a sim.solve_eigfreq method similar to solve_cw except that it returns an estimated eigenfrequency (and sets the fields to the corresponding eigenmode) based on shift-and-invert iteration around the source frequency.

Currently untested, just an initial stab at this. Needs:

• Tests
• Documentation

[stevengj]

Update: it seems to work!

I tried it on the holey waveguide tutorial to solve for the resonance mode, and I got eigfreq = 0.23454890203771478+0.0005581153769037464j (Q ≈ 210) to about 6–7 digits after 7 iterations of shift-and-invert (≈ 50000 timesteps), using a very crude shift (off by 6% from eigenfrequency).

(Caveat: I'm not sure why it's giving the eigenfrequency as having a positive imaginary part, when it should be negative, so I'll need to look into this. But at least it is converging, and plotting the field pattern confirms that it seems to be finding the resonant mode.)

import meep as mp

resolution = 20   # pixels/um

eps = 13          # dielectric constant of waveguide
w = 1.2           # width of waveguide
r = 0.36          # radius of holes
d = 1.4           # defect spacing (ordinary spacing = 1)
N = 5        # number of holes on either side of defect

sy = 6            # size of cell in y direction (perpendicular to wvg.)
pad = 2           # padding between last hole and PML edge
dpml = 1          # PML thickness

sx = 2*(pad+dpml+N)+d-1  # size of cell in x direction

cell = mp.Vector3(sx,sy,0)

blk = mp.Block(size=mp.Vector3(mp.inf,w,mp.inf), material=mp.Medium(epsilon=eps))
geometry = [blk]

for i in range(N):
        geometry.append(mp.Cylinder(r, center=mp.Vector3(d/2+i)))
        geometry.append(mp.Cylinder(r, center=mp.Vector3(-(d/2+i))))
        
pml_layers = [mp.PML(1.0)]

fcen = 0.25
df = 0.2

src = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                 component=mp.Hz,
                 center=mp.Vector3(0),
                 size=mp.Vector3(0,0))]

sym = [mp.Mirror(mp.X, phase=-1), mp.Mirror(mp.Y, phase=-1)]

sim = mp.Simulation(cell_size=cell, force_complex_fields=True,
                    geometry=geometry,
                    boundary_layers=pml_layers,
                    sources=src,
                    symmetries=sym,
                    resolution=resolution)

eigfreq = sim.solve_eigfreq(L=20, cwtol=1e-11, tol=1e-6)

[stevengj]

Hmm, now that I think about it, I'm not estimating the eigenvalue correctly, because the operator here is the discrete time-difference operator and not the ∂/∂t operator.

[stevengj]

Okay, I've pushed a patch that corrects the eigenvalue estimate for the time discretization. It still converges in 7 shift-and-invert iterations, but now the estimated frequency is 0.23453881435117674-3.0147532052022876e-05j (Q ≈ 3889.85), which has the correct sign of the imaginary part.

If I give it the initial guess frequency 0.2351 computed by Harminv in the tutorial, then it converges to nearly the same eigenvalue estimate 0.23453876201337692-3.014712616184812e-05j, but in only 3 shift-and-invert iterations.

(Note that the Q is higher than in the tutorial because I have more holes surrounding the cavity.)

[stevengj]

If I change the number of holes to N = 3 to match the tutorial, then I get a frequency of 0.23445413145107227-0.0003147776872265489j (Q = 372.4), versus 0.235109393214226--3.14979827803982e-4j (Q ≈ 373.2) from Harminv. (Lowering eigtol to 1e-8, I get the same eigenvalue to 9 digits.)

So, this is consistent with the Harminv result but is presumably more accurate.

[oskooi]

I have noticed that solve_eigfreq takes considerably longer for a parallel simulation than a serial run. The following holey-waveguide cavity example takes 41.9494 s using a single-core Kaby Lake 4.2 GHz but 320.6739 s using two cores. More significantly, the eigenfrequency results are different: 0.2344541320888902-0.0003147773149478817j (single core) vs. 0.22236758070320453-0.018193974796342j (two cores).

import meep as mp

w = 1.2           # width of waveguide                                                                                                             
r = 0.36          # radius of holes                                                                                                                
d = 1.4           # defect spacing (ordinary spacing = 1)                                                                                          
N = 3        # number of holes on either side of defect                                                                                            
sy = 6            # size of cell in y direction (perpendicular to wvg.)                                                                            
pad = 2           # padding between last hole and PML edge                                                                                         
dpml = 1          # PML thickness                                                                                                                  
sx = 2*(pad+dpml+N)+d-1  # size of cell in x direction                                                                                             

geometry = [mp.Block(size=mp.Vector3(mp.inf,w,mp.inf), material=mp.Medium(epsilon=13))]
for i in range(N):
    geometry.append(mp.Cylinder(r, center=mp.Vector3(d/2+i)))
    geometry.append(mp.Cylinder(r, center=mp.Vector3(-(d/2+i))))

fcen = 0.25
src = [mp.Source(mp.ContinuousSource(fcen),
                 component=mp.Hz,
                 center=mp.Vector3(0),
                 size=mp.Vector3(0,0))]

sim = mp.Simulation(cell_size=mp.Vector3(sx,sy),
                    force_complex_fields=True,
                    geometry=geometry,
                    boundary_layers=[mp.PML(1.0)],
                    sources=src,
                    symmetries=[mp.Mirror(mp.X, phase=-1), mp.Mirror(mp.Y, phase=-1)],
                    resolution=20)
sim.init_sim()
eigfreq = sim.solve_eigfreq(tol=1e-7)

if mp.am_master():
    print("eigfreq: {}".format(eigfreq))

[stevengj]

Regarding the results being different — this is a bug, I forgot to call sum_to_all in a couple of places. It should be fixed now.

(Hopefully that will fix the speed, too, since the bug probably messed up the convergence rate.)

[stevengj]

One of the tests is timing out, but that is because it was taking 49 minutes and 51 seconds before and so even a few seconds more exceeds the 50-minute time limit. We'll need to split up that test.

[stevengj]

Rebased; the tests should run faster now thanks to #1161.

[smartalecH]

Does this work in cylindrical coordinates too (e.g. for calculating bent waveguide modes?)

What's the advantage of this approach vs the typical harminv approach? If I set up a broadband simulation, I could even "capture" the mode profiles using a DFT monitor at the resonance frequency and normalize out the source (as long as it runs long enough to converge...), right?

[stevengj]

I think cylindrical coordinates should be fine, but I don't recall trying it.

The solve_eigfreq was kind of an experiment, and currently we haven't found any case where it was a clear advantage. Maybe it would be better if we implemented an Arnoldi method on top of it, but the basic difficulty is that the CW solver is slow to converge (though it could potentially be sped up too: #548).