Wrong results for simulation of a simple rectangular cross grating

[Samson-Ming]

Hi all,

I hope to simulate a simple cross grating made of a cylinder array on a substrate, and I just modify the attacted example "grating2d_triangular_lattice.py" and also add some extra consideration of symmetries, but I also obtain wrong results (R and T larger than 1). I wish someone could try my code and help me find what is problem. Thank you some much.

Besides, I use an image in docker.

Below is the code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Computes the diffraction orders of a 2D binary grating with
# rectangular lattice using a rectangular supercell and verifies
# that only the diffraction orders of the actual unit cell
# produce non-zero power (up to discretization error)
import math

import numpy as np

import meep as mp

resolution = 50  # pixels/渭m

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

nSi=3.361
Si = mp.Medium(index=nSi)

wvl = 1.55  # wavelength
fcen = 1 / wvl

# rectangular supercell
sx = 0.65
sy = 0.65

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

sz = dpml + dsub + hcyl + dair + dpml

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

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

# periodic boundary conditions
k_point = mp.Vector3(0,0,0)
symmetries=[mp.Mirror(mp.X),mp.Mirror(mp.Y)]

src_pt = mp.Vector3(0, 0, -0.5 * sz + dpml+0.3*dsub)
sources = [
    mp.Source(
        src=mp.GaussianSource(fcen, fwidth=0.1 * fcen),
        size=mp.Vector3(sx, sy, 0),
        center=src_pt,
        component=mp.Ex,
    )
]

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

input_pt = mp.Vector3(0, 0, 0.5 * sz - dpml-0.2*dair)
flux_mon= sim.add_mode_monitor(
    fcen, 0, 1, mp.ModeRegion(center=input_pt, size=mp.Vector3(sx, sy, 0))
)


sim.run(until_after_sources=mp.stop_when_fields_decayed(20, mp.Ex, src_pt, 1e-6))

input_flux = mp.get_fluxes(flux_mon)

print(input_flux[0])

#######################################################################
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,
    )
]

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

geometry = substrate + cyl_grating

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

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

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


sim.run(until_after_sources=mp.stop_when_fields_decayed(20, mp.Ex, src_pt, 1e-6))

# diffraction order of unit cell (rectangular lattice)
mx = 0
my = 0

# check: for diffraction orders 
nx = mx
for ny in range(4):
    kz2 = fcen**2- (nx / sx) ** 2 - (ny / sy) ** 2
    if kz2 > 0:
        res1 = sim.get_eigenmode_coefficients(
            tran_flux, mp.DiffractedPlanewave((nx, ny, 0), mp.Vector3(0, 1, 0), 1, 0)
        )
    kz_r2=(ng*fcen)**2  - (nx / sx) ** 2 - (ny / sy) ** 2
    if kz_r2 > 0:
        res2 = sim.get_eigenmode_coefficients(
            ref_flux, mp.DiffractedPlanewave((nx, ny, 0), mp.Vector3(0, 1, 0), 1, 0)
        )
        t_coeffs = res1.alpha
        T00 = abs(t_coeffs[0, 0, 0]) ** 2/input_flux[0]
        r_coeffs = res2.alpha
        R00 = abs(r_coeffs[0, 0, 0]) ** 2/input_flux[0]

        T_flux=mp.get_fluxes(tran_flux)
        T_tot=T_flux[0]/input_flux[0]
        R_flux=mp.get_fluxes(ref_flux)
        R_tot=R_flux[0]/input_flux[0]*(-1)
        
        print(f"order:, {nx}, {ny},{R00:.5f}")
        print(f"order:, {nx}, {ny},{T00:.5f}")
        print(input_flux[0])
        print(R_flux[0])
        print(T_flux[0])
        print(R_tot)
        print(T_tot)

[Samson-Ming]

I appreciate Professor Johnson's kind concern, and I would add some more details:

(1) The meep version I adopt is an image in docker from mochen4 (https://hubgw.docker.com/r/mochen4/mpdocker) which was updated 3 years ago.

(2) A converged reference result for this diffraction grating problem by other methods (FEM, RCWA, FIT, etc) is

 R=0.0823775905
 T=0.9176244730