bug in z-PML in cylindrical coordinates at r=0 for m=0, ±1

[oskooi]

Based on the discussion in #2148 (comment).

In cylindrical coordinates, for a simulation with $m=\pm 1$ in vacuum with PMLs on all sides and an $E_r$ source at $r=0$, the $E_r$ fields at $r=0$ are not absorbed by the PMLs in the $z$ direction. The fields therefore bounce around indefinitely and the total flux emitted by the source measured using a closed three-sided box of DFT monitors never converges with run time.

More details are provided in this gist.

[oskooi]

meep/src/dft.cpp

Lines 291 to 293 in 2aa9164

	complex<realnum> fc(f[0], f[1]);
	for (int i = 0; i < Nomega; ++i)
	dft[Nomega * idx_dft + i] += dft_phase[i] * fc;

Use IVEC_LOOP_LOC to identify location of fields and print out c, f[0], and f[1] and simulation time.

Make sure that one field value is not incorrect

[oskooi]

Here is a plot of the (magnitude of the complex) time-domain fields $E_\phi$ and $E_r$ versus simulation time at points near $r=0$ and $r \neq 0$. The $E_\phi$ (E_p) and $E_r$ (E_r) fields at points closest to $r=0$ do not converge to zero. This indicates a bug. The time-domain fields everywhere else in the cell do converge to zero, as expected.

1. Field monitor (blue dot) at (r,z) = (~0,2.5).

2. Field monitor (blue dot) at (r,z) = (~2.5,2.5).

For comparison, the same computation for fields at $r \neq 0$ shows that these fields do converge to zero.

diff --git a/src/dft.cpp b/src/dft.cpp
index 802687df..2e2c2b95 100644
--- a/src/dft.cpp
+++ b/src/dft.cpp
@@ -287,6 +287,12 @@ void dft_chunk::update_dft(double time) {
       for (int cmp = 0; cmp < numcmp; ++cmp)
         f[cmp] = w * fc->f[c][cmp][idx];
 
+    if (int(time*100) % 25 == 0) {
+          IVEC_LOOP_LOC(fc->gv, loc);
+          master_printf("dft:, %6.2f, %s, (%0.2f,%0.2f), %+0.12f%+0.12fj\n",
+                        time,component_name(c),loc.r(),loc.z(),f[0],f[1]);
+    }
+
     if (numcmp == 2) {
       complex<realnum> fc(f[0], f[1]);
       for (int i = 0; i < Nomega; ++i)

import meep as mp
import numpy as np


resolution = 30  # pixels/um                                                                 

s = 5.
dpml_r = 4.
dpml_z = 2.
cell_size = mp.Vector3(s+dpml_r,0,s+2*dpml_z)

pml_layers = [mp.PML(dpml_r,direction=mp.R),
              mp.PML(dpml_z,direction=mp.Z),]

fcen = 1.

sources = [mp.Source(
    src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
    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.01,0,0.5*s),
        size=mp.Vector3(0.02,0,0)))

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

    flux_plusz = mp.get_fluxes(flux_plus_z)[0]
    print(f"flux:, {sim.meep_time()}, {flux_plusz:.6f}")

[stevengj]

I suspect this has nothing to do with the PML. In fact, you should be able to remove the PML in the z direction completely since the radial absorbers should still cause the field to decay.

It seems like there is a problem with timestepping E_φ?

[stevengj]

It's a bit confusing to me that you show the fields going to a constant and somehow that's affecting the DTFT at a nonzero frequency? Can you look at the time dependence of E_φ more carefully?

[oskooi]

I suspect this has nothing to do with the PML. In fact, you should be able to remove the PML in the z direction completely since the radial absorbers should still cause the field to decay.

Removing the PML in the $z$ direction (but leaving the $r$-PML) shows the same behavior as before: the $E_\phi$ is not decaying away to zero at $r\approx0$.

It's a bit confusing to me that you show the fields going to a constant and somehow that's affecting the DTFT at a nonzero frequency? Can you look at the time dependence of E_φ more carefully?

In this example, the field oscillations are apparent unlike the previous example with PML in $r$ and $z$.

[oskooi]

Some additional information related to #2076 and #2082 which also involved an $E_r$ source at $r=0$. The reason this bug did not become apparent in Tutorial/Planar Cavity with Lossless Metallic Walls and the related unit test in test_ldos.py is because of the particular choice for the termination condition:

    sim.run(mp.dft_ldos(fcen,0,1),
            until_after_sources=mp.stop_when_fields_decayed(20,
                                                            mp.Er,
                                                            mp.Vector3(),
                                                            1e-6))

Here, it turns out that using decay_by=1e-6 does cause the simulation to terminate but using decay_by=1e-10 does not. Also, replacing the $E_r$ source with $E_p$ causes the simulation to never terminate regardless of the choice of decay_by.

[stevengj]

Would be good to review some papers on cylindrical FDTD and check that our equations match the ones in the literature for m = ±1.

[oskooi]

Would be good to review some papers on cylindrical FDTD and check that our equations match the ones in the literature for m = ±1.

For reference in this discussion, the (non-discretized) Maxwell's equations in cylindrical coordinates implemented in an early version of Meep including a figure showing the Yee cylindrical lattice are provided in Section 1.2 and Figure 1.1 of this March 2005 document by David Roundy: http://ab-initio.mit.edu/~meep/meep.pdf.

While that document does not contain any discussion of the special case of $r=0$ and $m=\pm 1$, the cylindrical step-curl equations which were reimplemented by @stevengj in c12ff75 in June 2008 including the special case of $r=0$ are shown in the comments of step_db.cpp:

meep/src/step_db.cpp

Lines 281 to 282 in 63cd51f

	// deal with annoying r=0 boundary conditions for m=0 and m=1
	if (gv.dim == Dcyl && gv.origin_r() == 0.0) DOCMP {

meep/src/step_db.cpp

Lines 313 to 315 in 63cd51f

	else if (fabs(m) == 1) {
	// D_stuff: d(Dp)/dt = d(Hr)/dz - d(Hz)/dr
	// B_stuff: d(Br)/dt = d(Ep)/dz - i*m*Ez/r

[oskooi]

Looks like there is no bug after all in the step-curl update at $r=0$. The issue here is rather with the quasi-PML in the r direction. This is because simply reducing the fwidth by a factor of 10 and doubling the cutoff of the pulsed source (which reduces the high-frequency components of the source which are near the Nyquist frequency of the grid that are poorly absorbed by the PML) produces flux values that are converged to 10 digits of accuracy. Doubling the PML thickness does not change the results.

Here is a comparison of two sources showing the flux values as a function of simulation time (after the source has turned off). The five data columns are: simulation time, flux +z, flux +r, flux -r, and total flux. Case 1 uses the same fwidth used to generate the results above and the default cutoff of 5.0. Case 2 uses the improved values.

1. fwidth=0.2*fcen,cutoff=5.0

flux:,  50.0, 0.0001461470, 0.0705033518, 0.0001461470, 0.0707956457
flux:, 150.0, 0.0001516641, 0.0705022176, 0.0001516641, 0.0708055459
flux:, 250.0, 0.0001462489, 0.0705008945, 0.0001462489, 0.0707933924
flux:, 450.0, 0.0001342017, 0.0704997094, 0.0001342017, 0.0707681129
flux:, 850.0, 0.0001519414, 0.0704996346, 0.0001519414, 0.0708035174

2. fwidth=0.02*fcen,cutoff=10.0

flux:, 1000.0, 0.0146318174, 7.0500586967, 0.0146318174, 7.0793223314
flux:, 1100.0, 0.0146318174, 7.0500586967, 0.0146318174, 7.0793223314
flux:, 1200.0, 0.0146318174, 7.0500586967, 0.0146318174, 7.0793223314
flux:, 1400.0, 0.0146318174, 7.0500586967, 0.0146318174, 7.0793223314
flux:, 1800.0, 0.0146318174, 7.0500586967, 0.0146318174, 7.0793223314

import meep as mp
import numpy as np


resolution = 30  # pixels/um                                                                                              

s = 5.
dpml = 4.
cell_size = mp.Vector3(s+dpml,0,s)

pml_layers = [mp.PML(dpml,direction=mp.R)]

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
    )
)

import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
sim.plot2D()
plt.savefig('cyl_sim.png',dpi=150,bbox_inches='tight')

for t in [0, 100, 200, 400, 800]:
    sim.run(until_after_sources=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'm not sure it's the quasi-PML either. The basic issue is that at the Nyquist frequency the group velocity is zero, so the fields have a hard time decaying and escaping.

I'm guessing that the problem at the Nyquist frequency is even worse in cylindrical coordinates, e.g. there might be a "waveguiding" effect that traps waves near r=0 around the Nyquist frequency, due to the r-varying discretization and numerical dispersion.

So, the solution is just to be extra careful of exciting Nyquist-frequency waves in cylindrical coordinates.

[oskooi]

For posterity, this was fixed by #2382.

[oskooi]

Whoops, I misspoke. This issue is not fixed by #2382 but rather by #2383?

[oskooi]

Looks like the bug is in the update $E$ from $D$ step of update_eh.cpp: $E=\varepsilon^{-1} D$ everywhere at $r=0$ regardless of the value of $m$. See line 218 in the code snippet below. As such, the time-stepping implementation is incorrect at $r=0$ for the $z$-PML. This could explain why we have observed the fields at $r=0$ to not decay away properly for $m = 0, \pm 1$ because the fields at $r=0$ are non-zero in those cases: (1) $E_z$ for $m=0$, and (2) $E_\phi$ and $H_r$ for $m=\pm 1$. See #2383 (comment).

The reason this bug has gone undetected until now is because it only affects the voxels at $r=0$ which therefore introduces a $\Delta r^2$ error.

meep/src/update_eh.cpp

Lines 200 to 226 in 6abb8b3

	/* Do annoying special cases for r=0 in cylindrical coords. Note
	that this only really matters for field output; the Ez and Ep
	components at r=0 don't usually affect the fields elsewhere
	because of the form of Maxwell's equations in cylindrical coords. */
	// (FIXME: handle Kerr case? Do we care about auxiliary PML fields here?)
	if (gv.dim == Dcyl && gv.origin_r() == 0.0) DOCMP FOR_FT_COMPONENTS(ft, ec) {
	if (f[ec][cmp] && (ec == Ep || ec == Ez || ec == Hr)) {
	component dc = field_type_component(ft2, ec);
	if (f[ec][cmp] == f[dc][cmp]) continue;
	const int yee_idx = gv.yee_index(ec);
	const int d_ec = component_direction(ec);
	const int sR = gv.stride(R), nZ = gv.num_direction(Z);
	realnum *E = f[ec][cmp];
	const realnum *D = f_minus_p[dc][cmp] ? f_minus_p[dc][cmp] : f[dc][cmp];
	const realnum *chi1inv = s->chi1inv[ec][d_ec];
	if (chi1inv)
	for (int iZ = 0; iZ < nZ; iZ++) {
	const int i = yee_idx + iZ - sR;
	E[i] = chi1inv[i] * D[i];
	}
	else
	for (int iZ = 0; iZ < nZ; iZ++) {
	const int i = yee_idx + iZ - sR;
	E[i] = D[i];
	}
	}
	}

[stevengj]

The comment claims "this really only matters for field output", i.e. the curl equations don't depend on those components?

[oskooi]

From Section 2.3 "Absorbing Boundary Condition" of Radio Science, Vol. 48, pp. 232-47, (2013):

"Here a generalized un-split field PML technique [Gedney, 1996; Taflove and Hagness, 2005] was employed to maximize the absorption of the outgoing waves in the PML regions in cylindrical coordinates for efficiently modeling BOR lenses within or on top of a layered medium. The detailed information and the derived formulas can be found in the reference [Taflove and Hagness, 2005]. Note that the choice of parameters used in the PML region depends on the local material properties."

There is no mention of special PML update equations at $r=0$.