[MITgcm-devel] changes in tracer vertical implicit scheme
Jean-Michel Campin
jmc at ocean.mit.edu
Thu Mar 8 10:41:22 EST 2012
Hi,
Just a short comment:
> What about defaults: I think that ALLOW_SOLVER_KLOOPINSIDE and potentially new i
> flags should be set outside of solve_pentadiagonal, but others might disagree.
the setting of ALLOW_SOLVER_KLOOPINSIDE cannot be moved to CPP_OPTIONS.h
without being renamed (because there is more than 1 solver in the code).
+ the "ALLOW_" part does not mean much (since once it's defined, the
KLOOPINSIDE code is used).
Cheers,
Jean-Michel
On Thu, Mar 08, 2012 at 04:16:26PM +0100, 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.
>
> 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.
>
> 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?
>
> 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?
>
> 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
More information about the MITgcm-devel
mailing list