Inconsistent result on different simulation resolutions

[radillus]

Hi, I'm simulating meta-len with pymeep 1.27.0(conda-forge). The illustration of experiment is exactly same as #2578 , the meta-len is placed at Design Region in figure, and a monitor plane gives field distribution after gaussian source passing the len.

It's highly stable when simulate with same metalen design and SAME resolution, however, when resolution changes, result becomes inconsistent even with very high resolution:

wavelength is 0.98 unit, and the size of pixels in design surface(left) is 0.05

here I demonstrate Ez's real and imag part(middle and right) under different simulation resolutions 170 136 150 135, about 150 pixel per wavelenth unit, 7 pixel per design unit(a pixel of design surface). As you can see, huge difference between different resolution makes results look unreliable. BTW I have tried using sim.run(until=..time..) and set a very long stop time and did not see any improvement.

Is this expected behavior or just I missed something?

---

here's full simulate result

https://github.com/NanoComp/meep/assets/81346885/b61fece9-3be5-4d1f-af1c-99f88e65d63d

here's code for simulation and visualization

import meep as mp
import numpy as np

DESIGN_RESOLUTION = 0.05  # μm(50nm)

def simulate(surface_refractive_mat:np.ndarray, simulation_resolution:float = 60) -> np.ndarray:
    dpml = 1.0       # PML层厚度 μm
    dsub = 0.5       # 基底厚度 μm
    dpad = 0.75      # 模型顶端与PML层距离 μm
    gh = 0.5         # 高度 μm
    lcen = 0.98      # 中心波长 μm

    gp = len(surface_refractive_mat) * DESIGN_RESOLUTION  # 周期 μm

    fcen = 1/lcen
    df = 0.2*fcen
    sx = dpml + dsub + gh + dpad + dpml
    sy = gp
    sz = gp

    pixel_len = DESIGN_RESOLUTION

    k_point = mp.Vector3(0, 0, 0)
    a_Si = mp.Medium(index=3.4)
    SiO2 = mp.Medium(index=1.5)
    pml_layers = [mp.PML(thickness=dpml, direction=mp.X)]

    src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub, 0, 0)
    mon_pt = mp.Vector3(0.5*sx-dpml-0.6*dpad, 0, 0)
    geometry = [mp.Block(
        material=SiO2,
        size=mp.Vector3(dpml+dsub, mp.inf, mp.inf),
        center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub), 0, 0),
    )]
    cell_size = mp.Vector3(sx, sy, sz)
    sources = [mp.Source(
        mp.GaussianSource(fcen, fwidth=df),
        component=mp.Ez,
        center=src_pt,
        size=mp.Vector3(0, sy, sz)
        )]
    for i,j in zip(*np.nonzero(surface_refractive_mat)):
        geometry.append(
            mp.Block(
                material=mp.Medium(index=surface_refractive_mat[i,j]),
                size=mp.Vector3(gh, pixel_len, pixel_len),
                center=mp.Vector3(
                    -0.5*sx+dpml+dsub+0.5*gh,
                    -gp / 2 + pixel_len / 2 + i * pixel_len,
                    -gp / 2 + pixel_len / 2 + j * pixel_len,
                    )
                )
            )

    sim = mp.Simulation(
        resolution=simulation_resolution,
        cell_size=cell_size,
        boundary_layers=pml_layers,
        geometry=geometry,
        k_point=k_point,
        sources=sources,
        )
    nonpml_vol = mp.Volume(center=mon_pt, size=mp.Vector3(
        y=sy,
        z=sz,
    ))
    dft_obj = sim.add_dft_fields([mp.Ex,mp.Ey,mp.Ez,mp.Hx,mp.Hy,mp.Hz], fcen, 0, 1, where=nonpml_vol)
    sim.run(until_after_sources=mp.stop_when_fields_decayed(3,mp.Ez,mon_pt,3e-4))

    Ex = sim.get_dft_array(dft_obj, mp.Ex, 0)
    Ey = sim.get_dft_array(dft_obj, mp.Ey, 0)
    Ez = sim.get_dft_array(dft_obj, mp.Ez, 0)
    Hx = sim.get_dft_array(dft_obj, mp.Hx, 0)
    Hy = sim.get_dft_array(dft_obj, mp.Hy, 0)
    Hz = sim.get_dft_array(dft_obj, mp.Hz, 0)

    field = np.empty((12,*Ex.shape))
    field[0,...] = Ex.real
    field[1,...] = Ex.imag
    field[2,...] = Ey.real
    field[3,...] = Ey.imag
    field[4,...] = Ez.real
    field[5,...] = Ez.imag
    field[6,...] = Hx.real
    field[7,...] = Hx.imag
    field[8,...] = Hy.real
    field[9,...] = Hy.imag
    field[10,...] = Hz.real
    field[11,...] = Hz.imag
    return field

# simulate under different resolution
RESOLUTION_FROM = 10
RESOLUTION_TO = 301
results = []
surface = np.random.rand(20,20) + 1

for resolution in range(RESOLUTION_FROM, RESOLUTION_TO):
    field = simulate(surface, resolution)
    results.append(field[4:6,:,:])  # Ez

# visualization
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

PLAY_SECONDS = 30

fig, axs = plt.subplots(1,3,figsize=(20,8),tight_layout=True)
axs[0].imshow(surface)
axs[0].set_title('surface')

def animate(i):
    fig.suptitle(f'resolution{i + RESOLUTION_FROM} size{len(surface)} to size{len(results[i])}')
    axs[1].clear()
    axs[1].imshow(np.real(results[i]))
    axs[1].set_title('field re')
    axs[2].clear()
    axs[2].imshow(np.imag(results[i]))
    axs[2].set_title('field im')

ani = FuncAnimation(fig, animate, frames=len(results), interval=int(1000*PLAY_SECONDS/len(results)), cache_frame_data=False)
try:
    ani.save(f'resolution_stability_test.mp4', writer='ffmpeg')
except BaseException:
    ani.save(f'resolution_stability_test.gif')
plt.close('all')

[stevengj]

Maybe you have a very sharp resonance? I would try looking at the transmission spectrum.