| 2 |
% $Name$ |
% $Name$ |
| 3 |
|
|
| 4 |
\section{Tracer equations} |
\section{Tracer equations} |
| 5 |
|
\label{sect:tracer_equations} |
| 6 |
|
|
| 7 |
The basic discretization used for the tracer equations is the second |
The basic discretization used for the tracer equations is the second |
| 8 |
order piece-wise constant finite volume form of the forced |
order piece-wise constant finite volume form of the forced |
| 9 |
advection-diussion equations. There are many alternatives to second |
advection-diffusion equations. There are many alternatives to second |
| 10 |
order method for advection and alternative parameterizations for the |
order method for advection and alternative parameterizations for the |
| 11 |
sub-grid scale processes. The Gent-McWilliams eddy parameterization, |
sub-grid scale processes. The Gent-McWilliams eddy parameterization, |
| 12 |
KPP mixing scheme and PV flux parameterization are all dealt with in |
KPP mixing scheme and PV flux parameterization are all dealt with in |
| 14 |
part of the tracer equations and the various advection schemes will be |
part of the tracer equations and the various advection schemes will be |
| 15 |
described here. |
described here. |
| 16 |
|
|
| 17 |
|
\subsection{Time-stepping of tracers: ABII} |
| 18 |
|
\label{sect:tracer_equations_abII} |
| 19 |
|
|
| 20 |
|
The default advection scheme is the centered second order method which |
| 21 |
|
requires a second order or quasi-second order time-stepping scheme to |
| 22 |
|
be stable. Historically this has been the quasi-second order |
| 23 |
|
Adams-Bashforth method (ABII) and applied to all terms. For an |
| 24 |
|
arbitrary tracer, $\tau$, the forced advection-diffusion equation |
| 25 |
|
reads: |
| 26 |
|
\begin{equation} |
| 27 |
|
\partial_t \tau + G_{adv}^\tau = G_{diff}^\tau + G_{forc}^\tau |
| 28 |
|
\end{equation} |
| 29 |
|
where $G_{adv}^\tau$, $G_{diff}^\tau$ and $G_{forc}^\tau$ are the |
| 30 |
|
tendencies due to advection, diffusion and forcing, respectively, |
| 31 |
|
namely: |
| 32 |
|
\begin{eqnarray} |
| 33 |
|
G_{adv}^\tau & = & \partial_x u \tau + \partial_y v \tau + \partial_r w \tau |
| 34 |
|
- \tau \nabla \cdot {\bf v} \\ |
| 35 |
|
G_{diff}^\tau & = & \nabla \cdot {\bf K} \nabla \tau |
| 36 |
|
\end{eqnarray} |
| 37 |
|
and the forcing can be some arbitrary function of state, time and |
| 38 |
|
space. |
| 39 |
|
|
| 40 |
|
The term, $\tau \nabla \cdot {\bf v}$, is required to retain local |
| 41 |
|
conservation in conjunction with the linear implicit free-surface. It |
| 42 |
|
only affects the surface layer since the flow is non-divergent |
| 43 |
|
everywhere else. This term is therefore referred to as the surface |
| 44 |
|
correction term. Global conservation is not possible using the |
| 45 |
|
flux-form (as here) and a linearized free-surface |
| 46 |
|
(\cite{griffies:00,campin:02}). |
| 47 |
|
|
| 48 |
|
The continuity equation can be recovered by setting |
| 49 |
|
$G_{diff}=G_{forc}=0$ and $\tau=1$. |
| 50 |
|
|
| 51 |
|
The driver routine that calls the routines to calculate tendencies are |
| 52 |
|
{\em S/R CALC\_GT} and {\em S/R CALC\_GS} for temperature and salt |
| 53 |
|
(moisture), respectively. These in turn call a generic advection |
| 54 |
|
diffusion routine {\em S/R GAD\_CALC\_RHS} that is called with the |
| 55 |
|
flow field and relevant tracer as arguments and returns the collective |
| 56 |
|
tendency due to advection and diffusion. Forcing is add subsequently |
| 57 |
|
in {\em S/R CALC\_GT} or {\em S/R CALC\_GS} to the same tendency |
| 58 |
|
array. |
| 59 |
|
|
| 60 |
|
\fbox{ \begin{minipage}{4.75in} |
| 61 |
|
{\em S/R GAD\_CALC\_RHS} ({\em pkg/generic\_advdiff/gad\_calc\_rhs.F}) |
| 62 |
|
|
| 63 |
|
$\tau$: {\bf tracer} (argument) |
| 64 |
|
|
| 65 |
|
$G^{(n)}$: {\bf gTracer} (argument) |
| 66 |
|
|
| 67 |
|
$F_r$: {\bf fVerT} (argument) |
| 68 |
|
|
| 69 |
|
\end{minipage} } |
| 70 |
|
|
| 71 |
|
The space and time discretization are treated separately (method of |
| 72 |
|
lines). Tendencies are calculated at time levels $n$ and $n-1$ and |
| 73 |
|
extrapolated to $n+1/2$ using the Adams-Bashforth method: |
| 74 |
|
\marginpar{$\epsilon$: {\bf AB\_eps}} |
| 75 |
|
\begin{equation} |
| 76 |
|
G^{(n+1/2)} = |
| 77 |
|
(\frac{3}{2} + \epsilon) G^{(n)} - (\frac{1}{2} + \epsilon) G^{(n-1)} |
| 78 |
|
\end{equation} |
| 79 |
|
where $G^{(n)} = G_{adv}^\tau + G_{diff}^\tau + G_{src}^\tau$ at time |
| 80 |
|
step $n$. The tendency at $n-1$ is not re-calculated but rather the |
| 81 |
|
tendency at $n$ is stored in a global array for later re-use. |
| 82 |
|
|
| 83 |
|
\fbox{ \begin{minipage}{4.75in} |
| 84 |
|
{\em S/R ADAMS\_BASHFORTH2} ({\em model/src/adams\_bashforth2.F}) |
| 85 |
|
|
| 86 |
|
$G^{(n+1/2)}$: {\bf gTracer} (argument on exit) |
| 87 |
|
|
| 88 |
|
$G^{(n)}$: {\bf gTracer} (argument on entry) |
| 89 |
|
|
| 90 |
|
$G^{(n-1)}$: {\bf gTrNm1} (argument) |
| 91 |
|
|
| 92 |
|
$\epsilon$: {\bf ABeps} (PARAMS.h) |
| 93 |
|
|
| 94 |
|
\end{minipage} } |
| 95 |
|
|
| 96 |
|
The tracers are stepped forward in time using the extrapolated tendency: |
| 97 |
|
\begin{equation} |
| 98 |
|
\tau^{(n+1)} = \tau^{(n)} + \Delta t G^{(n+1/2)} |
| 99 |
|
\end{equation} |
| 100 |
|
\marginpar{$\Delta t$: {\bf deltaTtracer}} |
| 101 |
|
|
| 102 |
|
\fbox{ \begin{minipage}{4.75in} |
| 103 |
|
{\em S/R TIMESTEP\_TRACER} ({\em model/src/timestep\_tracer.F}) |
| 104 |
|
|
| 105 |
|
$\tau^{(n+1)}$: {\bf gTracer} (argument on exit) |
| 106 |
|
|
| 107 |
|
$\tau^{(n)}$: {\bf tracer} (argument on entry) |
| 108 |
|
|
| 109 |
|
$G^{(n+1/2)}$: {\bf gTracer} (argument) |
| 110 |
|
|
| 111 |
|
$\Delta t$: {\bf deltaTtracer} (PARAMS.h) |
| 112 |
|
|
| 113 |
|
\end{minipage} } |
| 114 |
|
|
| 115 |
|
Strictly speaking the ABII scheme should be applied only to the |
| 116 |
|
advection terms. However, this scheme is only used in conjunction with |
| 117 |
|
the standard second, third and fourth order advection |
| 118 |
|
schemes. Selection of any other advection scheme disables |
| 119 |
|
Adams-Bashforth for tracers so that explicit diffusion and forcing use |
| 120 |
|
the forward method. |
| 121 |
|
|
| 122 |
|
|
| 123 |
|
|
| 124 |
|
|
| 125 |
|
\section{Linear advection schemes} |
| 126 |
|
\label{sect:tracer-advection} |
| 127 |
|
|
| 128 |
|
\begin{figure} |
| 129 |
|
\resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-lo.eps}} |
| 130 |
|
\caption{ |
| 131 |
|
Comparison of 1-D advection schemes. Courant number is 0.05 with 60 |
| 132 |
|
points and solutions are shown for T=1 (one complete period). |
| 133 |
|
a) Shows the upwind biased schemes; first order upwind, DST3, |
| 134 |
|
third order upwind and second order upwind. |
| 135 |
|
b) Shows the centered schemes; Lax-Wendroff, DST4, centered second order, |
| 136 |
|
centered fourth order and finite volume fourth order. |
| 137 |
|
c) Shows the second order flux limiters: minmod, Superbee, |
| 138 |
|
MC limiter and the van Leer limiter. |
| 139 |
|
d) Shows the DST3 method with flux limiters due to Sweby with |
| 140 |
|
$\mu=1$, $\mu=c/(1-c)$ and a fourth order DST method with Sweby limiter, |
| 141 |
|
$\mu=c/(1-c)$. |
| 142 |
|
\label{fig:advect-1d-lo} |
| 143 |
|
} |
| 144 |
|
\end{figure} |
| 145 |
|
|
| 146 |
|
\begin{figure} |
| 147 |
|
\resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-hi.eps}} |
| 148 |
|
\caption{ |
| 149 |
|
Comparison of 1-D advection schemes. Courant number is 0.89 with 60 |
| 150 |
|
points and solutions are shown for T=1 (one complete period). |
| 151 |
|
a) Shows the upwind biased schemes; first order upwind and DST3. |
| 152 |
|
Third order upwind and second order upwind are unstable at this Courant number. |
| 153 |
|
b) Shows the centered schemes; Lax-Wendroff, DST4. Centered second order, |
| 154 |
|
centered fourth order and finite volume fourth order and unstable at this |
| 155 |
|
Courant number. |
| 156 |
|
c) Shows the second order flux limiters: minmod, Superbee, |
| 157 |
|
MC limiter and the van Leer limiter. |
| 158 |
|
d) Shows the DST3 method with flux limiters due to Sweby with |
| 159 |
|
$\mu=1$, $\mu=c/(1-c)$ and a fourth order DST method with Sweby limiter, |
| 160 |
|
$\mu=c/(1-c)$. |
| 161 |
|
\label{fig:advect-1d-hi} |
| 162 |
|
} |
| 163 |
|
\end{figure} |
| 164 |
|
|
| 165 |
|
The advection schemes known as centered second order, centered fourth |
| 166 |
|
order, first order upwind and upwind biased third order are known as |
| 167 |
|
linear advection schemes because the coefficient for interpolation of |
| 168 |
|
the advected tracer are linear and a function only of the flow, not |
| 169 |
|
the tracer field it self. We discuss these first since they are most |
| 170 |
|
commonly used in the field and most familiar. |
| 171 |
|
|
| 172 |
\subsection{Centered second order advection-diffusion} |
\subsection{Centered second order advection-diffusion} |
| 173 |
|
|
| 174 |
The basic discretization, centered second order, is the default. It is |
The basic discretization, centered second order, is the default. It is |
| 175 |
designed to be consistant with the continuity equation to facilitate |
designed to be consistent with the continuity equation to facilitate |
| 176 |
conservation properties analogous to the continuum: |
conservation properties analogous to the continuum. However, centered |
| 177 |
|
second order advection is notoriously noisy and must be used in |
| 178 |
|
conjunction with some finite amount of diffusion to produce a sensible |
| 179 |
|
solution. |
| 180 |
|
|
| 181 |
|
The advection operator is discretized: |
| 182 |
\begin{equation} |
\begin{equation} |
| 183 |
{\cal A}_c \Delta r_f h_c \partial_\theta |
{\cal A}_c \Delta r_f h_c G_{adv}^\tau = |
| 184 |
+ \delta_i F_x |
\delta_i F_x + \delta_j F_y + \delta_k F_r |
|
+ \delta_j F_y |
|
|
+ \delta_k F_r |
|
|
= {\cal A}_c \Delta r_f h_c {\cal S}_\theta |
|
|
+ \theta {\cal A}_c \delta_k (P-E)_{r=0} |
|
| 185 |
\end{equation} |
\end{equation} |
| 186 |
where the area integrated fluxes are given by: |
where the area integrated fluxes are given by: |
| 187 |
\begin{eqnarray} |
\begin{eqnarray} |
| 188 |
F_x & = & U \overline{ \theta }^i |
F_x & = & U \overline{ \tau }^i \\ |
| 189 |
- \kappa_h \frac{\Delta y_g \Delta r_f h_w}{\Delta x_c} \delta_i \theta \\ |
F_y & = & V \overline{ \tau }^j \\ |
| 190 |
F_y & = & V \overline{ \theta }^j |
F_r & = & W \overline{ \tau }^k |
|
- \kappa_h \frac{\Delta x_g \Delta r_f h_s}{\Delta y_c} \delta_j \theta \\ |
|
|
F_r & = & W \overline{ \theta }^k |
|
|
- \kappa_v \frac{\Delta x_g \Delta y_g}{\Delta r_c} \delta_k \theta |
|
| 191 |
\end{eqnarray} |
\end{eqnarray} |
| 192 |
The quantities $U$, $V$ and $W$ are volume fluxes defined: |
The quantities $U$, $V$ and $W$ are volume fluxes defined: |
| 193 |
\marginpar{$U$: {\bf uTrans} } |
\marginpar{$U$: {\bf uTrans} } |
| 198 |
V & = & \Delta x_g \Delta r_f h_s v \\ |
V & = & \Delta x_g \Delta r_f h_s v \\ |
| 199 |
W & = & {\cal A}_c w |
W & = & {\cal A}_c w |
| 200 |
\end{eqnarray} |
\end{eqnarray} |
|
${\cal S}$ represents the ``parameterized'' SGS processes and physics |
|
|
and sources associated with the tracer. For instance, potential |
|
|
temperature equation in the ocean has is forced by surface and |
|
|
partially penetrating heat fluxes: |
|
|
\begin{equation} |
|
|
{\cal A}_c \Delta r_f h_c {\cal S}_\theta = \frac{1}{c_p \rho_o} \delta_k {\cal A}_c {\cal Q} |
|
|
\end{equation} |
|
|
while the salt equation has no real sources, ${\cal S}=0$, which |
|
|
leaves just the $P-E$ term. |
|
|
|
|
|
The continuity equation can be recovered by setting ${\cal Q}=0$, $\kappa_h = \kappa_v = 0$ and |
|
|
$\theta=1$. The term $\theta (P-E)_{r=0}$ is required to retain local |
|
|
conservation of $\theta$. Global conservation is not possible using |
|
|
the flux-form (as here) and a linearized free-surface |
|
|
(\cite{Griffies00,Campin02}). |
|
| 201 |
|
|
| 202 |
|
For non-divergent flow, this discretization can be shown to conserve |
| 203 |
|
the tracer both locally and globally and to globally conserve tracer |
| 204 |
|
variance, $\tau^2$. The proof is given in \cite{adcroft:95,adcroft:97}. |
| 205 |
|
|
| 206 |
|
\fbox{ \begin{minipage}{4.75in} |
| 207 |
|
{\em S/R GAD\_C2\_ADV\_X} ({\em gad\_c2\_adv\_x.F}) |
| 208 |
|
|
| 209 |
|
$F_x$: {\bf uT} (argument) |
| 210 |
|
|
| 211 |
|
$U$: {\bf uTrans} (argument) |
| 212 |
|
|
| 213 |
|
$\tau$: {\bf tracer} (argument) |
| 214 |
|
|
| 215 |
|
{\em S/R GAD\_C2\_ADV\_Y} ({\em gad\_c2\_adv\_y.F}) |
| 216 |
|
|
| 217 |
|
$F_y$: {\bf vT} (argument) |
| 218 |
|
|
| 219 |
|
$V$: {\bf vTrans} (argument) |
| 220 |
|
|
| 221 |
|
$\tau$: {\bf tracer} (argument) |
| 222 |
|
|
| 223 |
|
{\em S/R GAD\_C2\_ADV\_R} ({\em gad\_c2\_adv\_r.F}) |
| 224 |
|
|
| 225 |
|
$F_r$: {\bf wT} (argument) |
| 226 |
|
|
| 227 |
|
$W$: {\bf rTrans} (argument) |
| 228 |
|
|
| 229 |
|
$\tau$: {\bf tracer} (argument) |
| 230 |
|
|
| 231 |
|
\end{minipage} } |
| 232 |
|
|
| 233 |
|
|
| 234 |
|
\subsection{Third order upwind bias advection} |
| 235 |
|
|
| 236 |
|
Upwind biased third order advection offers a relatively good |
| 237 |
|
compromise between accuracy and smoothness. It is not a ``positive'' |
| 238 |
|
scheme meaning false extrema are permitted but the amplitude of such |
| 239 |
|
are significantly reduced over the centered second order method. |
| 240 |
|
|
| 241 |
|
The third order upwind fluxes are discretized: |
| 242 |
|
\begin{eqnarray} |
| 243 |
|
F_x & = & U \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^i |
| 244 |
|
+ \frac{1}{2} |U| \delta_i \frac{1}{6} \delta_{ii} \tau \\ |
| 245 |
|
F_y & = & V \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^j |
| 246 |
|
+ \frac{1}{2} |V| \delta_j \frac{1}{6} \delta_{jj} \tau \\ |
| 247 |
|
F_r & = & W \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^k |
| 248 |
|
+ \frac{1}{2} |W| \delta_k \frac{1}{6} \delta_{kk} \tau |
| 249 |
|
\end{eqnarray} |
| 250 |
|
|
| 251 |
|
At boundaries, $\delta_{\hat{n}} \tau$ is set to zero allowing |
| 252 |
|
$\delta_{nn}$ to be evaluated. We are currently examine the accuracy |
| 253 |
|
of this boundary condition and the effect on the solution. |
| 254 |
|
|
| 255 |
|
\fbox{ \begin{minipage}{4.75in} |
| 256 |
|
{\em S/R GAD\_U3\_ADV\_X} ({\em gad\_u3\_adv\_x.F}) |
| 257 |
|
|
| 258 |
|
$F_x$: {\bf uT} (argument) |
| 259 |
|
|
| 260 |
|
$U$: {\bf uTrans} (argument) |
| 261 |
|
|
| 262 |
|
$\tau$: {\bf tracer} (argument) |
| 263 |
|
|
| 264 |
|
{\em S/R GAD\_U3\_ADV\_Y} ({\em gad\_u3\_adv\_y.F}) |
| 265 |
|
|
| 266 |
|
$F_y$: {\bf vT} (argument) |
| 267 |
|
|
| 268 |
|
$V$: {\bf vTrans} (argument) |
| 269 |
|
|
| 270 |
|
$\tau$: {\bf tracer} (argument) |
| 271 |
|
|
| 272 |
|
{\em S/R GAD\_U3\_ADV\_R} ({\em gad\_u3\_adv\_r.F}) |
| 273 |
|
|
| 274 |
|
$F_r$: {\bf wT} (argument) |
| 275 |
|
|
| 276 |
|
$W$: {\bf rTrans} (argument) |
| 277 |
|
|
| 278 |
|
$\tau$: {\bf tracer} (argument) |
| 279 |
|
|
| 280 |
|
\end{minipage} } |
| 281 |
|
|
| 282 |
|
\subsection{Centered fourth order advection} |
| 283 |
|
|
| 284 |
|
Centered fourth order advection is formally the most accurate scheme |
| 285 |
|
we have implemented and can be used to great effect in high resolution |
| 286 |
|
simulation where dynamical scales are well resolved. However, the |
| 287 |
|
scheme is noisy like the centered second order method and so must be |
| 288 |
|
used with some finite amount of diffusion. Bi-harmonic is recommended |
| 289 |
|
since it is more scale selective and less likely to diffuse away the |
| 290 |
|
well resolved gradient the fourth order scheme worked so hard to |
| 291 |
|
create. |
| 292 |
|
|
| 293 |
|
The centered fourth order fluxes are discretized: |
| 294 |
|
\begin{eqnarray} |
| 295 |
|
F_x & = & U \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^i \\ |
| 296 |
|
F_y & = & V \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^j \\ |
| 297 |
|
F_r & = & W \overline{\tau - \frac{1}{6} \delta_{ii} \tau}^k |
| 298 |
|
\end{eqnarray} |
| 299 |
|
|
| 300 |
|
As for the third order scheme, the best discretization near boundaries |
| 301 |
|
is under investigation but currently $\delta_i \tau=0$ on a boundary. |
| 302 |
|
|
| 303 |
|
\fbox{ \begin{minipage}{4.75in} |
| 304 |
|
{\em S/R GAD\_C4\_ADV\_X} ({\em gad\_c4\_adv\_x.F}) |
| 305 |
|
|
| 306 |
|
$F_x$: {\bf uT} (argument) |
| 307 |
|
|
| 308 |
|
$U$: {\bf uTrans} (argument) |
| 309 |
|
|
| 310 |
|
$\tau$: {\bf tracer} (argument) |
| 311 |
|
|
| 312 |
|
{\em S/R GAD\_C4\_ADV\_Y} ({\em gad\_c4\_adv\_y.F}) |
| 313 |
|
|
| 314 |
|
$F_y$: {\bf vT} (argument) |
| 315 |
|
|
| 316 |
|
$V$: {\bf vTrans} (argument) |
| 317 |
|
|
| 318 |
|
$\tau$: {\bf tracer} (argument) |
| 319 |
|
|
| 320 |
|
{\em S/R GAD\_C4\_ADV\_R} ({\em gad\_c4\_adv\_r.F}) |
| 321 |
|
|
| 322 |
|
$F_r$: {\bf wT} (argument) |
| 323 |
|
|
| 324 |
|
$W$: {\bf rTrans} (argument) |
| 325 |
|
|
| 326 |
|
$\tau$: {\bf tracer} (argument) |
| 327 |
|
|
| 328 |
|
\end{minipage} } |
| 329 |
|
|
| 330 |
|
|
| 331 |
|
\subsection{First order upwind advection} |
| 332 |
|
|
| 333 |
|
Although the upwind scheme is the underlying scheme for the robust or |
| 334 |
|
non-linear methods given later, we haven't actually supplied this |
| 335 |
|
method for general use. It would be very diffusive and it is unlikely |
| 336 |
|
that it could ever produce more useful results than the positive |
| 337 |
|
higher order schemes. |
| 338 |
|
|
| 339 |
|
Upwind bias is introduced into many schemes using the {\em abs} |
| 340 |
|
function and is allows the first order upwind flux to be written: |
| 341 |
|
\begin{eqnarray} |
| 342 |
|
F_x & = & U \overline{ \tau }^i - \frac{1}{2} |U| \delta_i \tau \\ |
| 343 |
|
F_y & = & V \overline{ \tau }^j - \frac{1}{2} |V| \delta_j \tau \\ |
| 344 |
|
F_r & = & W \overline{ \tau }^k - \frac{1}{2} |W| \delta_k \tau |
| 345 |
|
\end{eqnarray} |
| 346 |
|
|
| 347 |
|
If for some reason, the above method is required, then the second |
| 348 |
|
order flux limiter scheme described later reduces to the above scheme |
| 349 |
|
if the limiter is set to zero. |
| 350 |
|
|
| 351 |
|
|
| 352 |
|
\section{Non-linear advection schemes} |
| 353 |
|
|
| 354 |
|
Non-linear advection schemes invoke non-linear interpolation and are |
| 355 |
|
widely used in computational fluid dynamics (non-linear does not refer |
| 356 |
|
to the non-linearity of the advection operator). The flux limited |
| 357 |
|
advection schemes belong to the class of finite volume methods which |
| 358 |
|
neatly ties into the spatial discretization of the model. |
| 359 |
|
|
| 360 |
|
When employing the flux limited schemes, first order upwind or |
| 361 |
|
direct-space-time method the time-stepping is switched to forward in |
| 362 |
|
time. |
| 363 |
|
|
| 364 |
|
\subsection{Second order flux limiters} |
| 365 |
|
|
| 366 |
|
The second order flux limiter method can be cast in several ways but |
| 367 |
|
is generally expressed in terms of other flux approximations. For |
| 368 |
|
example, in terms of a first order upwind flux and second order |
| 369 |
|
Lax-Wendroff flux, the limited flux is given as: |
| 370 |
|
\begin{equation} |
| 371 |
|
F = F_1 + \psi(r) F_{LW} |
| 372 |
|
\end{equation} |
| 373 |
|
where $\psi(r)$ is the limiter function, |
| 374 |
|
\begin{equation} |
| 375 |
|
F_1 = u \overline{\tau}^i - \frac{1}{2} |u| \delta_i \tau |
| 376 |
|
\end{equation} |
| 377 |
|
is the upwind flux, |
| 378 |
|
\begin{equation} |
| 379 |
|
F_{LW} = F_1 + \frac{|u|}{2} (1-c) \delta_i \tau |
| 380 |
|
\end{equation} |
| 381 |
|
is the Lax-Wendroff flux and $c = \frac{u \Delta t}{\Delta x}$ is the |
| 382 |
|
Courant (CFL) number. |
| 383 |
|
|
| 384 |
|
The limiter function, $\psi(r)$, takes the slope ratio |
| 385 |
|
\begin{eqnarray} |
| 386 |
|
r = \frac{ \tau_{i-1} - \tau_{i-2} }{ \tau_{i} - \tau_{i-1} } & \forall & u > 0 |
| 387 |
|
\\ |
| 388 |
|
r = \frac{ \tau_{i+1} - \tau_{i} }{ \tau_{i} - \tau_{i-1} } & \forall & u < 0 |
| 389 |
|
\end{eqnarray} |
| 390 |
|
as it's argument. There are many choices of limiter function but we |
| 391 |
|
only provide the Superbee limiter \cite{roe:85}: |
| 392 |
|
\begin{equation} |
| 393 |
|
\psi(r) = \max[0,\min[1,2r],\min[2,r]] |
| 394 |
|
\end{equation} |
| 395 |
|
|
| 396 |
|
\fbox{ \begin{minipage}{4.75in} |
| 397 |
|
{\em S/R GAD\_FLUXLIMIT\_ADV\_X} ({\em gad\_fluxlimit\_adv\_x.F}) |
| 398 |
|
|
| 399 |
|
$F_x$: {\bf uT} (argument) |
| 400 |
|
|
| 401 |
|
$U$: {\bf uTrans} (argument) |
| 402 |
|
|
| 403 |
|
$\tau$: {\bf tracer} (argument) |
| 404 |
|
|
| 405 |
|
{\em S/R GAD\_FLUXLIMIT\_ADV\_Y} ({\em gad\_fluxlimit\_adv\_y.F}) |
| 406 |
|
|
| 407 |
|
$F_y$: {\bf vT} (argument) |
| 408 |
|
|
| 409 |
|
$V$: {\bf vTrans} (argument) |
| 410 |
|
|
| 411 |
|
$\tau$: {\bf tracer} (argument) |
| 412 |
|
|
| 413 |
|
{\em S/R GAD\_FLUXLIMIT\_ADV\_R} ({\em gad\_fluxlimit\_adv\_r.F}) |
| 414 |
|
|
| 415 |
|
$F_r$: {\bf wT} (argument) |
| 416 |
|
|
| 417 |
|
$W$: {\bf rTrans} (argument) |
| 418 |
|
|
| 419 |
|
$\tau$: {\bf tracer} (argument) |
| 420 |
|
|
| 421 |
|
\end{minipage} } |
| 422 |
|
|
| 423 |
|
|
| 424 |
|
\subsection{Third order direct space time} |
| 425 |
|
|
| 426 |
|
The direct-space-time method deals with space and time discretization |
| 427 |
|
together (other methods that treat space and time separately are known |
| 428 |
|
collectively as the ``Method of Lines''). The Lax-Wendroff scheme |
| 429 |
|
falls into this category; it adds sufficient diffusion to a second |
| 430 |
|
order flux that the forward-in-time method is stable. The upwind |
| 431 |
|
biased third order DST scheme is: |
| 432 |
|
\begin{eqnarray} |
| 433 |
|
F = u \left( \tau_{i-1} |
| 434 |
|
+ d_0 (\tau_{i}-\tau_{i-1}) + d_1 (\tau_{i-1}-\tau_{i-2}) \right) |
| 435 |
|
& \forall & u > 0 \\ |
| 436 |
|
F = u \left( \tau_{i} |
| 437 |
|
- d_0 (\tau_{i}-\tau_{i-1}) - d_1 (\tau_{i+1}-\tau_{i}) \right) |
| 438 |
|
& \forall & u < 0 |
| 439 |
|
\end{eqnarray} |
| 440 |
|
where |
| 441 |
|
\begin{eqnarray} |
| 442 |
|
d_1 & = & \frac{1}{6} ( 2 - |c| ) ( 1 - |c| ) \\ |
| 443 |
|
d_2 & = & \frac{1}{6} ( 1 - |c| ) ( 1 + |c| ) |
| 444 |
|
\end{eqnarray} |
| 445 |
|
The coefficients $d_0$ and $d_1$ approach $1/3$ and $1/6$ respectively |
| 446 |
|
as the Courant number, $c$, vanishes. In this limit, the conventional |
| 447 |
|
third order upwind method is recovered. For finite Courant number, the |
| 448 |
|
deviations from the linear method are analogous to the diffusion added |
| 449 |
|
to centered second order advection in the Lax-Wendroff scheme. |
| 450 |
|
|
| 451 |
|
The DST3 method described above must be used in a forward-in-time |
| 452 |
|
manner and is stable for $0 \le |c| \le 1$. Although the scheme |
| 453 |
|
appears to be forward-in-time, it is in fact second order in time and |
| 454 |
|
the accuracy increases with the Courant number! For low Courant |
| 455 |
|
number, DST3 produces very similar results (indistinguishable in |
| 456 |
|
Fig.~\ref{fig:advect-1d-lo}) to the linear third order method but for |
| 457 |
|
large Courant number, where the linear upwind third order method is |
| 458 |
|
unstable, the scheme is extremely accurate |
| 459 |
|
(Fig.~\ref{fig:advect-1d-hi}) with only minor overshoots. |
| 460 |
|
|
| 461 |
|
\fbox{ \begin{minipage}{4.75in} |
| 462 |
|
{\em S/R GAD\_DST3\_ADV\_X} ({\em gad\_dst3\_adv\_x.F}) |
| 463 |
|
|
| 464 |
|
$F_x$: {\bf uT} (argument) |
| 465 |
|
|
| 466 |
|
$U$: {\bf uTrans} (argument) |
| 467 |
|
|
| 468 |
|
$\tau$: {\bf tracer} (argument) |
| 469 |
|
|
| 470 |
|
{\em S/R GAD\_DST3\_ADV\_Y} ({\em gad\_dst3\_adv\_y.F}) |
| 471 |
|
|
| 472 |
|
$F_y$: {\bf vT} (argument) |
| 473 |
|
|
| 474 |
|
$V$: {\bf vTrans} (argument) |
| 475 |
|
|
| 476 |
|
$\tau$: {\bf tracer} (argument) |
| 477 |
|
|
| 478 |
|
{\em S/R GAD\_DST3\_ADV\_R} ({\em gad\_dst3\_adv\_r.F}) |
| 479 |
|
|
| 480 |
|
$F_r$: {\bf wT} (argument) |
| 481 |
|
|
| 482 |
|
$W$: {\bf rTrans} (argument) |
| 483 |
|
|
| 484 |
|
$\tau$: {\bf tracer} (argument) |
| 485 |
|
|
| 486 |
|
\end{minipage} } |
| 487 |
|
|
| 488 |
|
|
| 489 |
|
\subsection{Third order direct space time with flux limiting} |
| 490 |
|
|
| 491 |
|
The overshoots in the DST3 method can be controlled with a flux limiter. |
| 492 |
|
The limited flux is written: |
| 493 |
|
\begin{equation} |
| 494 |
|
F = |
| 495 |
|
\frac{1}{2}(u+|u|)\left( \tau_{i-1} + \psi(r^+)(\tau_{i} - \tau_{i-1} )\right) |
| 496 |
|
+ |
| 497 |
|
\frac{1}{2}(u-|u|)\left( \tau_{i-1} + \psi(r^-)(\tau_{i} - \tau_{i-1} )\right) |
| 498 |
|
\end{equation} |
| 499 |
|
where |
| 500 |
|
\begin{eqnarray} |
| 501 |
|
r^+ & = & \frac{\tau_{i-1} - \tau_{i-2}}{\tau_{i} - \tau_{i-1}} \\ |
| 502 |
|
r^- & = & \frac{\tau_{i+1} - \tau_{i}}{\tau_{i} - \tau_{i-1}} |
| 503 |
|
\end{eqnarray} |
| 504 |
|
and the limiter is the Sweby limiter: |
| 505 |
|
\begin{equation} |
| 506 |
|
\psi(r) = \max[0, \min[\min(1,d_0+d_1r],\frac{1-c}{c}r ]] |
| 507 |
|
\end{equation} |
| 508 |
|
|
| 509 |
|
\fbox{ \begin{minipage}{4.75in} |
| 510 |
|
{\em S/R GAD\_DST3FL\_ADV\_X} ({\em gad\_dst3\_adv\_x.F}) |
| 511 |
|
|
| 512 |
|
$F_x$: {\bf uT} (argument) |
| 513 |
|
|
| 514 |
|
$U$: {\bf uTrans} (argument) |
| 515 |
|
|
| 516 |
|
$\tau$: {\bf tracer} (argument) |
| 517 |
|
|
| 518 |
|
{\em S/R GAD\_DST3FL\_ADV\_Y} ({\em gad\_dst3\_adv\_y.F}) |
| 519 |
|
|
| 520 |
|
$F_y$: {\bf vT} (argument) |
| 521 |
|
|
| 522 |
|
$V$: {\bf vTrans} (argument) |
| 523 |
|
|
| 524 |
|
$\tau$: {\bf tracer} (argument) |
| 525 |
|
|
| 526 |
|
{\em S/R GAD\_DST3FL\_ADV\_R} ({\em gad\_dst3\_adv\_r.F}) |
| 527 |
|
|
| 528 |
|
$F_r$: {\bf wT} (argument) |
| 529 |
|
|
| 530 |
|
$W$: {\bf rTrans} (argument) |
| 531 |
|
|
| 532 |
|
$\tau$: {\bf tracer} (argument) |
| 533 |
|
|
| 534 |
|
\end{minipage} } |
| 535 |
|
|
| 536 |
|
|
| 537 |
|
\subsection{Multi-dimensional advection} |
| 538 |
|
|
| 539 |
|
\begin{figure} |
| 540 |
|
\resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-lo-diag.eps}} |
| 541 |
|
\caption{ |
| 542 |
|
Comparison of advection schemes in two dimensions; diagonal advection |
| 543 |
|
of a resolved Gaussian feature. Courant number is 0.01 with |
| 544 |
|
30$\times$30 points and solutions are shown for T=1/2. White lines |
| 545 |
|
indicate zero crossing (ie. the presence of false minima). The left |
| 546 |
|
column shows the second order schemes; top) centered second order with |
| 547 |
|
Adams-Bashforth, middle) Lax-Wendroff and bottom) Superbee flux |
| 548 |
|
limited. The middle column shows the third order schemes; top) upwind |
| 549 |
|
biased third order with Adams-Bashforth, middle) third order direct |
| 550 |
|
space-time method and bottom) the same with flux limiting. The top |
| 551 |
|
right panel shows the centered fourth order scheme with |
| 552 |
|
Adams-Bashforth and right middle panel shows a fourth order variant on |
| 553 |
|
the DST method. Bottom right panel shows the Superbee flux limiter |
| 554 |
|
(second order) applied independently in each direction (method of |
| 555 |
|
lines). |
| 556 |
|
\label{fig:advect-2d-lo-diag} |
| 557 |
|
} |
| 558 |
|
\end{figure} |
| 559 |
|
|
| 560 |
|
\begin{figure} |
| 561 |
|
\resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-mid-diag.eps}} |
| 562 |
|
\caption{ |
| 563 |
|
Comparison of advection schemes in two dimensions; diagonal advection |
| 564 |
|
of a resolved Gaussian feature. Courant number is 0.27 with |
| 565 |
|
30$\times$30 points and solutions are shown for T=1/2. White lines |
| 566 |
|
indicate zero crossing (ie. the presence of false minima). The left |
| 567 |
|
column shows the second order schemes; top) centered second order with |
| 568 |
|
Adams-Bashforth, middle) Lax-Wendroff and bottom) Superbee flux |
| 569 |
|
limited. The middle column shows the third order schemes; top) upwind |
| 570 |
|
biased third order with Adams-Bashforth, middle) third order direct |
| 571 |
|
space-time method and bottom) the same with flux limiting. The top |
| 572 |
|
right panel shows the centered fourth order scheme with |
| 573 |
|
Adams-Bashforth and right middle panel shows a fourth order variant on |
| 574 |
|
the DST method. Bottom right panel shows the Superbee flux limiter |
| 575 |
|
(second order) applied independently in each direction (method of |
| 576 |
|
lines). |
| 577 |
|
\label{fig:advect-2d-mid-diag} |
| 578 |
|
} |
| 579 |
|
\end{figure} |
| 580 |
|
|
| 581 |
|
\begin{figure} |
| 582 |
|
\resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-hi-diag.eps}} |
| 583 |
|
\caption{ |
| 584 |
|
Comparison of advection schemes in two dimensions; diagonal advection |
| 585 |
|
of a resolved Gaussian feature. Courant number is 0.47 with |
| 586 |
|
30$\times$30 points and solutions are shown for T=1/2. White lines |
| 587 |
|
indicate zero crossings and initial maximum values (ie. the presence |
| 588 |
|
of false extrema). The left column shows the second order schemes; |
| 589 |
|
top) centered second order with Adams-Bashforth, middle) Lax-Wendroff |
| 590 |
|
and bottom) Superbee flux limited. The middle column shows the third |
| 591 |
|
order schemes; top) upwind biased third order with Adams-Bashforth, |
| 592 |
|
middle) third order direct space-time method and bottom) the same with |
| 593 |
|
flux limiting. The top right panel shows the centered fourth order |
| 594 |
|
scheme with Adams-Bashforth and right middle panel shows a fourth |
| 595 |
|
order variant on the DST method. Bottom right panel shows the Superbee |
| 596 |
|
flux limiter (second order) applied independently in each direction |
| 597 |
|
(method of lines). |
| 598 |
|
\label{fig:advect-2d-hi-diag} |
| 599 |
|
} |
| 600 |
|
\end{figure} |
| 601 |
|
|
| 602 |
|
|
| 603 |
|
|
| 604 |
|
In many of the aforementioned advection schemes the behavior in |
| 605 |
|
multiple dimensions is not necessarily as good as the one dimensional |
| 606 |
|
behavior. For instance, a shape preserving monotonic scheme in one |
| 607 |
|
dimension can have severe shape distortion in two dimensions if the |
| 608 |
|
two components of horizontal fluxes are treated independently. There |
| 609 |
|
is a large body of literature on the subject dealing with this problem |
| 610 |
|
and among the fixes are operator and flux splitting methods, corner |
| 611 |
|
flux methods and more. We have adopted a variant on the standard |
| 612 |
|
splitting methods that allows the flux calculations to be implemented |
| 613 |
|
as if in one dimension: |
| 614 |
|
\begin{eqnarray} |
| 615 |
|
\tau^{n+1/3} & = & \tau^{n} |
| 616 |
|
- \Delta t \left( \frac{1}{\Delta x} \delta_i F^x(\tau^{n}) |
| 617 |
|
+ \tau^{n} \frac{1}{\Delta x} \delta_i u \right) \\ |
| 618 |
|
\tau^{n+2/3} & = & \tau^{n} |
| 619 |
|
- \Delta t \left( \frac{1}{\Delta y} \delta_j F^y(\tau^{n+1/3}) |
| 620 |
|
+ \tau^{n} \frac{1}{\Delta y} \delta_i v \right) \\ |
| 621 |
|
\tau^{n+3/3} & = & \tau^{n} |
| 622 |
|
- \Delta t \left( \frac{1}{\Delta r} \delta_k F^x(\tau^{n+2/3}) |
| 623 |
|
+ \tau^{n} \frac{1}{\Delta r} \delta_i w \right) |
| 624 |
|
\end{eqnarray} |
| 625 |
|
|
| 626 |
|
In order to incorporate this method into the general model algorithm, |
| 627 |
|
we compute the effective tendency rather than update the tracer so |
| 628 |
|
that other terms such as diffusion are using the $n$ time-level and |
| 629 |
|
not the updated $n+3/3$ quantities: |
| 630 |
|
\begin{equation} |
| 631 |
|
G^{n+1/2}_{adv} = \frac{1}{\Delta t} ( \tau^{n+3/3} - \tau^{n} ) |
| 632 |
|
\end{equation} |
| 633 |
|
So that the over all time-stepping looks likes: |
| 634 |
|
\begin{equation} |
| 635 |
|
\tau^{n+1} = \tau^{n} + \Delta t \left( G^{n+1/2}_{adv} + G_{diff}(\tau^{n}) + G^{n}_{forcing} \right) |
| 636 |
|
\end{equation} |
| 637 |
|
|
| 638 |
|
\fbox{ \begin{minipage}{4.75in} |
| 639 |
|
{\em S/R GAD\_ADVECTION} ({\em gad\_advection.F}) |
| 640 |
|
|
| 641 |
|
$\tau$: {\bf Tracer} (argument) |
| 642 |
|
|
| 643 |
|
$G^{n+1/2}_{adv}$: {\bf Gtracer} (argument) |
| 644 |
|
|
| 645 |
|
$F_x, F_y, F_r$: {\bf af} (local) |
| 646 |
|
|
| 647 |
|
$U$: {\bf uTrans} (local) |
| 648 |
|
|
| 649 |
|
$V$: {\bf vTrans} (local) |
| 650 |
|
|
| 651 |
|
$W$: {\bf rTrans} (local) |
| 652 |
|
|
| 653 |
|
\end{minipage} } |
| 654 |
|
|
| 655 |
|
|
| 656 |
|
\section{Comparison of advection schemes} |
| 657 |
|
|
| 658 |
|
Figs.~\ref{fig:advect-2d-lo-diag}, \ref{fig:advect-2d-mid-diag} and |
| 659 |
|
\ref{fig:advect-2d-hi-diag} show solutions to a simple diagonal |
| 660 |
|
advection problem using a selection of schemes for low, moderate and |
| 661 |
|
high Courant numbers, respectively. The top row shows the linear |
| 662 |
|
schemes, integrated with the Adams-Bashforth method. Theses schemes |
| 663 |
|
are clearly unstable for the high Courant number and weakly unstable |
| 664 |
|
for the moderate Courant number. The presence of false extrema is very |
| 665 |
|
apparent for all Courant numbers. The middle row shows solutions |
| 666 |
|
obtained with the unlimited but multi-dimensional schemes. These |
| 667 |
|
solutions also exhibit false extrema though the pattern now shows |
| 668 |
|
symmetry due to the multi-dimensional scheme. Also, the schemes are |
| 669 |
|
stable at high Courant number where the linear schemes weren't. The |
| 670 |
|
bottom row (left and middle) shows the limited schemes and most |
| 671 |
|
obvious is the absence of false extrema. The accuracy and stability of |
| 672 |
|
the unlimited non-linear schemes is retained at high Courant number |
| 673 |
|
but at low Courant number the tendency is to loose amplitude in sharp |
| 674 |
|
peaks due to diffusion. The one dimensional tests shown in |
| 675 |
|
Figs.~\ref{fig:advect-1d-lo} and \ref{fig:advect-1d-hi} showed this |
| 676 |
|
phenomenon. |
| 677 |
|
|
| 678 |
|
Finally, the bottom left and right panels use the same advection |
| 679 |
|
scheme but the right does not use the multi-dimensional method. At low |
| 680 |
|
Courant number this appears to not matter but for moderate Courant |
| 681 |
|
number severe distortion of the feature is apparent. Moreover, the |
| 682 |
|
stability of the multi-dimensional scheme is determined by the maximum |
| 683 |
|
Courant number applied of each dimension while the stability of the |
| 684 |
|
method of lines is determined by the sum. Hence, in the high Courant |
| 685 |
|
number plot, the scheme is unstable. |
| 686 |
|
|
| 687 |
|
With many advection schemes implemented in the code two questions |
| 688 |
|
arise: ``Which scheme is best?'' and ``Why don't you just offer the |
| 689 |
|
best advection scheme?''. Unfortunately, no one advection scheme is |
| 690 |
|
``the best'' for all particular applications and for new applications |
| 691 |
|
it is often a matter of trial to determine which is most |
| 692 |
|
suitable. Here are some guidelines but these are not the rule; |
| 693 |
|
\begin{itemize} |
| 694 |
|
\item If you have a coarsely resolved model, using a |
| 695 |
|
positive or upwind biased scheme will introduce significant diffusion |
| 696 |
|
to the solution and using a centered higher order scheme will |
| 697 |
|
introduce more noise. In this case, simplest may be best. |
| 698 |
|
\item If you have a high resolution model, using a higher order |
| 699 |
|
scheme will give a more accurate solution but scale-selective |
| 700 |
|
diffusion might need to be employed. The flux limited methods |
| 701 |
|
offer similar accuracy in this regime. |
| 702 |
|
\item If your solution has shocks or propagating fronts then a |
| 703 |
|
flux limited scheme is almost essential. |
| 704 |
|
\item If your time-step is limited by advection, the multi-dimensional |
| 705 |
|
non-linear schemes have the most stability (up to Courant number 1). |
| 706 |
|
\item If you need to know how much diffusion/dissipation has occurred you |
| 707 |
|
will have a lot of trouble figuring it out with a non-linear method. |
| 708 |
|
\item The presence of false extrema is non-physical and this alone is the |
| 709 |
|
strongest argument for using a positive scheme. |
| 710 |
|
\end{itemize} |
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch