Open mochen4 opened 10 months ago
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: https://github.com/NanoComp/meep/blob/b065eae0491a7987983779c16115200acf0085d2/src/sources.cpp#L307 looks insufficient. If idx_vol < npts
, we should really resize the two arrays to be smaller (length = idx_vol
).
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
It also appears that more pixels need to be excluded (larger offset
) for larger $|m|$.
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
See also the accurate_fields_near_cylorigin
option in the Python API, which documents this behavior.
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.
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 https://github.com/NanoComp/meep/issues/2644#issuecomment-1732589039 and find the largest Courant factor at each m value.
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.
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.