[MITgcm-devel] mom_calc_visc.F
David Ferreira
dfer at mit.edu
Fri May 22 14:59:56 EDT 2009
Hong,
I'm not sure your problem is related to these bugs, but
if you want, I can send you the correct subroutine,
so you can try it before I check it in the code.
david
Hong Zhang wrote:
> Hi David,
> We use useAreaViscLength=.FALSE. in CS510 adjoint.
> The forward run is OK but we did have blowup at corner in backward
> integration.
> Have no idea it's related to this bug.
> But we capped the adjoint fields in the neighborhood of corners
> to get a stable adjoint run.
>
> cheers,
> hong
>
> David Ferreira wrote:
>> Again I reply to myself. But this time, I have a conclusion.
>> There are 2 bugs in mom_calc_visc.F:
>>
>> 1) with useAreaViscLength=.TRUE., the biharmonic viscosities
>> computed at the corner and center points differ by a factor 4 (the
>> corner ones
>> are wrong ones)
>>
>> 2) with useAreaViscLength=.FALSE., the computation of L4rdt is wrong
>> at both
>> positions (a 1 dropped in the formula and a 16 turned into a 6).
>>
>> I'm going to check in fixes for those 2 bugs. Both will affect the
>> runs with
>> variables viscosities, 1) probably more than 2). If you have a
>> complaint, now
>> is good time (cs510 people ?).
>> Cheers,
>> david
>>
>>
>> David Ferreira wrote:
>>> It probably means that the coefficient for the computation
>>> at the vorticity points should be changed similarly (from 1/8 to 1/32).
>>> We are left with the "1/32 or 1/20" problem depending on
>>> useAreaViscLength
>>> david
>>>
>>>
>>>
>>> chris hill wrote:
>>>> Darn - does that mean cnh has to sort it out! I'll take a look.
>>>>
>>>> As I said to David (when I did respond to his original question) -
>>>> there is no 100% _correct_ answer here. The optimal scaling of
>>>> effective dissipation/diffusion with grid spacing does not need to
>>>> follow the spacing exactly .e.g on a lat-lon grid typically the
>>>> meridional component of biharmonic mixing is scaled with a l^~2.5
>>>> not l^4 times coefficient factor. Nevertheless we should have a
>>>> consistent transfer function and be able to say what it is.
>>>>
>>>> cnh
>>>>
>>>> Jean-Michel Campin wrote:
>>>>> Hi David,
>>>>> In the CVS log-message, I found, from cnh (revision 1.25):
>>>>>> Change from 1/8 to 1/32 scaling from Baylor.
>>>>>> Need to check if there are any verif expts affected.
>>>>> May be this could help ?
>>>>>
>>>>> Jean-Michel
>>>>>
>>>>> On Mon, Mar 09, 2009 at 10:18:34PM -0400, David Ferreira wrote:
>>>>>> Since nobody answers my question, I'm doing my own
>>>>>> follow-up: Even more bugs in mom_calc_visc.F ?
>>>>>>
>>>>>> My first below concerned the computation of the viscosities
>>>>>> at the div points. There is an equivalent for at vorticity points:
>>>>>>
>>>>>> CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
>>>>>> C These are (powers of) length scales
>>>>>> IF (useAreaViscLength) THEN
>>>>>> L2=rAz(i,j,bi,bj)
>>>>>> L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
>>>>>> ELSE
>>>>>> L2=2. _d
>>>>>> 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
>>>>>> L4rdt=recip_dt/
>>>>>> & ( 6. _d
>>>>>> 0*(recip_DXV(I,J,bi,bj)**4+recip_DYU(I,J,bi,bj)**4)
>>>>>> & +8. _d
>>>>>> 0*((recip_DXV(I,J,bi,bj)*recip_DYU(I,J,bi,bj))**2))
>>>>>> ENDIF
>>>>>>
>>>>>> As for the divergence, the 2 methods for a carsetian square grid
>>>>>> converge for L2 but not L4rdt. Even weirder, the coefficients
>>>>>> are now 1/8 and 1/20 (compared to 1/32 and 1/20 at the div points).
>>>>>>
>>>>>> In the best case scenario, that's only 2 correct coeff out of 4,
>>>>>> no ?
>>>>>> Comments ?
>>>>>> david
>>>>>>
>>>>>>
>>>>>>
>>>>>> David Ferreira wrote:
>>>>>>> Hi,
>>>>>>> I've got a question about the variable bi-harmonic viscosities.
>>>>>>> In the
>>>>>>> MITgcm, there is the possibility to specify a non-dimensional
>>>>>>> horizontal
>>>>>>> bi-harmonic viscosity through viscA4Grid. The model then computes
>>>>>>> the actual viscosity taking into account the local size of the
>>>>>>> grid-point.
>>>>>>> This local grid-scale can be computed 2 different ways (through the
>>>>>>> flag useAreaViscLength):
>>>>>>>
>>>>>>> C These are (powers of) length scales
>>>>>>> IF (useAreaViscLength) THEN
>>>>>>> L2=rA(i,j,bi,bj)
>>>>>>> L4rdt=0.03125 _d 0*recip_dt*L2**2
>>>>>>> ELSE
>>>>>>> L2=2. _d
>>>>>>> 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))
>>>>>>> L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
>>>>>>> & +recip_DYF(I,J,bi,bj)**4)
>>>>>>> & +8. _d 0*((recip_DXF(I,J,bi,bj)
>>>>>>> & *recip_DYF(I,J,bi,bj))**2) )
>>>>>>> ENDIF
>>>>>>>
>>>>>>> The default (useAreaViscLength=.FALSE.) compute the grid
>>>>>>> factor directly from the length scales, otherwise it is done from
>>>>>>> the area of the grid-point (for the cube-sphere for example).
>>>>>>>
>>>>>>> I was expecting that the 2 formula above would give the same
>>>>>>> results
>>>>>>> for cartesian square grid-box (dx=dy). This is the case for L2
>>>>>>> but not for
>>>>>>> L4rdt.
>>>>>>> If useAreaViscLength is TRUE, then L4rdt = dx^4/(32*deltaT),
>>>>>>> otherwise it is L4rdt = dx^4/(20*deltaT).
>>>>>>>
>>>>>>> Does anyone know which coefficient is the correct one ?
>>>>>>> Thanks,
>>>>>>> david
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> 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