numerical stability in cylindrical coordinate

[mochen4]

RuntimeError: meep: simulation fields are NaN or Inf occured at the beginning of field updates when the source includes $r=0$ point at $|m| > 1$. I used a custom amplitude function such that the amplitude was 0 at $r=0$, and this check did pass.

This issue doesn't occur at $m = -1, 0, 1$. And since the amplitude should be 0 at $r=0$ for $|m| > 1$, we could just exclude the point when we define sources for now.

[stevengj]

One thing to try would be to add a line after this line:

if (amps_array[idx_vol] == 0) continue;

Then you have to be careful because this check:

meep/src/sources.cpp

Line 307 in b065eae

if (idx_vol > npts)

looks insufficient. If idx_vol < npts, we should really resize the two arrays to be smaller (length = idx_vol).

[mochen4]

Here is a simple script that shows the issue:

import numpy as np
import meep as mp

resolution = 25
frq, dfrq = 1, 0.2
dpml, dair = 1, 1
sr, sz = dair + dpml, dpml + dair + dpml
pml_layers = [mp.PML(thickness=dpml)]
cell_size = mp.Vector3(sr, 0, sz)

Jr = lambda v3: np.sin(v3.x + 0.5*dair)
offset = 0

sources = [mp.Source(mp.GaussianSource(frq, fwidth=dfrq), component=mp.Er,
        center=mp.Vector3(0.5 * dair, 0, -0.5 * sz + dpml + 0.25*dair), size=mp.Vector3(dair-offset),
        amp_func = Jr)]
sim = mp.Simulation(cell_size=cell_size, boundary_layers=pml_layers,
    resolution=resolution, sources=sources,dimensions=mp.CYLINDRICAL, m=2,)
sim.run(until_after_sources=10)

The issue goes aways if there is some offset between the source and r=0, e.g. offset = 0.15

[mochen4]

It also appears that more pixels need to be excluded (larger offset) for larger $|m|$.

[stevengj]

You could also just try a smaller Courant number following this comment: https://github.com/NanoComp/meep/blob/7e747d2c4de03329506dd2d31bfcc53e12460029/src/step_db.cpp#L436C30-L441

[stevengj]

See also the accurate_fields_near_cylorigin option in the Python API, which documents this behavior.

[mochen4]

Numerical instability does seem to be the issue, and smaller Courant factor indeed solves the problem. Maybe we need to change the default behavior for $|m|>1$: either zero more pixels near $r=0$ or use smaller Courant factors.

[mochen4]

Here is the paper I found https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=7696718, which suggested a suspiciously large small factor.

I will do some more numerical experiments based on #2644 (comment) and find the largest Courant factor at each m value.

[mochen4]

Here is what I found in Taflove's Computational Electrodynamic. He didn't give any theoretical proof, but his results are consistent with what I found with my test script.