Irregularities in DFT flux calculation in -z direction in cylindrical coordinates

[oskooi]

There seems to be a bug in the calculation of the Poynting flux in the $-z$ direction in cylindrical coordinates. This is demonstrated using a simple test case involving a point-dipole source for $E_r$ at $r \approx 0$ for $m = -1$. The radiated flux from the point dipole is computed using two different methods: (1) LDOS and (2) DFT flux monitors. The radiated flux should be the same (up to discretization error) using both approaches.

The test seems to be valid only when method (2) involves computing the flux in the $+z$ direction. Whenever the flux is computed in the $-z$ direction, the result is incorrect. This is demonstrated in several test cases shown below.

1. Dipole in vacuum above perfect reflector ground plane. [CORRECT]

DFT flux monitor in $r$ and $+z$. PML in $r$ and $+z$.

dipole_power:, 3.541996 (LDOS), 3.537343 (flux), 0.001315 (error)

2. Dipole in vacuum.

dipole_power:, 3.691985 (LDOS), 0.343770 (2a: flux-closed), 3.686632 (2b: flux-open-plus-z), -2.999091 (2c: flux-open-minus-z)

2a. DFT flux monitor in $r$, $+z$ and $-z$. PML in $r$, $+z$, and $-z$. [INCORRECT]

2b. DFT flux monitor in $r$ and $+z$. PML in $r$, $+z$, and $-z$. [CORRECT]

note: flux is multiplied by two because it is being collected in half the necessary area.

2c. DFT flux monitor in $r$ and $-z$. PML in $r$, $+z$, and $-z$. [INCORRECT]

note: flux is multiplied by two because it is being collected in half the necessary area.

import matplotlib.pyplot as plt
import meep as mp
import numpy as np

RESOLUTION_UM = 50
PML_UM = 1.0

sr = 5.0
sz = 3.0
cell_size = mp.Vector3(sr, 0, sz)

dipole_pos = mp.Vector3(1.5 / RESOLUTION_UM, 0, -0.5 * sz + 0.5)
frequency = 1.0
dipole_component = mp.Er
sources = [
    mp.Source(
        mp.GaussianSource(frequency, fwidth=0.1 * frequency),
        component=dipole_component,
        center=dipole_pos,
    )
]

boundary_layers = [
    mp.PML(thickness=PML_UM, direction=mp.R),
    mp.PML(thickness=PML_UM, direction=mp.Z, side=mp.High)
]

sim = mp.Simulation(
    resolution=RESOLUTION_UM,
    cell_size=cell_size,
    dimensions=mp.CYLINDRICAL,
    m=1,
    sources=sources,
    boundary_layers=boundary_layers,
    force_complex_fields=True
)

flux_mon = sim.add_flux(
    frequency,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(0.5 * (sr - PML_UM), 0, 0.5 * sz - PML_UM),
        size=mp.Vector3(sr - PML_UM, 0, 0),
    ),
    mp.FluxRegion(
        center=mp.Vector3(sr - PML_UM, 0, -0.5 * sz + 0.5 * (sz - PML_UM)),
        size=mp.Vector3(0, 0, sz - PML_UM),
    )
)

fig, ax = plt.subplots()
sim.plot2D(ax=ax)
if mp.am_master():
    fig.savefig("dipole_cyl_half_space.png", dpi=150, bbox_inches="tight")

sim.run(
    mp.dft_ldos(frequency, 0, 1),
    until_after_sources=mp.stop_when_fields_decayed(
    25, dipole_component, dipole_pos, 1e-8
    )
)

if dipole_pos.x == 0:
    delta_vol = np.pi / (RESOLUTION_UM**3)
else:
    delta_vol = 2 * np.pi * dipole_pos.x / (RESOLUTION_UM**2)

dipole_power_from_ldos = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * delta_vol
dipole_power_from_flux = mp.get_fluxes(flux_mon)[0]
err = abs(dipole_power_from_ldos - dipole_power_from_flux) / dipole_power_from_flux

print(f"dipole_power:, {dipole_power_from_ldos:.6f} (LDOS), "
      f"{dipole_power_from_flux:.6f} (flux), {err:.6f} (error)")

[oskooi]

Turns out there is no bug after all in the flux calculation in $-z$. I just needed to specify weight=-1.0 in the FluxRegion object for the $-z$ flux which had been left out by mistake.

With this change, the results are as expected: method (2), independent of the orientation of the flux monitors enclosing the point dipole, produces the same flux as (1). The relative error in the two methods seems to converge to zero with increasing resolution.

dipole_power:, 3.691984 (LDOS), 3.686698 (2a: flux-closed), 3.755602 (2b: flux-open-plus-z), 3.617795 (2c: flux-open-minus-z)

flux_mon_closed = sim.add_flux(
    frequency,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(0.5 * (sr - PML_UM), 0, 0.5 * sz - PML_UM),
        size=mp.Vector3(sr - PML_UM, 0, 0),
    ),
    mp.FluxRegion(
        center=mp.Vector3(sr - PML_UM, 0, 0),
        size=mp.Vector3(0, 0, sz - 2 * PML_UM),
    ),
    mp.FluxRegion(
        center=mp.Vector3(0.5 * (sr - PML_UM), 0, -0.5 * sz + PML_UM),
        size=mp.Vector3(sr - PML_UM, 0, 0),
        weight=-1.0 # -z direction
    ),
)

flux_mon_open_plus_z = sim.add_flux(
    frequency,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(0.5 * (sr - PML_UM), 0, 0.5 * sz - PML_UM),
        size=mp.Vector3(sr - PML_UM, 0, 0),
    ),
    mp.FluxRegion(
        center=mp.Vector3(sr - PML_UM, 0, 0.5 * sz - PML_UM - 0.125 * (sz - PML_UM)),
        size=mp.Vector3(0, 0, 0.25 * (sz - PML_UM)),
    ),
)

flux_mon_open_minus_z = sim.add_flux(
    frequency,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(0.5 * (sr - PML_UM), 0, -0.5 * sz + PML_UM),
        size=mp.Vector3(sr - PML_UM, 0, 0),
        weight=-1.0 # -z direction
    ),
    mp.FluxRegion(
        center=mp.Vector3(sr - PML_UM, 0, -0.5 * sz + PML_UM + 0.125 * (sz - PML_UM)),
        size=mp.Vector3(0, 0, 0.25 * (sz - PML_UM)),
    ),
)