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

Sun Apr 2 17:08:34 EDT 2006

Hi Martin,

I changed the dst3fl so that it is independant of the 
magnitude of the tracer field (changing the units of 
a tracer concentration has now no effect on the solution).
And It's around the same time that I moved Samar's version of 
the dst3 S/R to pkg/generic_advdiff.

The reason for not changing the default dst3 version (and keeping 
the right version under #ifdef ALLOW_MATRIX / if useMatrix ) is
that it's used with the adjoint (33 returns to 30 in the adjoint mode)
and would need to be tested.
I remember that at that time, I talked to Patrick about this.

Cheers,
Jean-Michel

On Sun, Apr 02, 2006 at 05:20:19PM +0200, 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