search

NanoComp / meep

free finite-difference time-domain (FDTD) software for electromagnetic simulations
GNU General Public License v2.0
1.26k stars 629 forks source link

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

Closed oskooi closed 6 years ago

oskooi commented 6 years ago

The kz 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 kz 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 kz. The k_point had originally been "in the plane" with kz=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 kz. In this example, kz is 1.0408. All modes computed by MPB via get_eigenmode_coefficients should have the same value. However, the kz 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 kz does yield modes with the same kz. 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 kz 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 commented 6 years ago

The initial portion of mpb.cpp:get_eigenmode used to set the kpoint (lines 238-242) only sets the ky 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). kz 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 kz.