[MITgcm-support] numerical stability criteria

[Matthew Mazloff]

Hi David,

Nothing is wrong....these stability formulae are guidelines.  What you 
are seeing is that the chosen scales and the chosen constants appended 
to the formulae are not exactly accurate for every problem.  The 
stability parameters give you a first guess at dt (and let you see what 
the model is most sensitive too).  Through trial and error you can push 
these stability limits, in a sense to discover the true scaling in the 
model.

I will just explain question (1), but the same goes for 2 and 3

1)  what is a stable dt for horizontal Laplacian friction?
Well if friction diffuses momentum (in a sense accelerates it) faster 
than the model can time step it....a shock can occur.  This is a very 
lame hand wavy explanation....I know there is a rigorous mathematical 
way of showing this but I can't remember it right now (Alistair, where 
did your teaching website go!).  Anyway, this is the basic idea.  So one 
tries to estimate the order of the diffusion by scaling the horizontal 
Laplacian as 1/X/X.   And we then choose the grid spacing to be X which 
makes great sense in a numerical finite difference problem.  Then it's 
assumed that one must keep 1/dt > Ah/X/X.  Through trial and error, a 
safe level below this value is determined and the result is the formula 
in the manual. 

What you have found is that in some problems the magnitude of Laplacian 
diffusion (and other things) is less than the value chosen from basic 
finite-difference scaling.  Thus, it is safe to increase dt.

I know this is a rather incomplete explanation (and I really hope its 
accurate)....anyway, its the basic idea as to why these stability 
parameters are not set in stone, rather they are guidelines.

hope this helps,
Matt


David Wang wrote:

> Hello all,
>
> I was checking the numerical stability of my setup against the 
> criteria given in the chapter 3.10.2.1 and found these calculations in 
> the manual don't seem to be all correct.
>
> 1. for the stability parameter to the horizontal Laplacian friction Sl 
> = 4*Ah*deltaTmom/delX^2,
> for this particular 4x4 gloal ocean example run described in the 
> manual, it's calculated to be 4*5E5*(40*60)/(77,000^2) = 0.81, which 
> is above the 0.3 upper limit. In 
> verification/global_ocean.90x40x15/input/data namelist, deltaTmom is 
> 1200, which however still gives Sl as large as 0.4, not 0.16 in the 
> manual.
>
> 2. for the stability parameter to the vertical Laplacian dissipation 
> Sl = 4*Az*deltaTmom/delZ^2 = 4*1E-3*(40*60)/(50^2) = 0.004, not 0.015. 
> Also, what its upper stability limit should be? It's not given in the 
> text.
>
> 3. the numerical stability for inertial oscillations Si = 
> f^2*deltaTmom^2 = (1.43E-4)^2 * (2400)^2 = 0.12, not 0.24
>
> I don't doubt the formulae, so what else could be wrong?
>
> Thanks,
> - D.W.
>