radiated flux from $E_r$ point source at $r=0$ in cylindrical coordinates depends on cell size

[oskooi]

#2182 demonstrated that for simulations involving an $E_r$ point source at $r=0$ in cylindrical coordinates via $m=\pm 1$, it is necessary for the GaussianSource object to specify: (1) a narrow bandwidth (via a small fwidth) and (2) a large cutoff. These settings are required to mitigate the high-frequency components generated when the source abruptly turns on/off. While this approach ensures that the radiated flux from the point source converges with runtime, there seems to be a separate unexpected phenomenon in which the radiated flux varies with the size of the cell.

Here is a demonstration using the same example in #2182: a point source in vacuum enclosed by flux monitors on three sides. A schematic of the setup and the full script is below. The flux in the $+z$, $+r$, and $-z$ directions are computed for two values of the non-PML (s) and PML thicknesses (dpml). The simulation uses $m=-1$. The four reported flux values are: (1) $+z$, (2) $+r$, (3) $-z$, and (4) the total (sum of 1-3). The total flux should be the same in all test cases.

1. s = 5.0, dpml = 1.0

2.2651640186, 2.5274938264, 2.2651640186, 7.0578218636

2. s = 5.0, dpml = 2.0

2.2659956514, 2.5234538578, 2.2659956514, 7.0554451607

3. s = 10.0, dpml = 1.0

2.2774459219, 2.5238767135, 2.2774459219, 7.0787685573

Note that the four flux values in each of the three test cases are different. For comparison, these values are the same in an equivalent simulation in 2d and 3d (results not shown).

There is another anomaly: for the same cell size, the radiated flux in the $\pm z$ directions (first and third values) are different for $m=+1$ and $m=-1$. They are expected to be the same. The results above were for $m=-1$. Here are the results for $m=+1$ for the same three test cases:

4. s = 5.0, dpml = 1.0

2.2326930716, 2.5274938264, 2.2326930716, 6.9928799697

5. s = 5.0, dpml = 2.0

2.2294222475, 2.5234538578, 2.2294222475, 6.9822983528

6. s = 10.0, dpml = 1.0

2.2643590373, 2.5238767135, 2.2643590373, 7.0525947882

Why does this matter? Because the radiated flux from the $r = 0$ source depends on the cell size, it is difficult to compute an accurate value for the constant scaling factor relating the $r = 0$ and $r \neq 0$ sources in #2108.

"""
Computes the radiated flux from an Er point source at r=0 in vacuum in cylindrical coordinates.
"""

import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import meep as mp


resolution = 30  # pixels/um                                                                                      

s = 5.
dpml = 1.
cell_size = mp.Vector3(s+dpml,0,s+2*dpml)

pml_layers = [mp.PML(dpml)]

fcen = 1.

sources = [
    mp.Source(
        src=mp.GaussianSource(fcen,fwidth=0.02*fcen,cutoff=10.0),
        center=mp.Vector3(),
        component=mp.Er,
    ),
]

sim = mp.Simulation(
    resolution=resolution,
    cell_size=cell_size,
    dimensions=mp.CYLINDRICAL,
    m=-1,
    sources=sources,
    boundary_layers=pml_layers
)

flux_plus_z = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
        center=mp.Vector3(0.5*s,0,0.5*s),
        size=mp.Vector3(s,0,0)
    )
)

flux_plus_r = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
	center=mp.Vector3(s,0,0),
	size=mp.Vector3(0,0,s)
    )
)

flux_minus_z = sim.add_flux(
    fcen,
    0,
    1,
    mp.FluxRegion(
	center=mp.Vector3(0.5*s,0,-0.5*s),
	size=mp.Vector3(s,0,0),
	weight=-1.0
    )
)

sim.plot2D()
if mp.am_master():
    plt.savefig('cyl_sim.png',dpi=150,bbox_inches='tight')

sim.run(until_after_sources=0)

for t in [0, 100, 200, 400, 800]:
    sim.run(until=t)

    flux_plusz = mp.get_fluxes(flux_plus_z)[0]
    flux_plusr = mp.get_fluxes(flux_plus_r)[0]
    flux_minusz = mp.get_fluxes(flux_minus_z)[0]
    flux_tot = flux_plusz + flux_plusr + flux_minusz
    print(f"flux:, {sim.meep_time():5.1f}, {flux_plusz:.10f}, "
          f"{flux_plusr:.10f}, {flux_minusz:.10f}, {flux_tot:.10f}")

[stevengj]

I'd like to see a plot of the flux vs dpml for fixed s. And maybe also at double the resolution.

[oskooi]

I'd like to see a plot of the flux vs dpml for fixed s. And maybe also at double the resolution.

The convergence results were accidentally posted in #2108 (comment). Those results suggest that there is a likely bug in the $z$-PML implementation in cylindrical coordinates (the $r$-PML is fine).

To investigate the $z$-PML further, here is a plot of the radiated flux from an $E_r$ point dipole source at $r=z=0$ vs. $r$-PML thickness but with the $z$-PML removed. This is a cavity with metallic sidewalls in the $z$ direction. This setup (shown in the schematic below) is identical to the LDOS tutorial. The flux results are identical for the $m=+1$ and $m=-1$ simulations, as expected. (Results were generated using this script.)

However, when we re-introduce the $z$-PML and plot the radiated flux vs. $z$-PML thickness for a fixed $r$-PML thickness, the results for the $m=+1$ and $m=-1$ simulations are different (by a non-trivial amount). It is unclear which of the two sets of results is correct.

[stevengj]

Well, ±m aren't exactly equivalent simulations (they aren't just complex-conjugates) because you have a complex-ω source with only the +ω term. The results should converge to one another as you increase both z and r PML thicknesses (and maybe the resolution too), I think.

If you want to make the ±m simulations exactly the same, differing only by a complex conjugation, then one way should be to add two sources, one at +fcen and one at -fcen.

[oskooi]

The results should converge to one another as you increase both z and r PML thicknesses (and maybe the resolution too), I think.

Increasing both $z$ and $r$ PML thicknesses simultaneously also demonstrates that the $m=+1$ and $m=-1$ simulations produce different results. This is true even when doubling the resolution from 50 to 100 as shown in the two figures. The horizontal axis labelled "PML thickness" refers to the $z$ and $r$ PML thicknesses.

The fact that the results for the $m=+1$ and $m=-1$ simulations are identical when the $z$ PML is absent (shown previously) does seem to indicate a likely bug with the $z$ PML.

Note that whenever the $z$ PML is present, we found that we needed to use a narrow bandwidth source with large cutoff to ensure the fields decay away. However, when the $z$ PML is absent, the pulsed source need not have these two properties. This could perhaps be another indication that something is wrong with the $z$ PML.

[oskooi]

Fixed by #2382.