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

Martin Losch mlosch at awi-bremerhaven.de
Tue Apr 18 12:39:31 EDT 2006 

Hi again,

after a short break I have turned to the dst3 issue again. I would  
like to change dst3 back to the "correct" version, if nobody objects.  
The adjoint should work even better with this version, as there are  
no divisions by small numbers etc (The multidimAdvection schemes 1,  
20, 30, work for me only as long as MultiDimAdvection = .false., is  
that correct?). However, I do see the point that the adjoint needs to  
be tested. If Patrick could tell me, which test are to be done I  
could try to do them.

Martin

On Apr 2, 2006, at 11:08 PM, Jean-Michel Campin wrote:

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

More information about the MITgcm-devel mailing list