[MITgcm-devel] (not so) funny things happen in seaice_lsr and pickups

[Jinlun Zhang]

Martin Losch wrote:
> Hi Jinlun,
>
> this is probably a question/remark for you.
>
> I have to admit, that I never looked at the details of the LSOR 
> solution technique, but as far as I can see, it probably works best, 
> wenn the "tridiagonal" system is diagonally dominant (so abs(B)>abs(A) 
> and abs(B)>abs(C), right?
Hi Martin, you are right that the matrix should be diagonally dominant 
all the time. If it is not, the coding might have a bug.
>
> I have now encoutered situations where one of the above conditions are 
> not met, and this is exactly when the solver does not converge. I get 
> these cases in physically plausible situations, e.g. at an ice edge 
> were P/zeta/eta is large in one grid point and nearly zero in the 
> neighboring gridpoint, or where P/zeta/eta are nearly zero in an 
> isolated gridpoint, which can easily happen when a closed ice cover 
> starts opening up due to wind stress, etc.
Whatever values P/zeta/eta have, it has to be diagonally dominant 
everywhere in open water or in thick ice areas, at any time.
>
> There is a part in the solver code, that has a comment (probably 
> originally yours): "SAFEGUARD AGAINST BAD FORCING ETC." which stops 
> the iteration if the innovation u(1)-u(3) increases with iteration 
> (instead of decreasing), indicating divergence. Please confirm or 
> correct my speculation, that this test anticipates that the solver 
> does not and will not converge in a few situations, so we end up with 
> potentially bad ice velocities in a few grid points, but since the ice 
> velocities are clipped (to 40cm/s), in the next time step, the forcing 
> situation may be different and the model recovers. A different 
> probably more drastic fix would be this: check if 
> abs(B)>abs(A),abs(C), and if not, reset everything to B=1,A=C=0,and 
> rhs=U, so that in this point the old ice velocity is used.
This is only for bad forcing when I used buoy derived geostrophic winds 
that had some bullseyes. If the wind forcing is right, things should be 
OK. The ice velocities capped at 40cm/s should be taken out from the 
code. I put the lines there just for a test and forgot to take them out. 
We should not cap the velocity like that. The cap may even cause blowup. 
I believe I said that before and would like to say again: let's take 
them out.
>
> In all my runs, where this problem occurs, the forcing changes very 
> slowly (once per day) compared to the model time step of 1200sec, so 
> that "BAD FORCING" can persist for quite a long time (order 72 time 
> steps). Maybe that's an issue, which could be solved by calling the 
> dynamics solver less frequently or increasing the forcing frequency.
My experience is that less frequent calling of dyn_solver or forcing 
input would not avoid blowup. The key is to make sure it is diagonally 
dominant. It is most likely due to coding bug. Recall that when I made 
that bug in B-grid, we would have some problems. After it was fixed, it 
does not blow up.
>
> Why the b-grid code seems to be more stable beats me, maybe the 
> coriolis terms adds a bit of fun here (as it always does on a C-grid), 
> but given that the blow-ups are so rare and seemingly co-incidental 
> (even a small deviation from the trajectory make them go away).
I think the solver in B-grid is now bug free after the bug fix. But I 
don't have experience with C-grid Coriolis.
>
> Any comments? When you first came up with this solver, you must have 
> had some problems with stability, otherwise this catch "SAFEGUARD 
> AGAINST ..." wouldn't be there, right? What is your experience with that?
>
> Martin
>
> PS. I don't think that the metric terms are really the issue, I have 
> now fixed them so that in my tests everything looks perfect ...
Is the code stable now?

Good luck, Jinlun