[MITgcm-devel] about dst3 (advection scheme 30)
Samar Khatiwala
spk at ldeo.columbia.edu
Sun Apr 2 12:15:28 EDT 2006
I agree!
Take a look at the email exchange on Jan 11, 2004 on the mitgcm-
support list.
I believe the signs in my version are correct. I faintly recall
checking them against the derivation
in Hundsdorfer.
Samar
On Apr 2, 2006, at 11:57 AM, Martin Losch wrote:
> 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
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel
More information about the MITgcm-devel
mailing list