[MITgcm-devel] lab_sea sensitivity between platforms/compilers

Jean-Michel Campin jmc at ocean.mit.edu
Fri Jul 31 12:43:52 EDT 2009 

i Martin,

Few remarks regarding this lab_sea "sensitivity" to compilers:

1) it does not seems that 1 compiler is doing something funny:
pgi results are also as different from g77 (with only 5 digits right)
than they are from ifort (only 5 digits right).

2) if I start from a pickup (let say, after 100 iter),
lab_sea (evp) diverge rapidly with ifort compared to g77 
in the same way as it does when I start from 0.
But lab_sea.lsr with ifort stay on track for few iterations (~3 or 4) 
and diverge after from g77.

3) ifort versus g77 comparaison with standard (evp) lab_sea test:
at the end of the 1rst iteration, UICE & VICE (output 64 bits)
are still identical but HEFF is already different.
This seems to indicate that seaice thermodynamics or/and advection
have something to do with this sensitivity problem.

To go further, I modify seaice_dynsolver.F by just loading 
from (64 bits) files UICE,VICE (+ call EXCH) just before SEAICE_OCEAN_STRESS 
call and loading from files FU,FV (+ EXCH) just after SEAICE_OCEAN_STRESS call

I took the files from the g77 run and checked that it did not change
g77 run. But the ifort run that I did with this modified code 
still diverge from the g77 at about the same point (may be I gained
1 or 2 more digits at the 2nd iter ? but the 3rd iter is as different
as it was with the original code).

My question: do you think reading UICE,VICE,FU,FV from file is enough 
to by-pass the seaice dynamics ?
If yes, it seems that the thermodynamics/advection is the 1rst thing to
check regarding this sensitivity problem.

Cheers,
Jean-Michel

On Mon, Jul 27, 2009 at 08:15:18PM -0400, Jean-Michel Campin wrote:
> Hi Martin,
> 
> I don't have a good method to check this.
> I was thinking of dumping intermediate state to disk and
> restarting from there. painful.
> 
> Cheers,
> Jean-Michel
> 
> On Sun, Jul 19, 2009 at 10:22:00AM +0200, Martin Losch wrote:
> > Hi Jean-Michel,
> >
> > I might be chasing a Friday afternoon blackout, but I ran lab_sea  
> > without the dynamics solver, and got very similar numbers on for my  
> > different platforms.
> > Turning off the thermodynamics did not have the same effect. But it may 
> > very well be a coincidence for my particular combination of  
> > platforms/compilers.
> >
> > How do you normally test these things?
> >
> > Martin
> >
> > On Jul 17, 2009, at 11:48 PM, Jean-Michel Campin wrote:
> >
> >> Hi Martin,
> >>
> >> Are you sure that this lab_sea "sensitivity" to compiler
> >> is not comming from the thermodynamics part ?
> >> If I look at the results of ifort/ifc (or pgi) test, compared
> >> with g77/gfortran, I still get a good results for
> >> global_ocean.cs32x15.icedyn (10 or 11 digits for cg2d)
> >> but all the lab_sea are failing with only 5 digits.
> >>
> >> I don't feel enough courage now to look in details in seaice_growth.F
> >> but do we still apply instantaneous melting when SST > t_freeze ?
> >> and do we still start from 1.m ice thickness everywhere
> >> in all those lab_sea test ?
> >>
> >> Cheers,
> >> Jean-Michel
> >>
> >> On Fri, Jul 17, 2009 at 10:00:16PM +0200, Martin Losch wrote:
> >>> Hi David (and others),
> >>>
> >>> your are right, SQRT(deltaC)/(SQRT(deltaC)+SEAICE_EPS) is order one  
> >>> or
> >>> zero (if deltaC=0), the correct expression would have been
> >>>>> deltaCreg = SQRT(deltaC)+SEAICE_EPS
> >>> to avoid divisions by zero.
> >>> And consequently, when I do that, then the sensitivity returns.
> >>> I found that when I fix zetaC to a small value (say 1.), the  
> >>> sensitivity
> >>> goes away, but since zetaC is a factor that scales everything, it  
> >>> just
> >>> means that the entire stress term is scaled down by orders of  
> >>> magnitude
> >>> (normally zetaC is order 1e10 and more), and the dynamics solver  
> >>> gives
> >>> nothing much. All that means, is that the dynamics solver is very
> >>> sensitive; and I still cannot figure out why (the same is true for  
> >>> LSR of
> >>> course), rats.
> >>>
> >>> Martin
> >>>
> >>> On Jul 17, 2009, at 4:32 PM, David Ferreira wrote:
> >>>
> >>>> Hi Martin,
> >>>>> In the EVP code there is a line
> >>>>>           deltaCreg = SQRT(MAX(deltaC,SEAICE_EPS_SQ))
> >>>>> If you replace this line with
> >>>>>           deltaCreg = SQRT(deltaC)/(SQRT(deltaC)+SEAICE_EPS)
> >>>> It doesn't look like these 2 are equivalent.
> >>>> Also, what about something like :
> >>>>
> >>>> IF deltaC >= SEAICE_EPS_SQ
> >>>>  deltaCreg =SQRT(deltaC)
> >>>> else
> >>>>  deltaCreg = SEAICE_EPS
> >>>> end
> >>>>
> >>>> it adds a "IF" but remove a "MAX" ... Not sure which one costs more
> >>>> david
> >>>>
> >>>>
> >>>>
> >>>>> the number of matching digits increase from 5 to 12 (for cg2d) for
> >>>>> the lab_sea/input experiment, between a
> >>>>> Linux csysl18 2.6.22.19-0.2-bigsmp #1 SMP 2008-12-18 10:17:03 +0100
> >>>>> i686 i686 i386
> >>>>> gfortran (gcc version 4.2.1)
> >>>>> and a
> >>>>> Darwin csysm15.dmawi.de 9.7.0 Darwin Kernel Version 9.7.0: Tue Mar
> >>>>> 31 22:52:17 PDT 2009; gfortran (gcc version 4.3.1)
> >>>>>
> >>>>> Both expressions give about the same results for large DeltaC, but
> >>>>> still I would like to find out what would be the correct  
> >>>>> replacement
> >>>>> of the SQRT(MAX(deltaC,SEAICE_EPS_SQ)). Here's what did *NOT* work:
> >>>>>           deltaCreg = max(sqrt(deltaC),seaice_eps)
> >>>>>           deltaCreg = sqrt(dmax1(deltaC,seaice_eps_sq))
> >>>>>           deltaCreg = SQRT(MAX(DBLE(deltaC),DBLE(SEAICE_EPS_SQ)))
> >>>>>           deltaCreg = SQRT(MAX(SEAICE_EPS_SQ,deltaC))
> >>>>>
> >>>>> Any idea, what might be going on? If I cannot solve this problem I
> >>>>> would suggest replacing this expression with the more stable one,
> >>>>> also in seaice_calc_viscosities (although it will change all
> >>>>> results).
> >>>>>
> >>>>> Martin
> >>>>> _______________________________________________
> >>>>> 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
> >
> > _______________________________________________
> > 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