[MITgcm-devel] changes in tracer vertical implicit scheme
Patrick Heimbach
heimbach at MIT.EDU
Thu Mar 8 07:47:02 EST 2012
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
More information about the MITgcm-devel
mailing list