bug in mode coefficients with DiffractedPlanewave and non-zero k_point

[oskooi]

#2103 described a procedure to compute the mode coefficients of a triangular lattice using a normally incident source for which k_point = 0. However, it is not possible to extend this procedure to an oblique source involving a nonzero k_point. This is because, as shown in the code snippet below, there is a bug in fields::get_eigenmode of src/mpb.cpp which occurs whenever the cell volume is modified internally from whatever the user had originally specified due to: (1) rounding of the cell dimensions to an integer number of pixels and/or (2) applying a symmetry operator which changes the cell dimensions (and may also involve padding the cell by a voxel).

Any change to the cell dimensions (v) in the periodic directions causes the check on line 339 against the dimensions of the DFT monitor (eig_vol) which extends the entire length of the unmodified cell in the periodic directions to fail: float(eig_vol.in_direction(dd)) == float(v.in_direction(dd). As a result, the variable kpoint is 0 rather than the k_point of the Simulation object. The wrong modes are therefore used to compute the mode coefficients.

meep/src/mpb.cpp

Lines 334 to 342 in 619c879

	// if the mode region extends over the full computational grid and we are bloch-periodic
	// in any direction, set the corresponding component of the eigenmode initial-guess
	// k-vector to be the (real part of the) bloch vector in that direction.
	vec kpoint(_kpoint);
	LOOP_OVER_DIRECTIONS(v.dim, dd) {
	if (dd != d && float(eig_vol.in_direction(dd)) == float(v.in_direction(dd)))
	if (boundaries[High][dd] == Periodic && boundaries[Low][dd] == Periodic)
	kpoint.set_direction(dd, real(k[dd]));
	}

To demonstrate this bug, we set up a test involving computing the diffraction orders of a 2D binary grating with triangular lattice. The test involves verifying that the sum of the power in all the reflected and transmitted diffraction orders (computing using mode decomposition) is equivalent to the total reflected and transmitted Poynting flux, respectively. There are two test cases: the incident planewave source is either (1) normal or (2) oblique. Results for (1) are correct but not (2).

The data shows the output of running the script below. There are two rows for the reflected (refl) and transmitted (tran) orders. The three numerical columns in each row are the (1) Poynting flux (flux), (2) sum of the flux of all the diffraction orders (orders), and (3) the relative error (error). The error is small for (1) but not (2).

1. theta = 0

refl:, 0.023371 (flux), 0.022366 (orders), 0.042999 (error)
tran:, 0.976509 (flux), 0.976532 (orders), 0.000024 (error)

2. theta = math.radians(50.0)

refl:, 0.021790 (flux), 0.009713 (orders), 0.554259 (error)
tran:, 0.977619 (flux), 0.056392 (orders), 0.942316 (error)

# Computes the diffraction orders of a 2D binary grating                                                                                                      
# with triangular lattice using a planewave at                                                                                                                
# oblique incidence.                                                                                                                                          

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


resolution = 80

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

wvl = 0.5    # wavelength                                                                                                                                     
fcen = 1/wvl

dpml = 1.0   # PML thickness                                                                                                                                  
dsub = 2.0   # substrate thickness                                                                                                                            
dair = 2.0   # air padding                                                                                                                                    
rcyl = 0.1   # cylinder radius                                                                                                                                
hcyl = 0.3   # cylinder height                                                                                                                                

# periodicity of triangular lattice                                                                                                                           
a = 0.6

# rectangular supercell                                                                                                                                       
sx = a
sy = a*np.sqrt(3)

sz = dpml+dsub+hcyl+dair+dpml

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

boundary_layers = [mp.PML(thickness=dpml,direction=mp.Z)]

# plane of incidence is yz                                                                                                                                    
# 0° is +z with CCW rotation about x                                                                                                                          
theta = math.radians(50.0)

if theta == 0:
    k = mp.Vector3()
else:
    k = mp.Vector3(0,0,fcen).rotate(mp.Vector3(1,0,0),theta)

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,0,-0.5*sz+dpml)
src_cmpt = mp.Ey # S_pol: Ex / P_pol: Ey                                                                                                                      
sources = [mp.Source(src=mp.GaussianSource(fcen,fwidth=0.1*fcen),
                     size=mp.Vector3(sx,sy,0),
                     center=src_pt,
                     component=src_cmpt,
                     amp_func=pw_amp(k,src_pt))]

symmetries = [mp.Mirror(direction=mp.X, phase=-1 if src_cmpt==mp.Ex else +1)]

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

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

stop_cond = mp.stop_when_fields_decayed(25,src_cmpt,src_pt,1e-6)
sim.run(until_after_sources=stop_cond)

input_flux = mp.get_fluxes(refl_flux)[0]
input_flux_data = sim.get_flux_data(refl_flux)

sim.reset_meep()

substrate = [mp.Block(size=mp.Vector3(mp.inf,mp.inf,dpml+dsub),
                      center=mp.Vector3(0,0,-0.5*sz+0.5*(dpml+dsub)),
                      material=glass)]

grating = [mp.Cylinder(center=mp.Vector3(0,0,-0.5*sz+dpml+dsub+0.5*hcyl),
                       radius=rcyl,
                       height=hcyl,
                       material=glass),
           mp.Cylinder(center=mp.Vector3(0.5*sx,0.5*sy,-0.5*sz+dpml+dsub+0.5*hcyl),
                       radius=rcyl,
                       height=hcyl,
                       material=glass),
           mp.Cylinder(center=mp.Vector3(-0.5*sx,0.5*sy,-0.5*sz+dpml+dsub+0.5*hcyl),
                       radius=rcyl,
                       height=hcyl,
                       material=glass),
           mp.Cylinder(center=mp.Vector3(0.5*sx,-0.5*sy,-0.5*sz+dpml+dsub+0.5*hcyl),
                       radius=rcyl,
                       height=hcyl,
                       material=glass),
           mp.Cylinder(center=mp.Vector3(-0.5*sx,-0.5*sy,-0.5*sz+dpml+dsub+0.5*hcyl),
                       radius=rcyl,
                       height=hcyl,
                       material=glass)]

geometry = substrate + grating

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

refl_flux = sim.add_mode_monitor(fcen,
                                 0,
                                 1,
                                 mp.ModeRegion(center=refl_pt,
                                               size=mp.Vector3(sx,sy,0)))

sim.load_minus_flux_data(refl_flux, input_flux_data)

tran_pt = mp.Vector3(0,0,0.5*sz-dpml)
tran_flux = sim.add_mode_monitor(fcen,
                                 0,
                                 1,
                                 mp.ModeRegion(center=tran_pt,
                                               size=mp.Vector3(sx,sy,0)))

sim.run(until_after_sources=stop_cond)

# loop over the orders of the supercell                                                                                                                       
# and sum the flux in all the reflected                                                                                                                       
# and transmitted orders of the unit cell                                                                                                                     
Rsum = 0
Tsum = 0
m = 5
tol = 1e-6
for nx in range(-m,m+1):
    for ny in range(-m,m+1):
        kz2 = (ng*fcen)**2-(k.x+nx/sx)**2-(k.y+ny/sy)**2
        if kz2 > 0:
            Rpol = 0
            for S_pol in [True, False]:
                res = sim.get_eigenmode_coefficients(refl_flux,
                                                     mp.DiffractedPlanewave((nx,ny,0),
                                                                            mp.Vector3(0,1,0),
                                                                            1 if S_pol else 0,
                                                                            0 if S_pol else 1))

                coeffs = res.alpha
                refl = abs(coeffs[0,0,1])**2 / input_flux

                if refl > tol:
                    # check if my is an integer to determine whether "real" order                                                                             
                    if (nx + ny) % 2 == 0:
                        Rpol += refl
                    else:
                        print("WARNING: artificial order contains nonzero power")
                    print("refl:, {}, {:2d}, {:2d}, {:.5f}".format('S' if S_pol else 'P',nx,ny,refl))

            Rsum += Rpol

        kz2 = fcen**2-(k.x+nx/sx)**2-(k.y+ny/sy)**2
        if kz2 > 0:
            Tpol = 0
            for S_pol in [True, False]:
                res = sim.get_eigenmode_coefficients(tran_flux,
                                                     mp.DiffractedPlanewave((nx,ny,0),
                                                                            mp.Vector3(0,1,0),
                                                                            1 if S_pol else 0,
                                                                            0 if S_pol else 1))
                coeffs = res.alpha
                tran = abs(coeffs[0,0,0])**2 / input_flux

                if tran > tol:
                    # check if my is an integer to determine whether "real" order                                                                             
                    if (nx + ny) % 2 == 0:
                        Tpol += tran
                    else:
                        print("WARNING: artificial order contains nonzero power")
                    print("tran:, {}, {:2d}, {:2d}, {:.5f}".format('S' if S_pol else 'P',nx,ny,tran))

            Tsum += Tpol


Rflux = -mp.get_fluxes(refl_flux)[0] / input_flux
err = abs(Rflux-Rsum)/Rflux
print("refl:, {:.6f} (flux), {:.6f} (orders), {:.6f} (error)".format(Rflux,Rsum,err))

Tflux = mp.get_fluxes(tran_flux)[0] / input_flux
err = abs(Tflux-Tsum)/Tflux
print("tran:, {:.6f} (flux), {:.6f} (orders), {:.6f} (error)".format(Tflux,Tsum,err))