[Mitgcm-support] third order flux near boundaries

[mitgcm-support at dev.mitgcm.org]

J-M,

I've work through the analysis for truncation error of the third order
schemes near a boundary. Remember that the 1/6th comes from considering
the truncation error of the evolution equation not the flux itself.
Using 1/8th in flux gives third order flux but only first order evolution.

It turns out that you *can* retain third order in the evolution equation
near the boundary. The third order upwind flux is normally:

(1) F(i) = U(i) [ -1 T(i-2) +  5 T(i-1) +  2 T(i)            ]/6

but if a boundary exists at the U(i-1) point then modifying the flux to:

(2) F(i) = U(i) [           + 94 T(i-1) + 43 T(i) + 1 T(i+1) ]/138

retains a cubic polynomial approximation for the evolution equation. The
above is rather "unnatural" and I've never seen anyone actually retain
third order accuracy near boundaries so I'm not advocating it - just 
noting that it's possible.

Instead, second order evolution can be retained near boundaries using:

(3) F(i) = U(i) [           +  2 T(i-1) +  5 T(i) - 1 T(i+1) ]/6

which is the downwind third order approximation. This is different
from anything that's been tried in the code because it's stencil
is as wide as the fourth order stencil.

Using the "K" method the third order interior flux is given as:

(4) F(i) = U(i) [ T(i-1) + (1-K)/4 { T(i-1)-T(i-2) } + (1+K)/4 { T(i) - T(i-1) }
]

with K=1/3. This is identical to (1) but by simply masking the first gradient
at the boundary gives:

(5) F(i) = U(i) [ T(i-1)                             + (1+K)/4 { T(i) - T(i-1) }
]

         = U(i) [           +  4 T(i-1) +  2 T(i)            ]/6

which is notable because 4/6 is not very different from 94/138. ie. (5) is close
to (2)
so could be considered quasi-third order accurate near boundaries.

To consider any of the discretizations we've actually tried in the model we
have to switch to the 1/8th method and consider the accuracy of the flux
(remember
we don't have high accuracy in the evolution equation for this).

The third order flux can be written F = U bar_i (T - 1/8 Tii) + 1/2 |U| del_i (
1/8 Tii )
where Tii is the second difference [1 -2 1]. For U>0 and setting Ti=0 on the
boundary gives:

(6) F(i) = U(i) [           +  5 T(i-1) +  3 T(i)             ]/8

This is structually similar to (5) but numerically further from (2) so that the
evolution equation is even less accurate. However, fitting a quadratic to T(i-1)
and
T(i) and using the boundary condition gives exactly (6) which means that the
flux (6) is still third order accurate.

Using this notation but with 1/6th: F = U bar_i (T - Tii/6 ) + 1/2 |U| del_i (
Tii/6 )
and masking Ti on the boundary is what we used to do in gad_u3_adv_*.F and this
recovers (5), the K method.

Currently, the masking is such as to reduce the flux to conventional second
order:

(7) F(i) = U(i) [           +  1 T(i-1) +  1 T(i)            ]/2

which has only first order accuracy in the evolution equation.

Chris's calc_gs.F uses (6) as far as I can tell?

I've no idea which is better. I don't recommend trying them all/any in the
model.
We'll test them in matlab to figure out which behaves best but there's quite
a few variations to work through...

A.