symmetry support for mode expansion

[stevengj]

Currently, the mode-expansion calculation doesn't work if there are symmetries bisecting the flux plane. (#293 fixes one problem with symmetry, but there seem to be others.)

• Until it is working, as a short-term fix we could exit with an error if you attempt mode expansion in a simulation with symmetries.

• Once mode-expansion supports parities (add parity argument to eigenmode coefficients #292), we might want to require (or at least document) that the parity should be consistent with any symmetry bisecting the flux plane.

[stevengj]

Probably the call to eigenmode_amplitude should pass S.transform(c, sn) (the logical component) instead of c (the physically stored component in this chunk).

[stevengj]

First things to try (after above fix).

1. Set up a test case, e.g. a 2d waveguide along the x axis from y = -1…1, and with an eigenmode source and flux plane from y=-5…5 at two x values, and an optional y=0 mirror symmetry plane. Goal: get the same answer with and without the mirror plane.

2. Remove the S.reduce call when adding a dft flux object. Hopefully with this change it will just work (at the cost of computing twice as many Fourier transforms).

3. Keep the S.reduce call but pass a parity to MPB in the mode decomposition, telling it to compute just EVEN_Y (or ODD_Y, depending on whether your mirror plane is even or odd) modes. Hopefully with this change it will just work (possibly modulo a factor of 2 in the normalization).

Once we figure out what is working, we can figure out how to expose this to the user. We want to both have the option to handle arbitrary symmetries (i.e. to drop the S.reduce call if needed) and also still retain the factor-of-2 speedup for flux computation (and for mode computation if mirror symmetries are passed to MPB).

[stevengj]

Partially solved by #417 (still lacking a Python API). Would be good to address point 3, above — verify that we can use the symmetry-reduced Fourier transforms as long as the MPB modes of only the corresponding symmetry are computed, at least for mirror planes. It would be good to have this option in performance-critical situations.

[oskooi]

Symmetry objects seem to be working now. The following test, based on the 2d waveguide example described above, produces the same output with and without the mirror symmetry plane.

Script

import meep as mp

sx = 10
sy = 5
cell_size = mp.Vector3(sx,sy,0)

geometry = [mp.Block(size=mp.Vector3(mp.inf,1,mp.inf), center=mp.Vector3(), material=mp.Medium(index=3.5))]

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

fcen = 1
sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen, fwidth=0.2*fcen),
                              component=mp.Ez,
                              size=mp.Vector3(0,sy-2*dpml,0),
                              center=mp.Vector3(-0.5*sx+dpml+0.2),
                              eig_match_freq=True,
                              eig_parity=mp.ODD_Z+mp.EVEN_Y)]

symmetries = [mp.Mirror(mp.Y,phase=1)]

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

xm = 0.5*sx-dpml-0.2
eig_mon = sim.add_mode_monitor(fcen, 0, 1, mp.FluxRegion(center=mp.Vector3(xm,0,0), size=mp.Vector3(0,sy-2*dpml,0)))

sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(xm), 1e-9))

coeffs, vgrps, kpoints = sim.get_eigenmode_coefficients(eig_mon, [1], eig_parity=mp.ODD_Z+mp.EVEN_Y)

print("tran:, {:.8f}".format(abs(coeffs[0,0,0])**2))

Output (with symmetry)

Using MPI version 3.0, 1 processes
-----------
Initializing structure...
Halving computational cell along direction y
Working in 2D dimensions.
Computational cell is 10 x 5 x 0 with resolution 30
     block, center = (0,0,0)
          size (1e+20,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (12.25,12.25,12.25)
time for set_epsilon = 0.0725322 s
-----------
MPB solved for omega_1(3.5,0,0) = 1.00876 after 13 iters
MPB solved for omega_1(3.46911,0,0) = 1 after 2 iters
field decay(t = 50.016666666666666): 0.0750753784221729 / 0.0750753784221729 = 1.0
field decay(t = 100.01666666666667): 0.13204088814763693 / 0.13204088814763693 = 1.0
field decay(t = 150.03333333333333): 1.5119788686314655e-13 / 0.13204088814763693 = 1.1450838371678468e-12
run 0 finished at t = 150.03333333333333 (9002 timesteps)
MPB solved for omega_1(3.5,0,0) = 1.00876 after 16 iters
MPB solved for omega_1(3.46911,0,0) = 1 after 2 iters
tran:, 1.23643620

Elapsed run time = 3.5922 s

Output (without symmetry)

Using MPI version 3.0, 1 processes
-----------
Initializing structure...
Working in 2D dimensions.
Computational cell is 10 x 5 x 0 with resolution 30
     block, center = (0,0,0)
          size (1e+20,1,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (12.25,12.25,12.25)
time for set_epsilon = 0.140093 s
-----------
MPB solved for omega_1(3.5,0,0) = 1.00876 after 13 iters
MPB solved for omega_1(3.46911,0,0) = 1 after 2 iters
on time step 2184 (time=36.4), 0.00183189 s/step
field decay(t = 50.016666666666666): 0.07507537841764578 / 0.07507537841764578 = 1.0
on time step 4176 (time=69.6), 0.00200858 s/step
field decay(t = 100.01666666666667): 0.13204088815118376 / 0.13204088815118376 = 1.0
on time step 6372 (time=106.2), 0.00182162 s/step
on time step 8575 (time=142.917), 0.0018161 s/step
field decay(t = 150.03333333333333): 1.5119856067368922e-13 / 0.13204088815118376 = 1.1450889401816987e-12
run 0 finished at t = 150.03333333333333 (9002 timesteps)
MPB solved for omega_1(3.5,0,0) = 1.00876 after 16 iters
MPB solved for omega_1(3.46911,0,0) = 1 after 2 iters
tran:, 1.23643620

Elapsed run time = 18.5075 s

[stevengj]

Great, would be good to at least document this as an optimization that people can do when it matters.

[stevengj]

Should work with 5666d2f