Differences between meep.adjoint.EigenmodeCoefficient and get_eigenmode_coefficients

[rerdo]

Dear all,

I have noticed a slight discrepancy when using get_eigenmode_coefficients and meep.adjoint.EigenmodeCoefficient to calculate the scattering terms, mainly in the reflectance. I think that the only difference between these objects is that meep.adjoint.EigenmodeCoefficient was written to be passed to the optimization method of mpa.OptimizationProblem. However, when comparing the results of the optimization, there is a slight difference from the output of get_eigenmode_coefficients computed from the simulation. I wrote a small test code showing the difference.

Am I comparing similar things here?

Thank you.

"""Test structure to compare the methods:
`mpa.EigenmodeCoefficient` and `sim.get_eigenmode_coefficients`.
"""

import meep as mp
import meep.adjoint as mpa
import numpy as np

w = 1.0
l = 1.0
pad = 1.0
dpml = 1.0

resolution = 20

core = mp.Medium(index=2)

sx = dpml + l + w + l + dpml
sy = dpml + pad + w + pad + dpml
cell_size = mp.Vector3(sx, sy, 0)

pml_layers = [mp.PML(thickness=dpml)]

eig_parity = mp.ODD_Y + mp.EVEN_Z
symmetries=[mp.Mirror(mp.Y, phase=-1)]
src_pt = mp.Vector3(-0.5 * sx + dpml)
src_size = mp.Vector3(y=sy)

wvl = 1.55
fcen = 1 / wvl
df = 0.01

refl_pt = mp.Vector3(-0.5 * sx + dpml + 0.5 * l)
tran_pt = mp.Vector3(0.5 * sx - dpml - 0.5 * l)

stop_cond = mp.stop_when_fields_decayed(50, mp.Ey, refl_pt, 1e-8)

sources = [
            mp.EigenModeSource(
                src=mp.GaussianSource(fcen, fwidth=df),
                size=src_size,
                center=src_pt,
                eig_band=1,
                eig_parity=eig_parity,
            ),
        ]


def straight_waveguide():

    geometry = [
            mp.Block(
                size=mp.Vector3(mp.inf, w, mp.inf),
                center=mp.Vector3(),
                material=core,
            )
        ]

    sim = mp.Simulation(
            resolution=resolution,
            default_material=mp.air,
            cell_size=cell_size,
            sources=sources,
            geometry=geometry,
            boundary_layers=pml_layers,
            symmetries=symmetries,
            k_point=mp.Vector3(),
        )

    refl_mon = sim.add_mode_monitor(
        fcen,
        df,
        1,
        mp.ModeRegion(center=refl_pt,
                      size=mp.Vector3(y=sy)),
        yee_grid=True,
    )

    sim.run(until_after_sources=stop_cond)

    res = sim.get_eigenmode_coefficients(
        refl_mon,
        [1],
        eig_parity=eig_parity,
    )

    coeffs = res.alpha
    input_flux = np.abs(coeffs[0, :, 0]) ** 2
    input_flux_data = sim.get_flux_data(refl_mon)

    return input_flux, input_flux_data


def waveguide_crossing(method,
                       input_flux,
                       input_flux_data):

    geometry = [
            mp.Block(
                size=mp.Vector3(mp.inf, w, mp.inf),
                center=mp.Vector3(),
                material=core,
            ),
            mp.Block(
                size=mp.Vector3(w, mp.inf, mp.inf),
                center=mp.Vector3(),
                material=core,
            ),
        ]
    
    sim = mp.Simulation(
            resolution=resolution,
            default_material=mp.air,
            cell_size=cell_size,
            sources=sources,
            geometry=geometry,
            boundary_layers=pml_layers,
            symmetries=symmetries,
            k_point=mp.Vector3(),
        )
    
    if method == "sim":

        refl_mon = sim.add_mode_monitor(
            fcen,
            df,
            1,
            mp.ModeRegion(center=refl_pt,
                          size=mp.Vector3(y=sy)),
            yee_grid=True,
        )

        tran_mon = sim.add_mode_monitor(
            fcen,
            df,
            1,
            mp.ModeRegion(center=tran_pt,
                          size=mp.Vector3(y=sy)),
            yee_grid=True,
        )

        sim.load_minus_flux_data(refl_mon, input_flux_data)

        sim.run(until_after_sources=stop_cond)

        refl_res = sim.get_eigenmode_coefficients(
            refl_mon,
            [1],
            eig_parity=eig_parity,
        )

        tran_res = sim.get_eigenmode_coefficients(
            tran_mon,
            [1],
            eig_parity=eig_parity,
            )

        refl_coeff = refl_res.alpha
        tran_coeff = tran_res.alpha

        R = np.abs(refl_coeff[0, :, 1]) ** 2 / input_flux
        T = np.abs(tran_coeff[0, :, 0]) ** 2 / input_flux
    
        return R, T
    
    elif method == "mpa":

        monitor_refl = mpa.EigenmodeCoefficient(
        sim,
        mp.Volume(
            center=refl_pt,
            size=mp.Vector3(y=sy),
        ),
        1,
        forward=False,
        eig_parity=eig_parity,
        subtracted_dft_fields=input_flux_data,
        )
       
        monitor_tran = mpa.EigenmodeCoefficient(
            sim,
            mp.Volume(
                center=tran_pt,
                size=mp.Vector3(y=sy),
            ),
            1,
            eig_parity=eig_parity,
            )

        monitor_refl.register_monitors([fcen])
        monitor_tran.register_monitors([fcen])

        sim.run(until_after_sources=stop_cond)

        r = monitor_refl()[0]
        t = monitor_tran()[0]

        R = np.power(np.abs(r), 2) / input_flux
        T = np.power(np.abs(t), 2) / input_flux

        return R, T


if __name__ == "__main__":

    # get normalization
    input_flux, input_flux_data = straight_waveguide()

    # run using sim
    R_sim, T_sim = waveguide_crossing("sim", 
                             input_flux,
                             input_flux_data,
                             )

    # run using mpa
    R_mpa, T_mpa = waveguide_crossing("mpa", 
                             input_flux,
                             input_flux_data,
                             )
    
    print(f"Reflectance using simulation method: {R_sim}")
    print(f"Transmitance using simulation method: {T_sim}")
    
    print(f"Reflectance using mpa: {R_mpa}")
    print(f"Transmitance using mpa: {T_mpa}")

Reflectance using simulation method: [0.00260735]
Transmitance using simulation method: [0.77770488]
Reflectance using mpa: [0.00280257]
Transmitance using mpa: [0.77770486]

[oskooi]

It seems the problem is related to the reflectance calculation since the transmittance values are nearly identical.

[rerdo]

When I increased the resolution to 100, the difference in reflectances decreased when comparing both methods. I obtained 0.00241624 using "sim" and 0.00242317 using "mpa." Although the difference may seem very little, I think we can still conclude they are different. It's bothering me as I try to understand why I'm getting different numbers when I compare the results obtained for topology optimization and the results I obtain when I just read the optimized geometry in a meep.MaterialGrid and run the simulation in a different code.

For topology optimization, I use a resolution of 50 for the simulation and 100 for the design region. When I run the optimized geometry in a different code, I use a resolution of 100 for the simulation, and this gives me a difference of around 1 %.