Understand FluxRegion

[Camaltra]

Hey Meep Community!

Intro

I'm currently trying to understand how FluxRegion work, especially the one that are been wrapped into Mode Monitor. I'll show some experiments I've made in order to understand, and then ask some question around them.

Experiments

I've made 3 different experiences in order to try to get and idea of the computation. In all the experience, the main objective will tend to compare the flux computed between FluxRegion and Mode Monitor.

Waveguide

Following the picture, here is my experiment. On the blue monitor there is actually 2 monitor, one is FluxRegion, the other one is a Mode Monitor. I extract the value of the flux in the following way, and give at each time the results associated. Smth to note, the simulation ran for 50 timestamps, on only one frequency (1 / 1.55). All the code of the simulation will be link at the end.

res = sim.get_eigenmode_coefficients(mode_monitor, [1])
coeff = res.alpha
mode_flux = np.abs(coeff[0, :, 0])
mode_squared_flux = mode_flux ** 2

print(f'Flux: {mode_flux}')
print(f'Squared Flux: {mode_squared_flux}')

>>> Flux: [4.03785533]
>>> Squared Flux: [16.3042757]

assert len(mp.get_flux_freqs(in_flux_monitor)) == 1
assert 1/WAVELENGTH == mp.get_flux_freqs(in_flux_monitor)[-1]

fluxes = mp.get_fluxes(in_flux_monitor)[-1]
print(f'Fluxes: {fluxes}')

>>> Fluxes: 15.508946127434934

Waveguide + Box

Following the picture, here is my second experiment. Same as before about the monitor, timestamps and frequencies. Here is the results, obtain in the same fashion that the previous example

>>> Flux: [8.44064548]
>>> Squared Flux: [71.2444961]
>>> Fluxes: 68.44586640859458

Waveguide + Box Full Simulation

A little tweak here is I last the simulation much longer of the Waveguide + Box in order to get some reflection. And here I get some pretty weird results I guess. I might understand this as the reflection came across monitors multiple time. But still I not sure of how I can understand this results.

>>> Flux: [18.17772316]
>>> Squared Flux: [330.42961913]
>>> Fluxes: 273.67158562865234

Question

So here is a list of my question in order to understand:

• Why there is a huge difference between the Flux computed thought Mode Monitor and the Flux Region? Especially even with and without the squared flux on the mode monitor:

  • If I sum all the flux found for each mode from 1 to 10, I still got something way lower.
  • If we compare with the squared value of the Mode Monitor, it it much larger with not make that much sense.

• When adding reflection, the Flux increase. As I think about it, it seems that the fields came across multiple time the monitor and can explain this phenomena.

• To get the forward only field that pass through the monitor, if I understand well, I just need to get it that ways?

mode_flux = np.abs(coeff[0, :, 0])

And for the backward

mode_flux = np.abs(coeff[0, :, 1])

All these ask and work is because I'm working on the Adjoint Method solver and try to put some monitors and so in order to get the best loss function.

Code

To note, the design region is 6, the waveguide width is 0.5, and the box is 2x2. The simulation happen only in 2D.

Wagevuide

def build_straight_waveguide(params: MirrorParams, plot: bool = False) -> float:
    dpml_x = 1.0
    dpml_y = 1.0
    dpml_z = 1.0

    sx = int(dpml_x * 2 + params.region_size_x + params.wg_length * 2)
    sy = int(dpml_y * 2 + params.region_size_y + DAIR)
    sz = 0 if not params.region_size_z else dpml_z * 2 + params.region_size_z + DAIR

    geometry = [
        mp.Block(
            size=mp.Vector3(mp.inf, params.wg_width, params.region_size_z),
            center=mp.Vector3(
                0,
                0,
                0,
            ),
            material=Si,
        ),
    ]

    sources = [
        mp.EigenModeSource(
            src=mp.GaussianSource(FCEN, fwidth=DF),
            center=mp.Vector3(-sx / 2 + dpml_x + 0.3, params.wg_center),
            size=mp.Vector3(0, params.wg_width * 2),
            eig_band=1,
        )
    ]

    boundary_layer = [mp.PML(dpml_x, direction=mp.X), mp.PML(dpml_y, direction=mp.Y)]
    if sz != 0:
        boundary_layer += [mp.PML(dpml_z, direction=mp.Z)]

    sim = mp.Simulation(
        resolution=RESOLUTION,
        cell_size=mp.Vector3(sx, sy, sz),
        boundary_layers=boundary_layer,
        geometry=geometry,
        sources=sources,
        k_point=mp.Vector3(),
    )

    mode_monitor = sim.add_mode_monitor(
        [1 / WAVELENGTH],
        mp.FluxRegion(
            center=mp.Vector3(-sx / 2 + dpml_x + 0.9),
            size=mp.Vector3(0, params.wg_width * 2, 0),
        ),
        yee_grid=True,
        decimation_factor=0,
    )

    in_region = mp.FluxRegion(
        center=mp.Vector3(-sx / 2 + dpml_x + 0.9),
        size=mp.Vector3(0, params.wg_width * 2, 0)
    )
    in_flux_monitor = sim.add_flux(1/WAVELENGTH, 0, 1, in_region)

    # sim.run(until_after_sources=mp.stop_when_fields_decayed(50.0, mp.Ez, refl_pt, 1e-6))
    sim.run(until=50)

    res = sim.get_eigenmode_coefficients(mode_monitor, [1])
    coeff = res.alpha
    mode_flux = np.abs(coeff[0, :, 0])
    mode_squared_flux = mode_flux ** 2

    print(f'Flux: {mode_flux}')
    print(f'Squared Flux: {mode_squared_flux}')

    assert len(mp.get_flux_freqs(in_flux_monitor)) == 1
    assert 1/WAVELENGTH == mp.get_flux_freqs(in_flux_monitor)[-1]

    fluxes = mp.get_fluxes(in_flux_monitor)[-1]
    print(f'Fluxes: {fluxes}')

    sim.plot2D(fields=mp.Ez)
    plt.show()

Wagevuide + Box

def build_mirror_test(params):
    dpml_x = 1.0
    dpml_y = 1.0
    dpml_z = 1.0

    sx = int(dpml_x * 2 + params.region_size_x + params.wg_length * 2)
    sy = int(dpml_y * 2 + params.region_size_y + DAIR)
    sz = 0 if not params.region_size_z else dpml_z * 2 + params.region_size_z + DAIR

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

    boundary_layer = [mp.PML(dpml_x, direction=mp.X), mp.PML(dpml_y, direction=mp.Y)]
    if sz != 0:
        boundary_layer += [mp.PML(dpml_z, direction=mp.Z)]

    geometries = [
        mp.Block(
            material=Si,
            center=mp.Vector3(-params.region_size_x / 2 - params.wg_length / 2, params.wg_center),
            size=mp.Vector3(params.wg_length, params.wg_width),
        ),
        mp.Block(
            material=Si,
            center=mp.Vector3(),
            size=mp.Vector3(params.region_size_x, params.region_size_y, params.region_size_z)
        )
    ]

    sources = [
        mp.EigenModeSource(
            src=mp.GaussianSource(FCEN, fwidth=DF),
            center=mp.Vector3(-sx / 2 + dpml_x + 0.3, params.wg_center),
            size=mp.Vector3(0, params.wg_width * 2),
            eig_band=1,
        )
    ]

    sim = mp.Simulation(
        resolution=RESOLUTION,
        cell_size=cell_size,
        boundary_layers=boundary_layer,
        geometry=geometries,
        sources=sources,
        default_material=SiO2,
        force_complex_fields=True,
    )

    mode_monitor = sim.add_mode_monitor(
        [1 / WAVELENGTH],
        mp.FluxRegion(
            center=mp.Vector3(-sx / 2 + dpml_x + 0.9),
            size=mp.Vector3(0, params.wg_width * 2, 0),
        ),
        yee_grid=True,
        decimation_factor=0,
    )

    in_region = mp.FluxRegion(
        center=mp.Vector3(-sx / 2 + dpml_x + 0.9),
        size=mp.Vector3(0, params.wg_width * 2, 0)
    )
    in_flux_monitor = sim.add_flux(1/WAVELENGTH, 0, 1, in_region)

    # sim.run(until_after_sources=mp.stop_when_fields_decayed(50.0, mp.Ez, refl_pt, 1e-6))
    sim.run(until=200)

    res = sim.get_eigenmode_coefficients(mode_monitor, [1])
    coeff = res.alpha
    mode_flux = np.abs(coeff[0, :, 0])
    mode_squared_flux = mode_flux ** 2

    print(f'Flux: {mode_flux}')
    print(f'Squared Flux: {mode_squared_flux}')

    assert len(mp.get_flux_freqs(in_flux_monitor)) == 1
    assert 1/WAVELENGTH == mp.get_flux_freqs(in_flux_monitor)[-1]

    fluxes = mp.get_fluxes(in_flux_monitor)[-1]
    print(f'Fluxes: {fluxes}')

    sim.plot2D(fields=mp.Ez)
    plt.show()

[stevengj]

The flux monitor includes all the power flux (the Poynting flux) through that line segment, whereas the mode monitor only includes the flux of the modes you specify. Unless there is reflected power, the former (total flux) should generally be larger than the latter (the flux of one or more modes).

(In a situation with reflections, the total flux may be smaller, since it includes subtracted terms from the reflected waves, whereas the mode monitor can separate out modes going in different directions.)

In your case, it looks like you are reporting a total flux that is smaller than your mode flux. I'm guessing that this is because you simply aren't running the simulation for long enough — the Fourier transform hasn't converged yet. (I would also make the flux regions longer, especially for the mode calculation … you want to be sure to fully capture the exponential tails of the modes.)

[Camaltra]

Hi Steven! Thank for you reply!

I did try with a much longer waveguide, increase the simulation time, so that I don't get any reflected waves, and made sure that all the fields went through monitors as you proposed. I also increase the size of both monitor as sy - 2 * dpml to catch any possible power.
It end up that I keep getting a higher results through Mode Monitor than through Fluxes Monitor , which is a bit strange... Maybe I'm doing something wrong within my code, but it seems good to me.

To get more context, as I mentioned, I work along the Adjoint Method, and try these two different approach to measure the output power of the source on a single and empty waveguide, as I saw them as well in the documentation and in some piece of code found on github.

As I understand, we should get a higher Total Fluxes than the single mode one, but also really close as the simulation in a waveguide should mostly support the TE mode.

I follow the Ez output picture, the output, and the code

sim.run(until_after_sources = mp.stop_when_dft_decayed(1e-5))

res = sim.get_eigenmode_coefficients(mode_monitor, [1])
coeff = res.alpha
mode_flux = np.abs(coeff[0, :, 0])
mode_squared_flux = mode_flux ** 2

print(f'Mode Flux: {mode_flux}')
print(f'Mode Squared Flux: {mode_squared_flux}')

>>> Mode Flux: [10.34006895]
>>> Mode Squared Flux: [106.91702595]

assert len(mp.get_flux_freqs(in_flux_monitor)) == 1
assert 1/WAVELENGTH == mp.get_flux_freqs(in_flux_monitor)[-1]

fluxes = mp.get_fluxes(in_flux_monitor)[-1]
print(f'Fluxes: {fluxes}')

>>> Fluxes: 102.5597203845114

Full Code

WAVELENGTH = 1.55
DF = 0.1

def build_straight_waveguide(params: MirrorParams, plot: bool = False) -> float:
    dpml_x = 1.0
    dpml_y = 1.0
    dpml_z = 1.0

    sx = int(dpml_x * 2 + params.region_size_x + params.wg_length * 2) + 20
    sy = int(dpml_y * 2 + params.region_size_y + DAIR)
    sz = 0 if not params.region_size_z else dpml_z * 2 + params.region_size_z + DAIR

    geometry = [
        mp.Block(
            size=mp.Vector3(mp.inf, params.wg_width, params.region_size_z),
            center=mp.Vector3(
                0,
                0,
                0,
            ),
            material=Si,
        ),
    ]

    sources = [
        mp.EigenModeSource(
            src=mp.GaussianSource(FCEN, fwidth=DF),
            center=mp.Vector3(-sx / 2 + dpml_x + 0.3, params.wg_center),
            size=mp.Vector3(0, params.wg_width * 2),
            eig_band=1,
        )
    ]

    boundary_layer = [mp.PML(dpml_x, direction=mp.X), mp.PML(dpml_y, direction=mp.Y)]
    if sz != 0:
        boundary_layer += [mp.PML(dpml_z, direction=mp.Z)]

    sim = mp.Simulation(
        resolution=RESOLUTION,
        cell_size=mp.Vector3(sx, sy, sz),
        boundary_layers=boundary_layer,
        geometry=geometry,
        sources=sources,
        k_point=mp.Vector3(),
    )

    mode_monitor = sim.add_mode_monitor(
        [1 / WAVELENGTH],
        mp.FluxRegion(
            center=mp.Vector3(-sx / 2 + dpml_x + 0.9),
            size=mp.Vector3(0, sy - 2 * dpml_y, 0),
        ),
        yee_grid=True,
        decimation_factor=0,
    )

    in_region = mp.FluxRegion(
        center=mp.Vector3(-sx / 2 + dpml_x + 0.9),
        size=mp.Vector3(0, sy - 2 * dpml_y, 0)
    )
    in_flux_monitor = sim.add_flux(1/WAVELENGTH, 0, 1, in_region)


    sim.run(until_after_sources = mp.stop_when_dft_decayed(1e-5))
    # sim.run(until=150)

    res = sim.get_eigenmode_coefficients(mode_monitor, [1])
    coeff = res.alpha
    mode_flux = np.abs(coeff[0, :, 0])
    mode_squared_flux = mode_flux ** 2

    print(f'Flux: {mode_flux}')
    print(f'Squared Flux: {mode_squared_flux}')

    assert len(mp.get_flux_freqs(in_flux_monitor)) == 1
    assert 1/WAVELENGTH == mp.get_flux_freqs(in_flux_monitor)[-1]

    fluxes = mp.get_fluxes(in_flux_monitor)[-1]
    print(f'Fluxes: {fluxes}')

    sim.plot2D(fields=mp.Ez)
    plt.show()

Again thanks for your reply!