[MITgcm-devel] changes in tracer vertical implicit scheme
Martin Losch
Martin.Losch at awi.de
Thu Mar 8 11:09:18 EST 2012
Hi there,
I just discussed with Patrick (via telefone) that I will leave this to you guys at MIT, and I will adjust the code to my beloved vector machine, once it is there.
Martin
On Mar 8, 2012, at 5:03 PM, Jean-Michel Campin wrote:
> Hi,
>
> Other comments:
> 1) the initial version was cheaper in term of number of 3-d arrays
> (6 now versus none before). Should we care about this ?
> 2) I would prefer to keep the full loop range within #ifdef TARGET_NEC_SX,
> but if we decide to always do full loop range, I would at least keep the
> other (DO j=jMin,jMax & DO i=iMin,iMax) commented out (and not completely
> removed). The reasons I see are:
> a) not so clear that the 4 input arguments iMin,iMax,jMin,jMax should be ignored.
> b) this introduces an other difference between the KLOOPINSIDE code
> (which will use iMin,iMax,jMin,jMax loop range) and the default
> full range loop.
>
> Cheers,
> Jean-Michel
>
> On Thu, Mar 08, 2012 at 10:32:28AM -0500, Patrick Heimbach wrote:
>> Hi Martin,
>>
>> On Mar 8, 2012, at 10:16 AM, Martin Losch wrote:
>>
>>> Hi Jean-Michel,
>>> thanks for the info (that I should have produced myself).
>>>
>>> Bottom line of my quick tests: Patrick's code and it's adjoint both vectorize. Unfortunately, global_ocean.90x40x15/input_ad with Patrick's data-parameters crashes after 26 timesteps, so that I cannot do a meaningful profiling of the code, but the compiler messages are good.
>>
>>
>> I removed a few things in the "standard" data file (it was also testing Leith, etc),
>> which I wanted to avoid. The data file I sent doesn't set any horiz. viscosity,
>> that could already explain things.
>> Alternatively, and actually preferable is to test impl. vert. adv. in global_ocean.cs32x15
>> since it's an ECCO-similar setup (takes slightly longer to run, but not too bad).
>>
>>> Still, I have a suggestion about moving the if-statements out of the ij-loops (not absolutely necessary) and make the loop-ranges go over the overlaps (for the case of TARGET_NEC_SX in analogy to impdiff.F, etc). But other than that, I am fine with the code.
>>
>> Can't we just have the loop-ranges go over the entire domain regardless,
>> such as to avoid the ugly #ifdef's?
>>
>>> What should we do now?
>>> 1. Keep version 1 but remove my previous additions?
>>> 2. Keep the KLOOPINSIDE-code as it?
>>> 3. Add the array-ified KLOOPINSIDE-code (KLOOPOUTSIDE) as a third option?
>>
>> See my previous mail.
>> Perhaps adapt your above suggested changes to that version.
>>
>>> What about defaults: I think that ALLOW_SOLVER_KLOOPINSIDE and potentially new flags should be set outside of solve_pentadiagonal, but others might disagree.
>>>
>>> What about solve_tridiagonal? Same problems? same procedure?
>>
>> Could be good to keep similar structure for both indeed.
>>
>> p.
>>
>>> Martin
>>>
>>> On Mar 8, 2012, at 4:05 PM, 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
>>>
>>>
>>> _______________________________________________
>>> 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
More information about the MITgcm-devel
mailing list