[MITgcm-devel] changes in tracer vertical implicit scheme

Patrick Heimbach heimbach at MIT.EDU
Thu Mar 8 10:32:28 EST 2012 

Hi Martin,

On Mar 8, 2012, at 10:16 AM, Martin Losch wrote:

> Hi Jean-Michel,
> thanks for the info (that I should have produced myself). 
> 
> Bottom line of my quick tests: Patrick's code and it's adjoint both vectorize. Unfortunately, global_ocean.90x40x15/input_ad with Patrick's data-parameters crashes after 26 timesteps, so that I cannot do a meaningful profiling of the code, but the compiler messages are good.


I removed a few things in the "standard" data file (it was also testing Leith, etc),
which I wanted to avoid. The data file I sent doesn't set any horiz. viscosity,
that could already explain things.
Alternatively, and actually preferable is to test impl. vert. adv. in global_ocean.cs32x15
since it's an ECCO-similar setup (takes slightly longer to run, but not too bad).

> Still, I have a suggestion about moving the if-statements out of the ij-loops (not absolutely necessary) and make the loop-ranges go over the overlaps (for the case of TARGET_NEC_SX in analogy to impdiff.F, etc). But other than that, I am fine with the code. 

Can't we just have the loop-ranges go over the entire domain regardless,
such as to avoid the ugly #ifdef's? 

> What should we do now?
> 1. Keep version 1 but remove my previous additions?
> 2. Keep the KLOOPINSIDE-code as it?
> 3. Add the array-ified KLOOPINSIDE-code (KLOOPOUTSIDE) as a third option?

See my previous mail.
Perhaps adapt your above suggested changes to that version.

> What about defaults: I think that ALLOW_SOLVER_KLOOPINSIDE and potentially new flags should be set outside of solve_pentadiagonal, but others might disagree.
> 
> What about solve_tridiagonal? Same problems? same procedure?

Could be good to keep similar structure for both indeed.

p.

> Martin
> 
> On Mar 8, 2012, at 4:05 PM, Jean-Michel Campin wrote:
> 
>> Hi Martin,
>> 
>> The implicit vertical advection is tested in:
>> advect_xz/input.nlfs/
>> I did not look to the KLOOPINSIDE code (Gael might give more
>> details on this piece of code he wrote), but in this experiment,
>> switching it on makes only "small" differences (with -ieee):
>> I get the same T,S monitor, but the full binary file are not
>> identical.
>> 
>> Cheers,
>> Jean-Michel
>> 
>> On Thu, Mar 08, 2012 at 02:44:22PM +0100, Martin Losch wrote:
>>> Hi Patrick et al.,
>>> 
>>> I'll have a look at how the code vectorizes.
>>> 
>>> More generally: Basically there are two versions of the code
>>> 1. the original solver (+ the not so good version with arrays that I implemented)
>>> 2. the KLOOPINSIDE part (and your array-fied version of it)
>>> My question, does these give the same results in a forward integration? I have the impression, that they are actaully implement different algorithms, is that true? At least, it's very difficult for me to do a comparison, because the order of operations does not seem to be the case.
>>> 
>>> Martin
>>> 
>>> On Mar 8, 2012, at 1:47 PM, Patrick Heimbach wrote:
>>> 
>>>> Hi Martin,
>>>> 
>>>> sorry I haven't had time to compare the two vectorized codes.
>>>> If the one I sent you vectorizes well on the NEC we could just keep that one,
>>>> or more precicesly:
>>>> 1. keep the original ALLOW_SOLVERS_KLOOPINSIDE code (less memory footprint)
>>>> 2. alternatively use the new one for vector which was derived from 1.
>>>> That would also make the relationship between the two versions more transparent.
>>>> 
>>>> I'd be reluctant to carry three pieces of code around ( = pkg/seaice ;)
>>>> Will try to look asap at the differences, but not sure when.
>>>> 
>>>> Cheers
>>>> -Patrick
>>>> 
>>>> On Mar 8, 2012, at 6:28 AM, Martin Losch wrote:
>>>> 
>>>>> Hi Patrick,
>>>>> 
>>>>> At first glance, you modifications look very much like what I have done, but I have to check more carefully. In the meantime I ran global_ocean.90x40x15 with your "data" file and the implvertadv flag set. I used version 1.4 of solve_pentadiagonal.F to create a reference and got this:
>>>>> 
>>>>> using v1.4 as a reference:
>>>>> update to solve_pentadiagonal.F v1.6 without #define ALLOW_SOLVERS_KLOOPINSIDE in CPP_OPTIONS.h
>>>>> Y Y Y Y 14> 0< 7 FAIL  global_ocean_test
>>>>> Start time:  Thu Mar  8 12:11:42 CET 2012
>>>>> End time:    Thu Mar  8 12:15:35 CET 2012
>>>>> [This is probably the problem that you described]
>>>>> 
>>>>> solve_pentadiagonal.F v1.6 + #define ALLOW_SOLVERS_KLOOPINSIDE in CPP_OPTIONS.h
>>>>> Y Y Y Y 16>16<16 pass  global_ocean_test
>>>>> Start time:  Thu Mar  8 11:52:52 CET 2012
>>>>> End time:    Thu Mar  8 11:55:39 CET 2012
>>>>> [Problem solved]
>>>>> 
>>>>> with TARGET_SX_NEC defined in solve_pentadiagonal.F
>>>>> Y Y Y Y 14> 0< 7 FAIL  global_ocean_test
>>>>> Start time:  Thu Mar  8 12:21:13 CET 2012
>>>>> End time:    Thu Mar  8 12:24:03 CET 2012
>>>>> [My code does not work, bummer]
>>>>> 
>>>>> with Patrick's routine
>>>>> Y Y Y Y 16>16<16 pass  global_ocean_test
>>>>> Start time:  Thu Mar  8 12:02:36 CET 2012
>>>>> End time:    Thu Mar  8 12:05:24 CET 2012
>>>>> [Problem also solved]
>>>>> 
>>>>> I admit that commenting the #define ALLOW_SOLVERS_KLOOPINSIDE was maybe not the wisest  move, but I did announce that in the devel-list (although not very clearly: <http://mitgcm.org/pipermail/mitgcm-devel/2011-October/005004.html>)
>>>>> As a quick fix, we could uncomment these line again and put an #ifndef TARGET_NEC_SX around them? And then I have to have a closer look at the vectorization of Patrick's code (but why doesn't my code work? It's so similar ...)
>>>>> 
>>>>> Martin
>>>>> 
>>>>> On Mar 7, 2012, at 6:07 PM, Patrick Heimbach wrote:
>>>>> 
>>>>>> Hi Martin,
>>>>>> 
>>>>>> doesn't work with your TARGET_NEC_SX flag either.
>>>>>> I didn't have time to look at your implementation more carefully.
>>>>>> What I did do though (easier+quicker) and what works is
>>>>>> to use the k-loop code, array-ify it,
>>>>>> switch loop orders (from k-inner to k-outer).
>>>>>> In this case the adjoint continues to work fine,
>>>>>> and perhaps this one vectorizes fine (similar idea as your addition of arrays).
>>>>>> Routine is attached for you to try.
>>>>>> A good chance it might be further improved.
>>>>>> 
>>>>>> Use global_ocean.90x40x15 for testing.
>>>>>> In CPP_OPTIONS.h set
>>>>>> #define INCLUDE_IMPLVERTADV_CODE
>>>>>> 
>>>>>> Use attached "data" file.
>>>>>> 
>>>>>> p.
>>>>>> 
>>>>>> <data><solve_pentadiagonal.F>
>>>>>> 
>>>>>> On Mar 7, 2012, at 3:29 AM, Martin Losch wrote:
>>>>>> 
>>>>>>> Oh, me again?
>>>>>>> 
>>>>>>> I went through the change of Oct13, and they are (almost) all within #ifdef TARGET_NEC_SX except for 
>>>>>>> lines 197,200,205-7 (introduced a temporary value, should be unproblematic I think).
>>>>>>> lines 219 (added a CADJ loop=sequential, that is probably not without effect)
>>>>>>> and most importantly I removed the hard coding of the ALLOW_SOLVERS_KLOOPINSIDE.
>>>>>>> 
>>>>>>> Can you try to use current code with 
>>>>>>> #define ALLOW_SOLVERS_KLOOPINSIDE
>>>>>>> That will turn on Gael's (?) code for the adjoint. But you could also try to use the vectorizing code (either by removing the TARGET_NEC_SX flags, or by setting this flag). If you give me verification experiment, I'll test that, too.
>>>>>>> 
>>>>>>> I am glad, that I am so far way and out of the reach of the relentless MITgcm-jurisdiction.
>>>>>>> 
>>>>>>> Martin
>>>>>>> 
>>>>>>> On Mar 7, 2012, at 7:16 AM, Patrick Heimbach wrote:
>>>>>>> 
>>>>>>>> Hi Jean-Michel,
>>>>>>>> 
>>>>>>>> it appears that the broken impl.vert.adv. adjoint is not related to your changes below,
>>>>>>>> but happened much earlier (Oct 13, 2011), and we didn't notice since it isn't tested yet
>>>>>>>> (the culprit left the country just in time ;)
>>>>>>>> 
>>>>>>>> Reverting to solve_pentadiagonal.F revision 1.4
>>>>>>>> fixes the problems with the adjoint.
>>>>>>>> 
>>>>>>>> The changes to vectorize this one need to be revisited.
>>>>>>>> We should also start testing this in global_ocean.cs32x15
>>>>>>>> 
>>>>>>>> -Patrick
>>>>>>>> 
>>>>>>>> On Feb 1, 2012, at 5:20 PM, Jean-Michel Campin wrote:
>>>>>>>> 
>>>>>>>>> Hi Patrick,
>>>>>>>>> 
>>>>>>>>> The 2 majors changes (on Dec 01, 2011) are:
>>>>>>>>> 1) new argument "recip_hFacNew" computed in thermodynamics.F
>>>>>>>>> and passed to gad_implicit_r.F , and then passed to the selected
>>>>>>>>> gad_*_impl_r.F S/R ;
>>>>>>>>> 2) the term: tracer* del_k(wtrans)* deltaT has been moved from
>>>>>>>>> gad_implicit_r.F (in version 1.14i of gad_implicit_r.F, lines 210 - 213):
>>>>>>>>>    gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
>>>>>>>>> &      + deltaTLev(1)*recip_rA(i,j,bi,bj)
>>>>>>>>> &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
>>>>>>>>> &       *tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))*rkSign
>>>>>>>>> to gad_calc_rhs.F :
>>>>>>>>> in gad_calc_rhs.F, now we have line 164:
>>>>>>>>> IF (implicitAdvection) rAdvFac = rkSign
>>>>>>>>> but before it was:
>>>>>>>>> IF (implicitAdvection) rAdvFac = 0. _d 0
>>>>>>>>> 
>>>>>>>>> What could be tried is to have only 1 of the 2 changes, and see 
>>>>>>>>> how this goes.
>>>>>>>>> 
>>>>>>>>> Cheers,
>>>>>>>>> Jean-Michel
>>>>>>>>> 
>>>>>>>>> On Wed, Feb 01, 2012 at 12:01:16PM -0500, Patrick Heimbach wrote:
>>>>>>>>>> 
>>>>>>>>>> Jean-Michel,
>>>>>>>>>> 
>>>>>>>>>> for the record, are you aware that your changes to tracer vertical implicit scheme
>>>>>>>>>> break the adjoint (which we now use in several production calculations)?
>>>>>>>>>> 
>>>>>>>>>> Easy way to test  (I had thought it's already tested by default) is in
>>>>>>>>>> verification/global_ocean.cs32x15/input_ad/ 
>>>>>>>>>> to set
>>>>>>>>>> tempImplVertAdv=.TRUE.,
>>>>>>>>>> saltImplVertAdv=.TRUE.,
>>>>>>>>>> tempVertAdvScheme=3,
>>>>>>>>>> saltVertAdvScheme=3,
>>>>>>>>>> 
>>>>>>>>>> If you point me to likely culprit routines, I'll take a look.
>>>>>>>>>> 
>>>>>>>>>> Thanks
>>>>>>>>>> -Patrick
>>>>>>>>>> 
>>>>>>>>>> On Dec 1, 2011, at 9:43 AM, Jean-Michel Campin wrote:
>>>>>>>>>> 
>>>>>>>>>>> Hi,
>>>>>>>>>>> 
>>>>>>>>>>> Some fix/changes related to tracer vertical direction implicit method:
>>>>>>>>>>> 
>>>>>>>>>>> 1) implicit vertical advection conservation with Adams-Bashforth 
>>>>>>>>>>> and/or non-lin free-surf:
>>>>>>>>>>> the term Tr*d/dz(w) was not going through AB and neither was
>>>>>>>>>>> rescaled (NonLin-FreeSurf) and this was breaking tracer conservation.
>>>>>>>>>>> I remove Tr*d/dz(w) in gad_implicit_r.F and add it in gad_calc_rhs.F;
>>>>>>>>>>> this fix conservation in both cases.
>>>>>>>>>>> 2) implicit vertical diffusion or advection with non-lin free-surf:
>>>>>>>>>>> Was using current recip_hFacC instead of future recip_hFacC.
>>>>>>>>>>> This lead to small error when using NonLin-FreeSurf and also break 
>>>>>>>>>>> conservation with NonLin-FreeSurf in r-coordinate (but not with r*).
>>>>>>>>>>> The right recip_hFac is now computed in thermodynamics.F and passed
>>>>>>>>>>> to all vertical implicit S/R.
>>>>>>>>>>> 
>>>>>>>>>>> I also added diagnostics for gtNm1 & gsNm1, so that it's easier
>>>>>>>>>>> to check tracer conservation with NonLin-FreeSurf + AB (the quantity
>>>>>>>>>>> being conserved is Tr + (0.5+abEps)*deltaT*gtNm1).
>>>>>>>>>>> And will add a secondary test to advect_xz to test implicit vertical
>>>>>>>>>>> advection (+ NonLin-FreeSurf).
>>>>>>>>>>> 
>>>>>>>>>>> I don't know why those changes affect bottom_ctrl_5x5 (adjoint gradient)
>>>>>>>>>>> and break agreement with forward gradient.
>>>>>>>>>>> I removed all the implicit vertical code in thermodynamics.F (checked-in),
>>>>>>>>>>> and it's better, but still not as good than what it used to be.
>>>>>>>>>>> If someone has time to take a look at this experiment.
>>>>>>>>>>> 
>>>>>>>>>>> Cheers,
>>>>>>>>>>> Jean-Michel
>>>>>>>>>>> 
>>>>>>>>>>> _______________________________________________
>>>>>>>>>>> MITgcm-devel mailing list
>>>>>>>>>>> MITgcm-devel at mitgcm.org
>>>>>>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>>>>>> 
>>>>>>>>>> ---
>>>>>>>>>> Patrick Heimbach | heimbach at mit.edu | http://www.mit.edu/~heimbach
>>>>>>>>>> MIT | EAPS 54-1420 | 77 Massachusetts Ave | Cambridge MA 02139 USA
>>>>>>>>>> FON +1-617-253-5259 | FAX +1-617-253-4464 | SKYPE patrick.heimbach
>>>>>>>>>> 
>>>>>>>>>> 
>>>>>>>>>> 
>>>>>>>>>> _______________________________________________
>>>>>>>>>> MITgcm-devel mailing list
>>>>>>>>>> MITgcm-devel at mitgcm.org
>>>>>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>>>>> 
>>>>>>>>> _______________________________________________
>>>>>>>>> MITgcm-devel mailing list
>>>>>>>>> MITgcm-devel at mitgcm.org
>>>>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>>>> 
>>>>>>>> ---
>>>>>>>> Patrick Heimbach | heimbach at mit.edu | http://www.mit.edu/~heimbach
>>>>>>>> MIT | EAPS 54-1420 | 77 Massachusetts Ave | Cambridge MA 02139 USA
>>>>>>>> FON +1-617-253-5259 | FAX +1-617-253-4464 | SKYPE patrick.heimbach
>>>>>>>> 
>>>>>>>> 
>>>>>>>> 
>>>>>>>> _______________________________________________
>>>>>>>> MITgcm-devel mailing list
>>>>>>>> MITgcm-devel at mitgcm.org
>>>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>>> 
>>>>>>> 
>>>>>>> _______________________________________________
>>>>>>> MITgcm-devel mailing list
>>>>>>> MITgcm-devel at mitgcm.org
>>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>>> 
>>>>>> ---
>>>>>> Patrick Heimbach | heimbach at mit.edu | http://www.mit.edu/~heimbach
>>>>>> MIT | EAPS 54-1420 | 77 Massachusetts Ave | Cambridge MA 02139 USA
>>>>>> FON +1-617-253-5259 | FAX +1-617-253-4464 | SKYPE patrick.heimbach
>>>>>> 
>>>>>> 
>>>>>> _______________________________________________
>>>>>> MITgcm-devel mailing list
>>>>>> MITgcm-devel at mitgcm.org
>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>> 
>>>>> 
>>>>> _______________________________________________
>>>>> MITgcm-devel mailing list
>>>>> MITgcm-devel at mitgcm.org
>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>> 
>>>> ---
>>>> Patrick Heimbach | heimbach at mit.edu | http://www.mit.edu/~heimbach
>>>> MIT | EAPS 54-1420 | 77 Massachusetts Ave | Cambridge MA 02139 USA
>>>> FON +1-617-253-5259 | FAX +1-617-253-4464 | SKYPE patrick.heimbach
>>>> 
>>>> 
>>>> 
>>>> _______________________________________________
>>>> MITgcm-devel mailing list
>>>> MITgcm-devel at mitgcm.org
>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>> 
>>> 
>>> _______________________________________________
>>> MITgcm-devel mailing list
>>> MITgcm-devel at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>> 
>> _______________________________________________
>> MITgcm-devel mailing list
>> MITgcm-devel at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
> 
> 
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel

---
Patrick Heimbach | heimbach at mit.edu | http://www.mit.edu/~heimbach
MIT | EAPS 54-1420 | 77 Massachusetts Ave | Cambridge MA 02139 USA
FON +1-617-253-5259 | FAX +1-617-253-4464 | SKYPE patrick.heimbach

More information about the MITgcm-devel mailing list