[MITgcm-devel] changes in tracer vertical implicit scheme
Jean-Michel Campin
jmc at ocean.mit.edu
Thu Mar 8 11:03:13 EST 2012
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
More information about the MITgcm-devel
mailing list