[MITgcm-devel] changes in tracer vertical implicit scheme

Martin Losch Martin.Losch at awi.de
Thu Mar 8 08:44:22 EST 2012 

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

More information about the MITgcm-devel mailing list