[MITgcm-devel] bug: division by zero

Thu May 31 10:28:16 EDT 2007

Hi Jean-Michel,

initializing with 1e10 does not change the cs-verification  
experiments, nor does it change my llc-configuration results.

Should we check this in then (with the appropriate comments), or  
should I rather keep this as a local solution?

Martin

On 31 May 2007, at 15:24, Jean-Michel Campin wrote:

> 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
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel