[MITgcm-devel] about dst3 (advection scheme 30)

[Martin Losch]

Hi Samar,
when you have a look at the history of, e.g, gad_dst3_adv_x.F, for  
example here:
> http://dev.mitgcm.org/cgi-bin/viewcvs.cgi/MITgcm/pkg/ 
> generic_advdiff/gad_dst3_adv_x.F?only_with_tag=MAIN
you can see that the first version (1.1) of gad_dst3_adv_x.F is  
exactly what you (re-)coded (reproduced) for the your matrix scheme.  
Version 1.2 is already "non-linear" to resemble the dst3fl version,  
but now dst3fl has evolved independently so I don't see any reason  
for keeping this slightly non-linear version. But maybe there are  
reasons that I cannot see?

Martin


On Apr 2, 2006, at 4:25 PM, Samar Khatiwala wrote:

> Martin, some history here...
>
> I don't know about any "original version", but the code within  
> ifdef ALLOW_MATRIX / endif was
> something I wrote. I wanted a linear advection scheme (as DST3  
> should be), and the way the code
> was previously implemented made it nonlinear. (A bit odd for a  
> scheme that was NOT flux limiting.)
> I got tired of using my own mods every time I wanted to use DST3,  
> so Patrick kindly incorporated my
> modifications in gad_dst3_adv_* within the ALLOW_MATRIX block. As  
> currently implemented, the
> (correct) linear version of DST3 is only executed when using the  
> matrix package.
>
> Samar
>
> On Apr 2, 2006, at 7:56 AM, Martin Losch wrote:
>
>> Hi there,
>>
>> while searching for a possible bug that causes the noise in my  
>> 1/6th degree simulations, I came across a (probably historical)  
>> oddity in the DST3 routines gad_dst3_adv_x/y/r.F:
>> The original implementation was (only for the horizontal x- 
>> direction):
>>>         Rjp=(tracer(i+1,j)-tracer(i,j))*maskW(i+1,j,k,bi,bj)
>>>         Rj =(tracer(i,j)-tracer(i-1,j))*maskW(i,j,k,bi,bj)
>>>         Rjm=(tracer(i-1,j)-tracer(i-2,j))*maskW(i-1,j,k,bi,bj)
>>>
>>>         cfl=uVel(i,j,k,bi,bj)*deltaT*recip_dxc(i,j,bi,bj)
>>>         d0=(2.-abs(cfl))*(1.-abs(cfl))*oneSixth
>>>         d1=(1.-cfl)*(1.+cfl)*oneSixth
>>>         uT(i,j)=
>>>      &   0.5*(uTrans(i,j)+abs(uTrans(i,j)))
>>>      &      *( Tracer(i-1,j) + d0*Rj + d1*Rjm )
>>>      &  +0.5*(uTrans(i,j)-abs(uTrans(i,j)))
>>>      &      *( Tracer( i ,j) - d0*Rj + d1*Rjp )
>>
>> Then Alistair changed it to make it look more like the flux  
>> limited version DST3FL (advection Scheme 33):
>>>         Rjp=(tracer(i+1,j)-tracer(i,j))*maskW(i+1,j,k,bi,bj)
>>>         Rj =(tracer(i,j)-tracer(i-1,j))*maskW(i,j,k,bi,bj)
>>>         Rjm=(tracer(i-1,j)-tracer(i-2,j))*maskW(i-1,j,k,bi,bj)
>>>
>>>         cfl=abs(uVel(i,j,k,bi,bj)*deltaT*recip_dxc(i,j,bi,bj))
>>>         d0=(2.-cfl)*(1.-cfl)*oneSixth
>>>         d1=(1.-cfl*cfl)*oneSixth
>>> c       thetaP=0.
>>> c       IF (Rj.NE.0.) thetaP=Rjm/Rj
>>>         thetaP=Rjm/(1.D-20+Rj)
>>>         psiP=d0+d1*thetaP
>>> c       psiP=max(0.,min(min(1.,psiP),(1.-cfl)/(1.D-20+cfl)*thetaP))
>>>         thetaM=Rjp/(1.D-20+Rj)
>>> c       thetaM=0.
>>> c       IF (Rj.NE.0.) thetaM=Rjp/Rj
>>>         psiM=d0+d1*thetaM
>>> c       psiM=max(0.,min(min(1.,psiM),(1.-cfl)/(1.D-20+cfl)*thetaM))
>>>         uT(i,j)=
>>>      &   0.5*(uTrans(i,j)+abs(uTrans(i,j)))
>>>      &      *( Tracer(i-1,j) + psiP*Rj )
>>>      &  +0.5*(uTrans(i,j)-abs(uTrans(i,j)))
>>>      &      *( Tracer( i ,j) - psiM*Rj )
>>
>> The CVS comment is "Slight re-write to help debug flux limited form."
>> This version involves divisions by Rj, the case of Rj=0 is caught  
>> by adding a small number (1D-20), Patrick then changed this code  
>> again to make it adjoinable, but making the small number  
>> adjoinable. As the debugging of the flux limited form is probably  
>> finished, wouldn't it make more sense to revert the code to the  
>> orginal version which did not involve any division and should be  
>> adjoinable much more easily. Also, the original version is  
>> actually there, ifdef ALLOW_MATRIX. It all doesn't make too much  
>> sense to me. I would suggest reverting to the orgininal version,  
>> remove the ALLOW_MATRIX-ifdefs, and maybe add a flag such as  
>> ALLOW_DEBUGGING_DST3FL which contains the current version of the  
>> code. This modification should not affect any verification  
>> experiment, as the only one using advScheme=30 is matrix_example,  
>> where ALLOW_MATRIX is define, and even there it is "only" in the  
>> passive tracer package:
>>> csysm3::verification> grep 30, */input/data*
>>> matrix_example/input/data.ptracers: PTRACERS_advScheme(1)=30,
>>
>> What do you think about this?
>>
>> 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