improve stability of lorentzian susceptibility

[oskooi]

This is a simple fix suggested by Felix Schwarz (@fesc3555) a while ago on the mailing list to improve the stability of the lorentzian susceptibility.

Replaces and closes #238.

For reference, the original posting is reposted here along with @stevengj's reply.

original post: 3/15/2018

Dear Meep community,

I sent a problem a couple of days ago and tracked it down a bit. the
problem was, that if you define a non-diagonal sigma in your
lorentzian_susceptibility, the calculation becomes an unstable and
fields and polarizations grow exponentially.

The problem, I think, is a bug in the library and has something to do
with the averaging on the yee-grid. The polarization explodes at the
boundary of the regions where sigma!=0.  To see that in detail, I added
I diagnostic output line inside lorentzian_susceptibility::update_P in
susceptibility.cpp, around line 223. The Line (with a bit of context) is:
LOOP_OVER_VOL_OWNED(gv, c, i) {
          realnum pcur = p[i];
          p[i] = gamma1inv * (pcur * (2 - omega0dtsqr_denom)
                              - gamma1 * pp[i]
                              + omega0dtsqr * (s[i] * w[i]
                                               + OFFDIAG(s1,w1,is1,is)));
          pp[i] = pcur;
          if(p[i] > 1e30){ std::cout << i << " " << s[i] << " " <<
s[i+is1] << " " << s[i-is] << " " << s[i-is+is1] << std::endl; exit(1);}
}

it gets triggered, and the output is all zeros. However, s1 (the
non-diagonal part) is non-zero, which is, I guess, because of the
yee-grid (my material_function sets diagonal and non-diagonal elements
at the same time).

I do not know, why p explodes, though. In principle, even if a p is set
somewhere outside the metal region because of some yee-grid shift, it
should decay as fast as anywhere else, or not?

The whole problem, btw, disappears if you just put
if(s[i] != 0){...
inside the LOOP_OVER_VOL_OWNED, but that seems like a shitty solution
(just ignoring all the nice yee-grid stuff)

If you want to reproduce the Error, a minimal(ish) example is on my github:
https://github.com/fesc3555/meep_nondiagonal_sigma

Any input on this would be apprechiated.

Thanks!

reply: 3/15/2018

> I do not know, why p explodes, though. In principle, even if a p is set
> somewhere outside the metal region because of some yee-grid shift, it
> should decay as fast as anywhere else, or not?
> 
> The whole problem, btw, disappears if you just put
> if(s[i] != 0){...
> inside the LOOP_OVER_VOL_OWNED, but that seems like a shitty solution
> (just ignoring all the nice yee-grid stuff)

I don't really understand why it's unstable at first glance, but I agree that 
the source of the problem here is the Yee grid: you are apparently at a Yee 
grid point that is outside your metal, but an adjacent off-diagonal component 
is inside the metal.

The "if (s[i] != 0)" workaround seems fine; it won't degrade the accuracy since 
we don't have subpixel smoothing for Lorentzian susceptibilities anyway.  If 
you want to submit a pull request with this fix that would be great.