[Mitgcm-support] Re: GM c48

Wed Jul 9 15:55:46 EDT 2003

With the latest code in repository
- the following options seem to work with the adjoint
   'clipping', 'dm95', 'gkw91'
- the following options don't work
   'linear', 'gkw91', 'ldd97'

Yes, the constant exponents were bugs,
not sure about the km1 issue.

Surprising to me to need Large_SlopeSqr = 1.e48
even surprising to me that they can be as large as 1.e26
Doesn't that hint at another problem somewhere in GM?

Good idea to replace the GMREDI_OPTIONS by GM_EXCLUDE_...

Splitting into the multiplication by inverse helps to
avoid recomputations. I hadn't looked at the slope_psi part
in detail yet (cf. tag-index).

-Patrick

Jean-Michel Campin wrote:
> Hi guys,
> 
> I started to look at c48 version of GM. 
> 
> o The fact that the 4 experiments that used GM fail the testscript 
>   in both checkpoint47g_post & checkpoint48 is annoying (for us and
>   also for others that download the code, run an example and
>   cannot reproduce the standard results).
> 
> o Things are not simple. gmredi_slope_limit.F is becoming almost
>   as hard to read as pieces of code that I usually write myself.
>   Let's concentrate on the 2 S/R that are responsible for the output 
>   differences: gmredi_slope_limit.F & gmredi_calc_tensor.F
> 
> ----------------------
> gmredi_slope_limit.F 
>  a) a 2nd tapering scheme has been added over all the others.
>    This is enough to change the results:
> Y Y Y Y  8 13 10 10  9 16 12 12  8  7  9  8  7  8  9  7  7 FAIL  global_ocean.90x40x15
> Y Y Y Y  8 13 10 10  9 16 11 13  9  7  9  8  8  8  9  7  8 FAIL  global_with_exf  
>    This is not a compiler Optimization issue, this is a real change in
>    the code. Therefore, it would nice to put the 2nd slope_max (=Large_SlopeSqr)
>    in data.gmredi so that one can recover the original code & results
>    setting e.g. Large_SlopeSqr = 1.e48 (I checked that it works). 
>    And for the adjoint, Large_SlopeSqr= 1.e8.
> 
>  b) A 2 steps computation (multiplication by the inverse) has been replaced
>    by a single division:
>    c48, line 250,251 :
>            SlopeX(i,j) = -dSigmaDx(i,j)/dSigmaDrReal(i,j)
>            SlopeY(i,j) = -dSigmaDy(i,j)/dSigmaDrReal(i,j) 
>    47f_post:
>            dRdSigmaLtd(i,j) = 1./(dSigmaDrReal(i,j))
>            SlopeX(i,j)=-dSigmaDx(i,j)*dRdSigmaLtd(i,j)
>            SlopeY(i,j)=-dSigmaDy(i,j)*dRdSigmaLtd(i,j)
>    This is enough to change the results:
> Y Y Y Y  9 16 16 16 16 16 16 16 16 13 13 12 16 12 11 10 13 FAIL  front_relax 
> Y Y Y Y 12 16 16 16 16 16 16 16 16 13 13 12 13 13 13 12 13 FAIL  lab_sea 
>    This is not sensitive to the Optimization level (I get differences also
>    if I use no optimization at all).
>    We can decide to change the output of those 2 exp. But I would prefer
>    to change at the same time both the tapering part and the clipping part
>    (right now, only the tapering part has been simplified).
>    What I propose is:
>     1) return to the old form.
>     2) make a tag
>     3) replace the *dRdSigmaLtd by /dSigmaDrReal in BOTH the tapering part
>        and the clipping part. (if it work with TAF for one, it should work
>        also for the other ?). 
>     4) change the output that are affected. put some clear comments that
>        explain what causes the differences.
>     5) make a tag.
> 
>  c) GM_Small_Number=1.D-12 (but not used in gmredi_slope_limit.F)
>    Small_Number=1.D-20 (before was Small_Number=1.D-12)
>    Again, this is enough to change the results: 
> Y Y Y Y 12 16 16 13 13 16 16 16 13 12 12 12 12 11 12 12 12 FAIL  global_ocean.90x40x15
> Y Y Y Y 12 16 16 13 13 16 16 16 13 12 12 12 12 11 12 12 12 FAIL  global_with_exf 
> 
>  d) Cspd=2. _d 3 (before was  Cspd=2.) in the ldd97  tapering scheme.
>    This looks like a bug ?
>    this affects the results of lab_sea:
> Y Y Y Y  1  3  4  5  4  5  5  6  4  3  4  2  3  3  4  2  3 FAIL  lab_sea 
>  
> ----------------------
> gmredi_calc_tensor.F
>    c48, line 271:
>      &          *_maskW(i,j,k,bi,bj)*maskp1
>    c48, line 344:
>      &          *_maskS(i,j,k,bi,bj)*maskp1
>    *maskp1 was not there before and it should not be there.
>    I think that this breaks the "adiabatic property' of the Redi tensor.
>    Anyway, this is enough to change the results:
> Y Y Y Y  6 16  8  8  7 16  9 12  6  6  6  7  5  5  6  5  5 FAIL  global_ocean.90x40x15
> Y Y Y Y  6 16  7  8  7 16  9 12  6  6  6  7  5  5  6  5  5 FAIL  global_with_exf
>    
>    Patrick: could you check that the adjoint is OK with the right 
>      form, i.e., with no *maskp1
> ----------------------
> 
>   Concerning the advective form of GM (with GM_AdvSeparate T or F),
>   I can reproduce the results of c46. 
>  
> Question for Patrick:
>   In gmredi_calc_tensor.F, the TAF specific initializations are done
>   in the 1rst k loop, DO k=2,Nr ; why not in a specific block, before 
>   this K=2,Nr loop, and this time from k=1,Nr ? 
> 
> what I will do now: try to reach point 2 of note.b above
>   
> See you,
> 
> Jean-Michel

-- 
_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
Patrick Heimbach     Massachusetts Institute of Technology
FON: +1/617/253-5259                    EAPS, Room 54-1518
FAX: +1/617/253-4464               77 Massachusetts Avenue
mailto:heimbach at mit.edu                 Cambridge MA 02139
http://www.mit.edu/~heimbach/                       U.S.A.