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

[Martin Losch]

Hi Alistair,

I still cannot recieve emails, but I read your comments in the "devel- 
archives". Here's the clarification: in an earlier email I put the  
differences (see a copy below). The original version ( which is  
essentially the same as the version within #ifdef ALLOW_MATRIX,  
except for a sign error) does not involve ANY ratios. That's why I  
think that this version will  work much better with the adjoint.

Martin

>>>>>> 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.