non-zero kz not passed to get_eigenmode_coefficients for 2d cell

[oskooi]

The k_z component of the k_point is ignored by get_eigenmode_coefficients (taken to be as just 0) for a 2d cell (in x and y) even though k_z is non-zero and the simulation itself is 3d.

As a demonstration of this bug, the binary grating example is modified in the script below to include a non-zero k_z. The k_point had originally been "in the plane" with k_z=0 resulting in a 2d simulation.

import meep as mp
import math
import cmath
import numpy as np

resolution = 40        # pixels/μm                                                                                                            

dpml = 1.0             # PML thickness                                                                                                        
dsub = 1.0             # substrate thickness                                                                                                  
dpad = 1.0             # length of padding between grating and pml                                                                            
gp = 10.0              # grating period                                                                                                       
gh = 0.5               # grating height                                                                                                       
gdc = 0.5              # grating duty cycle                                                                                                   

sx = dpml+dsub+gh+dpad+dpml
sy = gp

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

# replace anisotropic PML with isotropic Absorber to attenuate parallel-directed fields of oblique source                                     
abs_layers = [mp.Absorber(thickness=dpml,direction=mp.X)]

wvl = 0.5              # center wavelength                                                                                                    
fcen = 1/wvl           # center frequency                                                                                                     
df = 0.05*fcen         # frequency width                                                                                                      

ng = 1.5
glass = mp.Medium(index=ng)

# out-of-plane (XZ) rotation angle of incident planewave; CCW about Y axis, 0 degrees along in-plane direction                                
theta_out = math.radians(20.3)

# in-plane (XY) rotation angle of incident planewave; CCW about Z axis, 0 degrees along +X axis                                               
theta_in = math.radians(10.7)

# k (in source medium) with correct length (plane of incidence: XY)                                                                           
k = mp.Vector3(math.cos(theta_in)*math.cos(theta_out),math.sin(theta_in)*math.cos(theta_out),math.sin(theta_out)).scale(fcen*ng)

symmetries = []
eig_parity = mp.NO_PARITY
if theta_out == 0:
  eig_parity = mp.ODD_Z
  if theta_in == 0:
    k = mp.Vector3(0,0,0)
    symmetries = mp.Mirror(mp.Y)
    eig_parity += mp.EVEN_Y

def pw_amp(k,x0):
  def _pw_amp(x):
    return cmath.exp(1j*2*math.pi*k.dot(x+x0))
  return _pw_amp

src_pt = mp.Vector3(-0.5*sx+dpml+0.3*dsub,0,0)
sources = [mp.Source(mp.GaussianSource(fcen,fwidth=df),
                     component=mp.Ez,
                     center=src_pt,
                     size=mp.Vector3(0,sy,0),
                     amp_func=pw_amp(k,src_pt))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=abs_layers,
                    k_point=k,
                    default_material=glass,
                    sources=sources,
                    symmetries=symmetries)

refl_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0)
refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))

sim.run(until_after_sources=200)

input_flux = mp.get_fluxes(refl_flux)
input_flux_data = sim.get_flux_data(refl_flux)

sim.reset_meep()

geometry = [mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub),0,0)),
            mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,0,0))]

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=abs_layers,
                    geometry=geometry,
                    k_point=k,
                    sources=sources,
                    symmetries=symmetries)

refl_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=refl_pt, size=mp.Vector3(0,sy,0)))
sim.load_minus_flux_data(refl_flux,input_flux_data)

tran_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0)
tran_flux = sim.add_flux(fcen, 0, 1, mp.FluxRegion(center=tran_pt, size=mp.Vector3(0,sy,0)))

sim.run(until_after_sources=300)

nm_r = np.floor((math.sqrt(math.pow(fcen*ng,2)-math.pow(k.z,2))-k.y)*gp)-np.ceil((-math.sqrt(math.pow(fcen*ng,2)-math.pow(k.z,2))-k.y)*gp) # \
number of reflected orders                                                                                                                    
if theta_in == 0 and theta_out == 0:
  nm_r = nm_r*0.5 # since eig_parity removes degeneracy in y-direction                                                                        
Rsum = 0
for nm in range(int(nm_r)):
  res = sim.get_eigenmode_coefficients(refl_flux, [nm+1], eig_parity=eig_parity)
  r_coeffs = res.alpha
  r_kdom = res.kdom[0]
  Rmode = (abs(r_coeffs[0,0,1])**2)/input_flux[0]
  r_angle = np.sign(r_kdom.y)*math.acos(r_kdom.x/(ng*fcen))
  print("refl:, {}, {:.2f}, {:.8f}".format(nm,math.degrees(r_angle),Rmode))
  Rsum += Rmode

nm_t = np.floor((math.sqrt(math.pow(fcen,2)-math.pow(k.z,2))-k.y)*gp)-np.ceil((-math.sqrt(math.pow(fcen,2)-math.pow(k.z,2))-k.y)*gp) # number\
 of transmitted orders                                                                                                                        
if theta_in == 0 and theta_out == 0:
  nm_t = nm_t*0.5 # since eig_parity removes degeneracy in y-direction                                                                        
Tsum = 0
for nm in range(int(nm_t)):
  res = sim.get_eigenmode_coefficients(tran_flux, [nm+1], eig_parity=eig_parity)
  t_coeffs = res.alpha
  t_kdom = res.kdom[0]
  Tmode = (abs(t_coeffs[0,0,0])**2)/input_flux[0]
  t_angle = np.sign(t_kdom.y)*math.acos(t_kdom.x/fcen)
  print("tran:, {}, {:.2f}, {:.8f}".format(nm,math.degrees(t_angle),Tmode))
  Tsum += Tmode

print("mode-coeff:, {:.6f}, {:.6f}, {:.6f}".format(Rsum,Tsum,Rsum+Tsum))

r_flux = mp.get_fluxes(refl_flux)
t_flux = mp.get_fluxes(tran_flux)
Rflux = -r_flux[0]/input_flux[0]
Tflux =  t_flux[0]/input_flux[0]
print("poynting-flux:, {:.6f}, {:.6f}, {:.6f}".format(Rflux,Tflux,Rflux+Tflux))

The output indicates that the cell is 3d which is correct.

-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 4.5 x 10 x 0.025 with resolution 40
     block, center = (-1.25,0,0)
          size (2,1e+20,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (2.25,2.25,2.25)
     block, center = (0,0,0)
          size (0.5,5,1e+20)
          axes (1,0,0), (0,1,0), (0,0,1)
          dielectric constant epsilon diagonal = (2.25,2.25,2.25)
time for set_epsilon = 0.336185 s
time for set_conductivity = 0.006742 s
time for set_conductivity = 0.00620389 s
time for set_conductivity = 0.00629306 s
time for set_conductivity = 0.00618911 s
time for set_conductivity = 0.00618219 s
time for set_conductivity = 0.0062089 s
-----------
Meep: using complex fields.
...

Note that the eig_parity is set to NO_PARITY due to the non-zero k_z. In this example, k_z is 1.0408. All modes computed by MPB via get_eigenmode_coefficients should have the same value. However, the k_z values of these modes are always 0.

..
Dominant planewave for band 4: (1.998494,-0.077596,0.000000)
...
Dominant planewave for band 26: (1.900688,0.622404,0.000000)
..

This indicates that get_eigenmode_coefficients is not correctly inferring the dimensions of the cell.

As an additional note, calling get_eigenmode separately with a non-zero k_z does yield modes with the same k_z. This is demonstrated in the following script.

import meep as mp
import math

resolution = 50        # pixels/μm                                                                                                            

dpml = 2               # PML thickness                                                                                                        
dmat = 10              # length of material                                                                                                   

sx = dpml+dmat+dpml
sy = 10

cell_size = mp.Vector3(sx,sy,0)
pml_layers = [mp.Absorber(thickness=dpml,direction=mp.X)]

wvl = 0.5              # center wavelength                                                                                                    
fcen = 1/wvl           # center frequency                                                                                                     

ng = 1.5
glass = mp.Medium(index=ng)

num_band = 6

# out-of-plane (XZ) rotation angle of incident planewave; CCW about Y axis, 0 degrees along in-plane direction                                
theta_out = math.radians(20.3)

# in-plane (XY) rotation angle of incident planewave; CCW about Z axis, 0 degrees along +X axis                                               
theta_in = math.radians(10.7)

# k (in source medium) with correct length (plane of incidence: XY)                                                                           
k = mp.Vector3(math.cos(theta_in)*math.cos(theta_out),math.sin(theta_in)*math.cos(theta_out),math.sin(theta_out)).scale(fcen*ng)

sim = mp.Simulation(resolution=resolution,
                    cell_size=cell_size,
                    boundary_layers=pml_layers,
                    k_point=k)

sim.init_sim()

mon_pt = mp.Vector3(0.5*sx-dpml-0.2*dmat,0,0)
mon_EigenmodeData = sim.get_eigenmode(fcen, mp.X, mp.Volume(center=mon_pt,size=mp.Vector3(0,sy,0)), num_band, k, verbose=True)

The output shows a KPOINT with the correct k_z value.

-----------
Initializing structure...
Working in 3D dimensions.
Computational cell is 14 x 10 x 0.02 with resolution 50
time for set_epsilon = 0.914184 s
time for set_conductivity = 0.0185812 s
time for set_conductivity = 0.0202799 s
time for set_conductivity = 0.018465 s
time for set_conductivity = 0.0202739 s
time for set_conductivity = 0.0190289 s
time for set_conductivity = 0.019552 s
-----------
Meep: using complex fields.
KPOINT: 2.76475, 0.522404, 1.04081
    near maximum in trace
...

[oskooi]

The initial portion of mpb.cpp:get_eigenmode used to set the kpoint (lines 238-242) only sets the k_y component to 0.52240. This is expected since eig_vol has size 0 x 10 x 0 (i.e., the size in the y-direction is equal to that of the cell size). k_z is not set here since the cell size is 0.025.

It seems that the discrepancy in the two results for get_eigenmode_coefficients and get_eigenmode is due to the _kpoint parameter passed to get_eigenmode during get_eigenmode_coefficients: the _kpoint is always (0,0,0) whereas when calling get_eigenmode directly it is (2.76475,0.52240,1.04081).

Thus, it seems that get_eigenmode_coefficients is calling get_eigenmode with the wrong kpoint for the case of a 2d cell with non-zero k_z.