[MITgcm-devel] bug: division by zero

Thu May 31 09:24:09 EDT 2007

Hi Martin,

Thanks for the feedback.
There are few good things here (among them, global_ocean.cs32x15 
test is "clean").
the useAreaViscLength flag was added mainly for CS-grid, 
but I though it was more for a corner problem (where
some distances are not so well defined), and since you
have land at the corner, it would make sence for you
to be able to run with useAreaViscLength=F.

Now, with your lat-lon-cap grid, and keeping useAreaViscLength=F,
does the results change when you initialise both dxV,dyU to a 
different values, e.g., 1.e+10 instead of 1 ?
(I can do similar test with cs-32, but corners are not
 dry, so it's different).
- if no, then it's a valid solution.
- if yes, then we need to add the missing exch (BTW, I think
  it's the only one missing) and use it to exch dxV,dyU.
 but in this case and in the short term, you still have the 
 possibility to use useAreaViscLength=T.

Cheers,
Jean-Michel

On Thu, May 31, 2007 at 10:26:49AM +0200, Martin Losch wrote:
> Hi Jean-Michel,
> 
> thanks for the suggestions/clarification.
> 
> a. My viscosity parameters are
> viscC2Leith=2.0,
> viscC2LeithD=2.0,
> viscA4grid = 0.01,
> 
> I have not tried useAreaViscLength=.true., which will  most likely  
> get around the problem (because there is no division involved)
> 
> b. I have initialized dxV = dyU = 1 in ini_curvilinear grid and guess  
> what, global_ocean.cs32x15 does not change, in fact all experiments  
> in question pass (except for fizhi, but the fail with g77 anyway). So  
> it looks like that these viscosities on the halos are not really used  
> for anything.
> 
> What could be a sensible solution?
> I guess, the best is to have b-grid exchanges, but problably the most  
> difficult one, right?
> Initializing dxV = 1, and dyU = 1 seems a little odd to me ...
> 
> Martin
> 
> On 30 May 2007, at 20:38, Jean-Michel Campin wrote:
> 
> >Hi Martin,
> >
> >Few things:
> >a) we still don't have a regular test with a compiler that
> >check for uninitialised variable (Constantinos find one,
> >but we don't have the licence anymore). And I am reluctant to
> >start changing range of loops without a serious way to
> >check that it's ok.
> >But it's not the first time we have this kind of problems with
> >variable viscosity, and I am not sure that we do not need
> >an accurate value of horizontal viscosity in (at least part of)
> >the overlap (this answers partly point 2).
> >Which viscosity parameters are you using ?
> >b) in relation with point 1, I don't want to use the wrong EXCH
> >for those 2 grid-variables, because they are overwriting
> >valid values of dxV,dyU with wrong things.
> >May be it's time to write an exch_uv for B-grid variables, if it's
> >really necessary.
> >In the mean time, can you try to initialise those quantities to,
> >e.g. one instead of zero (you may even be able to check if it  
> >change the
> >results), and see how it goes.
> >
> >Jean-Michel
> >
> >On Wed, May 30, 2007 at 06:40:44PM +0200, Martin Losch wrote:
> >>Hi there,
> >>
> >>do you have an opinion on this? What should/can I do about it?
> >>
> >>Martin
> >>
> >>On 25 May 2007, at 09:16, Martin Losch wrote:
> >>
> >>>Hi,
> >>>
> >>>this is probably a bug:
> >>>
> >>>on our new sx8 I get a division by zeros in mom_calc_visc.F for
> >>>curvilinear grids on the cubed sphere/llc and variable viscosity.
> >>>The culprits are recip_dxV and recip_dyU, because they are the only
> >>>ones that are not exchanged in ini_curvilinear_grid.F
> >>>Why are they not exchanged? I gather from the comments in
> >>>ini_curvilinear_grid, that for cube sphere exchanges this has not
> >>>been sorted out and some sort of hack has been applied which is now
> >>>catching up with me. What can we do?
> >>>1. use the wrong exchanges so that dxV/dyU and their reciprocals
> >>>have values other than zero in the overlaps and keep our fingers
> >>>crossed that no-one will notice (I assume that the viscosity values
> >>>in the overlaps are not that important (o:)
> >>>2. modify mom_calc_visc, so that it does not compute anything in
> >>>the overlaps (will that work?)
> >>>
> >>>solution 1 changes the results of testreport for
> >>>global_ocean.cs32x15, if I uncomment these two lines in
> >>>ini_curvilinear_grid.F:
> >>>c     CALL EXCH_Z_3D_RS( dxV, 1, myThid )
> >>>c     CALL EXCH_Z_3D_RS( dyU, 1, myThid )
> >>>I could do this:
> >>>#ifdef TARGET_SX_NEC
> >>>    CALL EXCH_Z_3D_RS( dxV, 1, myThid )
> >>>    CALL EXCH_Z_3D_RS( dyU, 1, myThid )
> >>>#endif
> >>>which will fix the problem for me but as soon as you run into a
> >>>compiler that is as picky as the sxf90 you'll have the problem  
> >>>again.
> >>>
> >>>(BTW. sxf90 has options to initial all stack/heap with nan. If I do
> >>>that the model stops in mom_vi_hdissip.F, havent' figured out yet,
> >>>why.)
> >>>
> >>>What do you think?
> >>>
> >>>Martin
> >>>
> >>>PS. I hate going to new machines, all these problems that you
> >>>discover ...
> >>>
> >>>_______________________________________________
> >>>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