[MITgcm-devel] changes in tracer vertical implicit scheme

Patrick Heimbach heimbach at MIT.EDU
Thu Mar 8 10:26:40 EST 2012 

Hi there,

if it is agreed that Gael's implementation is acceptable,
and if it's confirmed that the array-fied version vectorizes well,
we could use attached cleaned version, which
is exactly the same algorithm, but once with k-loop inside, once outside.
(and only one CPP flag left).

-Patrick

-------------- next part --------------
A non-text attachment was scrubbed...
Name: solve_pentadiagonal.F
Type: application/octet-stream
Size: 8437 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-devel/attachments/20120308/36517186/attachment-0001.obj>
-------------- next part --------------


On Mar 8, 2012, at 10:05 AM, 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

---
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