bad interpolation from Meep grid to MPB during eigenmode calculation

[oskooi]

The following example tries to launch an eigenmode source in a holey waveguide in a 2d simulation but aborts with an error due to MPB. The same error occurs when trying to use mode decomposition for a 2d eigenmode. Neither EigenModeSource or get_eigenmode_coefficients has previously been tested for this type of use case. All the examples in Tutorial/Mode Decomposition and Tutorial/Eigenmode Source are in 1d.

The schematic below (obtained using plot2D) shows the 2d source (red hatches) and mode monitor (blue hatches).

additional info: removing the symmetries from the Simulation constructor does not solve the problem. However, changing the size of the source/monitor from a 2d rectangle to a line in y (i.e. size=mp.Vector3(0,sy)) works fine.

import meep as mp
import numpy as np

resolution = 50 # pixels/μm                                                                                                                                          

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

eps = 13          # dielectric constant of waveguide                                                                                                                 
w = 1.2           # width of waveguide                                                                                                                               
r = 0.36          # radius of holes                                                                                                                                  
N = 11            # number of holes in waveguide                                                                                                                     

sx = N
sy = 10
cell_size = mp.Vector3(sx,sy)

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

for n in range(N):
    geometry.append(mp.Cylinder(radius=r,
                                center=mp.Vector3(-0.5*sx+0.5+n),
                                height=mp.inf,
                                material=mp.air))

f_bandedge = 0.20025  # from MPB                                                                                                                                     
fsrc = 0.9*f_bandedge # 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                                                                                                                              

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

sources = [mp.EigenModeSource(src=mp.GaussianSource(fsrc,fwidth=df),
                              center=mp.Vector3(0,0),
                              size=mp.Vector3(1,sy),
                              direction=mp.NO_DIRECTION,
                              eig_parity=mp.ODD_Y+mp.EVEN_Z,
                              eig_kpoint=mp.Vector3(kx),
                              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)

sim.run(until_after_sources=100)

output

Using MPI version 3.1, 1 processes
-----------
Initializing structure...
Halving computational cell along direction y
time for choose_chunkdivision = 0.000737682 s
Working in 2D dimensions.
Computational cell is 11 x 10 x 0 with resolution 50
     block, center = (0,0,0)
          size (1e+20,1.2,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (13,13,13)
     cylinder, center = (-5,0,0)
          radius 0.36, height 1e+20, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (1,1,1)
     cylinder, center = (-4,0,0)
          radius 0.36, height 1e+20, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (1,1,1)
     cylinder, center = (-3,0,0)
          radius 0.36, height 1e+20, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (1,1,1)
     cylinder, center = (-2,0,0)
          radius 0.36, height 1e+20, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (1,1,1)
     cylinder, center = (-1,0,0)
          radius 0.36, height 1e+20, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (1,1,1)
     cylinder, center = (0,0,0)
          radius 0.36, height 1e+20, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (1,1,1)
     cylinder, center = (1,0,0)
          radius 0.36, height 1e+20, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (1,1,1)
     cylinder, center = (2,0,0)
          radius 0.36, height 1e+20, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (1,1,1)
     cylinder, center = (3,0,0)
          radius 0.36, height 1e+20, axis (0, 0, 1)
          dielectric constant epsilon diagonal = (1,1,1)
     ...(+ 2 objects not shown)...
time for set_epsilon = 0.499155 s
-----------
Traceback (most recent call last):
  File "holey-wvg-eigsource-pulse-bug.py", line 54, in <module>
    sim.run(until_after_sources=100)
  File "/usr/local/lib/python3.5/site-packages/meep/simulation.py", line 2286, in run
    self.init_sim()
  File "/usr/local/lib/python3.5/site-packages/meep/simulation.py", line 1338, in init_sim
    self.add_source(s)
  File "/usr/local/lib/python3.5/site-packages/meep/simulation.py", line 1624, in add_source
    add_eig_src()
  File "/usr/local/lib/python3.5/site-packages/meep/__init__.py", line 3868, in add_eigenmode_source
    return _meep.fields_add_eigenmode_source(self, *args)
RuntimeError: meep: invalid dielectric function for MPB

[stevengj]

You need direction=mp.X, but I'm not sure why it's giving the "invalid dielectric function" error.

[oskooi]

Even with direction=mp.X in the EigenModeSource definition, the same error occurs:

  File "/usr/local/lib/python3.5/site-packages/meep/__init__.py", line 3868, in add_eigenmode_source
    return _meep.fields_add_eigenmode_source(self, *args)
RuntimeError: meep: invalid dielectric function for MPB

[stevengj]

I don't think it likes the source extending all the way to the edge. Try size=mp.Vector3(1,sy-2*dpml)

[oskooi]

No, same error again unfortunately.

[oskooi]

If I remove the Cylinder objects from the geometry then it seems to work.

[stevengj]

It seems to be a bad interaction with subpixel averaging. Subpixel averaging gives an anisotropic ε on the Yee grid, but we have to interpolate these ε tensor components to get the corresponding MPB values, and that interpolation seems to be destroying positive-definiteness.

[stevengj]

Note that in the case where this does not abort, then there is no need to worry — at worst, some isolated pixels will be messed up, but it will still converge with resolution.

[stevengj]

To quantify the mode mismatch between MPB and Meep, you could compute the backward-propagating flux from an eigenmode source (at the center frequency), which should go to zero with resolution.

[oskooi]

As shown in the results below for the mode at the pulse center frequency of 0.15 in a single mode waveguide, the backward-propagating flux does not seem to be converging to zero with increasing resolution. Rather, it seems to be converging to a constant value. (The results are converged with PML thickness and runtime as well as eig_resolution and eig_tolerance for the EigenModeSource.) Turning off subpixel smoothing reduces the oscillations in the data.

Is there some other kind of discretization effect involved when calling MPB from Meep that we are overlooking?

[stevengj]

You might need to increase the runtime and the cutoffs for the gaussians to avoid spectral leakage, since you are trying to resolve pretty small backward fluxes here.

The difference between Meep and MPB should converge with resolution here, although you'll eventually hit a noise floor due to the eigensolver tolerance, the eigensolver source size, spectral leakage, PML reflections, floating-point roundoff, and other errors.

(It's not clear that you've hit a noise floor here since it's still going down, albeit slowly.)

[stevengj]

This code from #919 looks suspicious to me (cc @smartalecH). It takes all of the tensor components from the same idx despite the fact that they are different points in the Yee grid.

Try changing this line to if (1) in order to disable the dispersive code and see if that improves the behavior above.

[oskooi]

There does not seem to be any change in the results to either the RuntimeError: meep: invalid dielectric function for MPB reported for the holey waveguide (with subpixel smoothing) or the backward-propagating flux for the single-mode ridge waveguide when changing src/monitor.cpp:261 with if (1).

[smartalecH]

It takes all of the tensor components from the same idx despite the fact that they are different points in the Yee grid

I agree this is a problem... I totally forgot about this issue at the time.

It seems odd to me that we are sampling an already sampled yee grid to generate the MPB grid. This not only makes it incredibly difficult to handle the staggering, but it means that "upsampling" the MPB resolution (e.g. in get_eigenmode_coeffcients()) requires extensive linear interpolation on a lower resolution grid.

Why not sample directly from the geometry tree like we do in meep-geom? This would dramatically simplify the code (no more ugly interpolation across grid points by looping through chunks) and it lets us more accurately upsample the mpb grid if desired.

It would be easy to parallelize the grid initialization, as we discussed before. Doing things this way may even solve the weird issue with sub-pixel averaging.

[stevengj]

Why not sample directly from the geometry tree like we do in meep-geom?

mpb.cpp can't assume that the Meep geometry came from a geometry tree; it might have been set directly with the C++ interface, for example.

I suppose we could have optional support for meepgeom in mpb.cpp in cases where this could work.