[MITgcm-devel] about dst3 (advection scheme 30)
Martin Losch
mlosch at awi-bremerhaven.de
Sun Apr 2 11:57:23 EDT 2006
Hi Samar,
this is your code, taken from the up-to-date version of
gad_dst3_adv_x (version 1.6):
> 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) )
and this is what is in the first version of gad_dst3_adv_x:
> 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 )
You can see that it is the same except for a SIGN ERROR! You tell me
which version is right (o: (I guess yours is correct).
In any case, in this implementation, half of the code in
gad_dst3_adv_x is not necessary and there is no division on involved.
The misunderstanding about divisions by zero appears, because for
version 1.2 the routine was rewritten to resemble gad_dst3fl_adv_x.
Then the only difference between dst3 and the flux limited version is
the flux limiter itself.
And this is exactly my point: as the current version of the flux
limited dst3 has diverged from this form anyway, there is no reason
for all of these superfluous computations of thetaM/P and psiM/P in
the gad_dst3_adv_x anymore, which make it non-linear. Removing these
would leave us with your code which is identical to the original
except for the sign error (in version/revision 1.1). Do you agree?
Martin
On Apr 2, 2006, at 5:40 PM, Samar Khatiwala wrote:
> 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
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list