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

[oskooi]

Closes #2182.

It turns out there is a bug in the $z$-PML at $r=0$ for $m=\pm 1$. This is due to the current setup of step-db which involves manually zeroing the $D_z$ fields at $r=0$:

meep/src/step_db.cpp

Lines 342 to 347 in d370f1a

	if (ft == D_stuff) {
	ZERO_Z(f[Dz][cmp]);
	if (f_cond[Dz][cmp]) ZERO_Z(f_cond[Dz][cmp]);
	if (f_u[Dz][cmp]) ZERO_Z(f_u[Dz][cmp]);
	}
	}

The $B_r$ fields had been omitted. This meant that the $z$-PML at $r=0$ was incorrect and thus the incident fields were bouncing around indefinitely as reported in #2182.

This PR also includes a new unit test based on the results in #2182. The test verifies that the flux radiated by an $E_r$ point source at $r=z=0$ in vacuum (schematic below) converges with the simulation runtime. This test is currently failing on master.

[stevengj]

Doesn't the Yee grid for Bz put the Bz component at r=Δr/2? Since it is at r≠0, why should it be set to zero?

[oskooi]

Doesn't the Yee grid for Bz put the Bz component at r=Δr/2? Since it is at r≠0, why should it be set to zero?

That's true but we are forcing the $B_r$ fields to be zero at $r=0$, not $B_z$.

[stevengj]

No, still doesn't make sense to me — for m = ±1, the radial component can be nonzero at r=0.

By setting Br to zero at r=0, you are introducing a first-order error — it should still converge with resolution, and may coincidentally be circumventing some other problem with the z PML, but it doesn't seem like the right solution to the latter problem.

[stevengj]

If the flux in the R direction is not going to zero with increasing radial cell size, I'd like to see a picture of the fields for a large-radius simulation — how are the waves even reaching the R boundary? The fields of a point source should decay with 1/r (so that the flux goes like 1/r^2)

[oskooi]

I'd like to see a picture of the fields for a large-radius simulation

Looks like the $E_r$ fields due to an $E_r$ source at $(r,z) = (0,0)$ are decaying to zero in the radial direction for an empty cell in cylindrical coordinates.

[stevengj]

You could also look at Ez as a function of r for z ≠ 0 (since it =0 at z=0 by symmetry), since that is what creates the Poynting flux. It would also be good to simply plot the Poynting flux (through a plane) vs r.

[oskooi]

You could also look at Ez as a function of r for z ≠ 0 (since it =0 at z=0 by symmetry), since that is what creates the Poynting flux. It would also be good to simply plot the Poynting flux (through a plane) vs r.

For an empty cell and a CW source, the steady-state $E_z$ and $P_r$ (flux in $r$ direction) at $z \neq 0$ both decay with $r$ as expected.

$P_r$ decays quadratically with $r$ as is evident from the second log-log plot below.

For reference, the scripts used to generate these results are provided in two gists: 1 (flux) and 2 (fields).

Summary

When using a pulsed source, the only way to ensure that $P_z$ (flux in the $z$ direction) at $r = 0$ given either (1) a source at $r = 0$ with $m = \pm 1$, or (2) a source at $r \neq 0$ with $|m| > 0 $, is converged with respect to the runtime is, as previously noted, to use a time profile with two properties: (1) narrow bandwidth (fwidth $=\Delta f < ~0.1 f_{cen}$ ) and (2) long cutoff (cutoff $> ~10$ ). This requirement is necessary for the two test cases of the empty cell as well as the LED structure. Note that the poor $z$-PML performance is independent of the source location and value of $m$: sources at $r = 0$ and $r \neq 0$ both require a long-duration pulse. (For a source at $r > 0$, the $m=0$ case does not have this restriction.)

Unfortunately, this imposes a major computational cost because computing the flux for an $r \neq 0$ point source which requires multiple simulations for $m = 0, \pm 1, \pm 2, ...$ each $m$ of $|m| > 0$ with a long runtime means that the total runtime is necessarily long. Even though the different $m$ simulations can be run in parallel, the total cost of the 2D computation (number of CPU cycles, memory) can still be significant.

Unless we figure out how to improve the $z$-PML performance at $r = 0$, it seems that the computational savings provided by setting up the problem in cylindrical coordinates rather than 3D may not be sufficiently large to prefer cylindrical coordinates over 3D.

[oskooi]

A useful reference from Prof. Douglas H. Werner's group at Penn State EE which provides the discretized equations for FDTD in cylindrical coordinates is Radio Science, Vol. 48, pp. 232-47, (2013).

The Yee grid shown in Figure 1 of this paper is identical to our Yee grid.

Section 2.2 "Singularity Issues of Electric and Magnetic Fields on the Axis of Symmetry" states:

"Due to the singularity of the field components ( $E_\phi$, $E_z$, $H_r$ ) located on the axis of symmetry, special boundary conditions are required for the update equations in the BOR (body of revolution) FDTD method. For instance, the $E_z$ component is zero for any mode with $m > 0$, while $E_\phi$ and $H_r$ are zero for any mode with $m \neq 1$ ."

This seems to be consistent with our implementation in step_db.cpp. This also means that the change proposed in this PR is unnecessary.

Additionally, Appendix B of the paper provides the update equations for the special case of $r=0$:

It would be good to compare these equations with our implementation in step_db.cpp to make sure we do not have a bug.

[oskooi]

edit (5/26/2023): the $m = 0$ case was actually never tested with the changes in this PR. See comments in #2538.

Looks like all that is required to fix this bug is to remove the special case of $r=0$ in the update $E$ from $D$ step.

With this one change, the $z$-PML correctly attenuates fields at $r=0$ for simulations with $m=0, \pm 1$. There is no longer any need to use a narrow bandwidth pulse with long cutoff.

[stevengj]

Can you check that the field output for Er and Ep looks okay at r=0 after this PR? The comment on the code says this was originally for visualization, but maybe it was superseded by some later change.

[oskooi]

Can you check that the field output for Er and Ep looks okay at r=0 after this PR?

The ($E_r$, $E_\phi$, $E_z$) fields at $r=0$ seem to be fine with this PR.

Here are plots of the steady-state fields at $r=0$ from an $E_r$ CW point source at $(r,z)=(0,0)$ in vacuum for a $m=1$ simulation (script). The only difference in the results between this PR and master is the $E_\phi$ fields at $r=0$ (middle inset of the two figures). The fields from this PR seem to be correct whereas those from the master do not because this PR uses the correct $z$-PML at $r=0$. The $E_r$ and $E_z$ fields for this PR and master are the same.

Based on these results, I think we can go ahead and merge this PR.

[oskooi]

An addendum in case the results in this PR are referenced in the future. It is important to note that $E_z$ at $r=0$ for a $m=1$ simulation should be zero. This is one of the special cases of the field-update equations shown above. However, the $E_z$ plots (right inset) of the previous comment show that $E_z$ is nonzero at $r=0$. This is not a bug. It is an artifact of get_array which always interpolates the fields to the center of the voxel in 2d. The fields at $r=0$ are therefore interpolated to $r=0.5\Delta r$.

As a check, I verified that the raw values of the $E_z$ fields array exactly at $r=0$ within get_array_slice_chunk_loop of array_slice.cpp is indeed zero. I did a similar check for (1) the $E_\phi$ fields at $r=0$ for a $m=0$ simulation and (2) the $E_\phi$ and $E_z$ fields at $r=0$ for a $m=2$ simulation.