[MITgcm-devel] changes in tracer vertical implicit scheme

Martin Losch Martin.Losch at awi.de
Thu Mar 8 06:28:52 EST 2012 

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

More information about the MITgcm-devel mailing list