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

[Samar Khatiwala]

Martin,

I implemented the code as it appears in the manual and the book by  
Hundsdorfer et al.
(with some minor changes and correcting for some mistakes in the  
earlier documentation). I brought
this up many times with JMC and Steph and was told that the code was  
written so as to avoid
division by a small number. That this made DST3 slightly nonlinear  
was not considered important.
This version still remains as the default for purely legacy reasons.  
Personally, I think the linear
version -already implemented within the ALLOW_MATRIX block - should  
be the default. I use
DST3 extensively, and I always have to comment out the nonlinear  
version.

Samar

On Apr 2, 2006, at 11:20 AM, Martin Losch wrote:

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