[MITgcm-devel] changes in tracer vertical implicit scheme

Patrick Heimbach heimbach at MIT.EDU
Wed Mar 7 12:07:22 EST 2012 

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.

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


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

More information about the MITgcm-devel mailing list