Lorentzian check reports instability even when stable

[FilipDominec]

I would like to open the broad question of FDTD numerical stability. By trial-and-error I isolated two independent criteria of how a definition of a material can make the simulation unstable (the third mechanism of unstable PML for grazing incidence wave has been resolved in MEEP 1.2). In the treatise that follows, I tried to describe these two criteria.

In last 3 years I tested different materials thoroughly and these rules seem to hold so far (though I have no theory why).

What does not always hold, however, is the corresponding function `lorentzian_unstable' (in src/susceptibility.cpp). MEEP checks for the position of the oscillator pole in the frequency complex plane, based on some mathematics that I admit not to understand yet, and sometimes reports false positives. In the cases near the edge of stability it reports instability and aborts.
I experimentally verified that by disabling the function, recompiling and running the same simulation again, no instability developed even in very long time.

Of course, there are also false negatives, due to the second criterion: The simulation goes unstable whenever the permittivity at high frequencies is too low.

I suggest that first of all, the `lorentzian_unstable' should be changed to a warning only, it should be called only once after the initialisation (not in n-times every step), and no matter what the outcome, the simulation should continue. In the future, we should thoroughly test out the hypothesis attached below. This starts to be a bit experimental research instead of programming, but with so many people computing some plasmonics etc, it is important to make it clear.

Looking forward to your comments!
Filip

[stevengj]

@FilipDominec, wouldn't the normal way to proceed be to perform a Von Neumann stability analysis for a homogeneous material? Trial-and-error does not inspire confidence, although it can be a useful first step.

I haven't thought about this for a while, but the lorentzian_unstable is probably overly conservative because it considers only the P equation by itself, not coupled to Maxwell. It is supposed to be just a standard z-transform analysis of a linear-ODE discretization, although of course it's possible there is a bug.

Because the check is performed outside the LOOP_OVER_VOL_OWNED, there is no significant performance penalty to performing it on any iteration, and it is easier than adding a flag to the structure to record whether the check has been performed.

[FilipDominec]

1. My original post was motivated by practical needs to compute realistic materials, rather than to investigate mathematics behind FDTD (I would be happy to do so, but I can not make too many detours in my research). I was surprised that my results, based on trial-and-error and repeatedly verified in practice, seem to contradict the theory built in MEEP, and I wanted to share this experience.

2. The scheme interface for MEEP implements a very convenient feature, that is, the definition of a material, which then can be placed into some subset of the volume. With this, both stability rules noted above can be checked. I guess this would eliminate 99 % of cases when one makes the simulation unstable. (I had to re-implement such a material definition in Python, as I wanted to seamlessly integrate the simulation with the power of scipy/matplotlib.)

3. I can see that the performance impact of `lorentzian_unstable' is negligible, but if a warning is to be issued, it should not be printed out several times each time step.

[stevengj]

I'm just reluctant to implement a stability criterion in Meep which is purely empirical, because I don't know how trustworthy it is. I'd prefer to be overly conservative than to be overly permissive.

[FilipDominec]

@stevengj I can understand your point. It is indeed beyond my capabilities (and allocated time) to mathematically prove the two rules posted above, and it may be they have no mathematical background at all or later turn out to be just an approximation.

But the currently implemented criterion is known to systematically give wrong results, and it does so in a restrictive way (which eventually made me to wipe it out of the sources and recompile meep to proceed).

I believe it can not be of any harm if this stability check for the 1st rule, whatever its implemenation be, just displays a warning and lets one run the simulation at their own risk.

[FilipDominec]

In fact, both rules I realized by blind experimentation are more or less contained in the MEEP documentation.

The first rule, limiting the frequency of an oscillator, is mentioned in 2nd paragraph here:
http://ab-initio.mit.edu/wiki/index.php/Materials_in_Meep#Numerical_stability, where the oscillator frequency is simply required to be below 1/(π Δ_t_).

The second rule is hidden e.g. in the definition of the Courant factor http://ab-initio.mit.edu/wiki/index.php/Meep_Reference, if one rewrites minimum refractive index to arguably more exact formulation minimum real part of refractive index, measured at the frequency 1/(π Δ_t_).
This is an important difference, as I ran into stability troubles mostly when I defined dispersive materials.

I think this is all in all a beautiful result; everything suggests the necessary conditions to avoid numerical instabilities in FDTD are already contained in the docs.

[FilipDominec]

@stevengj Two years later, my experience only confirms what I suggested above. The built-in stability check in MEEP should be replaced for the benefit of MEEP users - no matter that it might have been mathematically derived in any way.

It is often falsely positive (simulation is aborted when it would be perfectly stable otherwise), and falsely negative (one can design materials that pass the test, but make the simulation unstable.

[stevengj]

For now, the stability check has been disabled. If people notice a simulation blowing up, it is up to them to reduce the Courant factor.