/[MITgcm]/manual/s_algorithm/text/tracer.tex
ViewVC logotype

Diff of /manual/s_algorithm/text/tracer.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.3 by adcroft, Tue Sep 25 20:13:42 2001 UTC revision 1.22 by jmc, Wed Jun 28 17:01:34 2006 UTC
# Line 2  Line 2 
2  % $Name$  % $Name$
3    
4  \section{Tracer equations}  \section{Tracer equations}
5    \label{sect:tracer_equations}
6    \begin{rawhtml}
7    <!-- CMIREDIR:tracer_equations: -->
8    \end{rawhtml}
9    
10  The basic discretization used for the tracer equations is the second  The basic discretization used for the tracer equations is the second
11  order piece-wise constant finite volume form of the forced  order piece-wise constant finite volume form of the forced
12  advection-diussion equations. There are many alternatives to second  advection-diffusion equations. There are many alternatives to second
13  order method for advection and alternative parameterizations for the  order method for advection and alternative parameterizations for the
14  sub-grid scale processes. The Gent-McWilliams eddy parameterization,  sub-grid scale processes. The Gent-McWilliams eddy parameterization,
15  KPP mixing scheme and PV flux parameterization are all dealt with in  KPP mixing scheme and PV flux parameterization are all dealt with in
# Line 14  part of the tracer equations and the var Line 18  part of the tracer equations and the var
18  described here.  described here.
19    
20  \subsection{Time-stepping of tracers: ABII}  \subsection{Time-stepping of tracers: ABII}
21    \label{sect:tracer_equations_abII}
22    \begin{rawhtml}
23    <!-- CMIREDIR:tracer_equations_abII: -->
24    \end{rawhtml}
25    
26  The default advection scheme is the centered second order method which  The default advection scheme is the centered second order method which
27  requires a second order or quasi-second order time-stepping scheme to  requires a second order or quasi-second order time-stepping scheme to
# Line 41  only affects the surface layer since the Line 49  only affects the surface layer since the
49  everywhere else. This term is therefore referred to as the surface  everywhere else. This term is therefore referred to as the surface
50  correction term. Global conservation is not possible using the  correction term. Global conservation is not possible using the
51  flux-form (as here) and a linearized free-surface  flux-form (as here) and a linearized free-surface
52  (\cite{Griffies00,Campin02}).  (\cite{griffies:00,campin:02}).
53    
54  The continuity equation can be recovered by setting  The continuity equation can be recovered by setting
55  $G_{diff}=G_{forc}=0$ and $\tau=1$.  $G_{diff}=G_{forc}=0$ and $\tau=1$.
56    
57  The driver routine that calls the routines to calculate tendancies are  The driver routine that calls the routines to calculate tendencies are
58  {\em S/R CALC\_GT} and {\em S/R CALC\_GS} for temperature and salt  {\em S/R CALC\_GT} and {\em S/R CALC\_GS} for temperature and salt
59  (moisture), respectively. These in turn call a generic advection  (moisture), respectively. These in turn call a generic advection
60  diffusion routine {\em S/R GAD\_CALC\_RHS} that is called with the  diffusion routine {\em S/R GAD\_CALC\_RHS} that is called with the
61  flow field and relevent tracer as arguments and returns the collective  flow field and relevant tracer as arguments and returns the collective
62  tendancy due to advection and diffusion. Forcing is add subsequently  tendency due to advection and diffusion. Forcing is add subsequently
63  in {\em S/R CALC\_GT} or {\em S/R CALC\_GS} to the same tendancy  in {\em S/R CALC\_GT} or {\em S/R CALC\_GS} to the same tendency
64  array.  array.
65    
66  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
# Line 66  $F_r$: {\bf fVerT} (argument) Line 74  $F_r$: {\bf fVerT} (argument)
74    
75  \end{minipage} }  \end{minipage} }
76    
77    The space and time discretization are treated separately (method of
78  The space and time discretizations are treated seperately (method of  lines). Tendencies are calculated at time levels $n$ and $n-1$ and
79  lines). The Adams-Bashforth time discretization reads:  extrapolated to $n+1/2$ using the Adams-Bashforth method:
80  \marginpar{$\epsilon$: {\bf AB\_eps}}  \marginpar{$\epsilon$: {\bf AB\_eps}}
 \marginpar{$\Delta t$: {\bf deltaTtracer}}  
81  \begin{equation}  \begin{equation}
82  \tau^{(n+1)} = \tau^{(n)} + \Delta t \left(  G^{(n+1/2)} =
83  (\frac{3}{2} + \epsilon) G^{(n)} - (\frac{1}{2} + \epsilon) G^{(n-1)}  (\frac{3}{2} + \epsilon) G^{(n)} - (\frac{1}{2} + \epsilon) G^{(n-1)}
 \right)  
84  \end{equation}  \end{equation}
85  where $G^{(n)} = G_{adv}^\tau + G_{diff}^\tau + G_{src}^\tau$ at time  where $G^{(n)} = G_{adv}^\tau + G_{diff}^\tau + G_{src}^\tau$ at time
86  step $n$.  step $n$. The tendency at $n-1$ is not re-calculated but rather the
87    tendency at $n$ is stored in a global array for later re-use.
88    
89  Strictly speaking the ABII scheme should be applied only to the  \fbox{ \begin{minipage}{4.75in}
90  advection terms. However, this scheme is only used in conjuction with  {\em S/R ADAMS\_BASHFORTH2} ({\em model/src/adams\_bashforth2.F})
91  the standard second, third and fourth order advection  
92  schemes. Selection of any other advection scheme disables  $G^{(n+1/2)}$: {\bf gTracer} (argument on exit)
93  Adams-Bashforth for tracers so that explicit diffusion and forcing use  
94  the forward method.  $G^{(n)}$: {\bf gTracer} (argument on entry)
95    
96    $G^{(n-1)}$: {\bf gTrNm1} (argument)
97    
98    $\epsilon$: {\bf ABeps} (PARAMS.h)
99    
100    \end{minipage} }
101    
102    The tracers are stepped forward in time using the extrapolated tendency:
103    \begin{equation}
104    \tau^{(n+1)} = \tau^{(n)} + \Delta t G^{(n+1/2)}
105    \end{equation}
106    \marginpar{$\Delta t$: {\bf deltaTtracer}}
107    
108  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
109  {\em S/R TIMESTEP\_TRACER} ({\em model/src/timestep\_tracer.F})  {\em S/R TIMESTEP\_TRACER} ({\em model/src/timestep\_tracer.F})
110    
111  $\tau$: {\bf tracer} (argument)  $\tau^{(n+1)}$: {\bf gTracer} (argument on exit)
112    
113  $G^{(n)}$: {\bf gTracer} (argument)  $\tau^{(n)}$: {\bf tracer} (argument on entry)
114    
115  $G^{(n-1)}$: {\bf gTrNm1} (argument)  $G^{(n+1/2)}$: {\bf gTracer} (argument)
116    
117  $\Delta t$: {\bf deltaTtracer} (PARAMS.h)  $\Delta t$: {\bf deltaTtracer} (PARAMS.h)
118    
119  \end{minipage} }  \end{minipage} }
120    
121    Strictly speaking the ABII scheme should be applied only to the
122    advection terms. However, this scheme is only used in conjunction with
123    the standard second, third and fourth order advection
124    schemes. Selection of any other advection scheme disables
125    Adams-Bashforth for tracers so that explicit diffusion and forcing use
126    the forward method.
127    
128    
129    
130    
131    \section{Linear advection schemes}
132    \label{sect:tracer-advection}
133    \begin{rawhtml}
134    <!-- CMIREDIR:linear_advection_schemes: -->
135    \end{rawhtml}
136    
137  \begin{figure}  \begin{figure}
138  \resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-lo.eps}}  \resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-lo.eps}}
139  \caption{  \caption{
# Line 136  $\mu=c/(1-c)$. Line 171  $\mu=c/(1-c)$.
171  }  }
172  \end{figure}  \end{figure}
173    
 \section{Linear advection schemes}  
   
174  The advection schemes known as centered second order, centered fourth  The advection schemes known as centered second order, centered fourth
175  order, first order upwind and upwind biased third order are known as  order, first order upwind and upwind biased third order are known as
176  linear advection schemes because the coefficient for interpolation of  linear advection schemes because the coefficient for interpolation of
# Line 148  commonly used in the field and most fami Line 181  commonly used in the field and most fami
181  \subsection{Centered second order advection-diffusion}  \subsection{Centered second order advection-diffusion}
182    
183  The basic discretization, centered second order, is the default. It is  The basic discretization, centered second order, is the default. It is
184  designed to be consistant with the continuity equation to facilitate  designed to be consistent with the continuity equation to facilitate
185  conservation properties analogous to the continuum. However, centered  conservation properties analogous to the continuum. However, centered
186  second order advection is notoriously noisey and must be used in  second order advection is notoriously noisy and must be used in
187  conjuction with some finite amount of diffusion to produce a sensible  conjunction with some finite amount of diffusion to produce a sensible
188  solution.  solution.
189    
190  The advection operator is discretized:  The advection operator is discretized:
# Line 177  W & = & {\cal A}_c w Line 210  W & = & {\cal A}_c w
210    
211  For non-divergent flow, this discretization can be shown to conserve  For non-divergent flow, this discretization can be shown to conserve
212  the tracer both locally and globally and to globally conserve tracer  the tracer both locally and globally and to globally conserve tracer
213  variance, $\tau^2$. The proof is given in \cite{Adcroft95,Adcroft97}.  variance, $\tau^2$. The proof is given in \cite{adcroft:95,adcroft:97}.
214    
215  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
216  {\em S/R GAD\_C2\_ADV\_X} ({\em gad\_c2\_adv\_x.F})  {\em S/R GAD\_C2\_ADV\_X} ({\em gad\_c2\_adv\_x.F})
# Line 225  F_r & = & W \overline{\tau - \frac{1}{6} Line 258  F_r & = & W \overline{\tau - \frac{1}{6}
258  \end{eqnarray}  \end{eqnarray}
259    
260  At boundaries, $\delta_{\hat{n}} \tau$ is set to zero allowing  At boundaries, $\delta_{\hat{n}} \tau$ is set to zero allowing
261  $\delta_{nn}$ to be evaluated. We are currently examing the accuracy  $\delta_{nn}$ to be evaluated. We are currently examine the accuracy
262  of this boundary condition and the effect on the solution.  of this boundary condition and the effect on the solution.
263    
264    \fbox{ \begin{minipage}{4.75in}
265    {\em S/R GAD\_U3\_ADV\_X} ({\em gad\_u3\_adv\_x.F})
266    
267    $F_x$: {\bf uT} (argument)
268    
269    $U$: {\bf uTrans} (argument)
270    
271    $\tau$: {\bf tracer} (argument)
272    
273    {\em S/R GAD\_U3\_ADV\_Y} ({\em gad\_u3\_adv\_y.F})
274    
275    $F_y$: {\bf vT} (argument)
276    
277    $V$: {\bf vTrans} (argument)
278    
279    $\tau$: {\bf tracer} (argument)
280    
281    {\em S/R GAD\_U3\_ADV\_R} ({\em gad\_u3\_adv\_r.F})
282    
283    $F_r$: {\bf wT} (argument)
284    
285    $W$: {\bf rTrans} (argument)
286    
287    $\tau$: {\bf tracer} (argument)
288    
289    \end{minipage} }
290    
291  \subsection{Centered fourth order advection}  \subsection{Centered fourth order advection}
292    
293  Centered fourth order advection is formally the most accurate scheme  Centered fourth order advection is formally the most accurate scheme
294  we have implemented and can be used to great effect in high resolution  we have implemented and can be used to great effect in high resolution
295  simultation where dynamical scales are well resolved. However, the  simulation where dynamical scales are well resolved. However, the
296  scheme is noisey like the centered second order method and so must be  scheme is noisy like the centered second order method and so must be
297  used with some finite amount of diffusion. Bi-harmonic is recommended  used with some finite amount of diffusion. Bi-harmonic is recommended
298  since it is more scale selective and less likely to diffuse away the  since it is more scale selective and less likely to diffuse away the
299  well resolved gradient the fourth order scheme worked so hard to  well resolved gradient the fourth order scheme worked so hard to
# Line 248  F_r & = & W \overline{\tau - \frac{1}{6} Line 307  F_r & = & W \overline{\tau - \frac{1}{6}
307  \end{eqnarray}  \end{eqnarray}
308    
309  As for the third order scheme, the best discretization near boundaries  As for the third order scheme, the best discretization near boundaries
310  is under investigation but currenlty $\delta_i \tau=0$ on a boundary.  is under investigation but currently $\delta_i \tau=0$ on a boundary.
311    
312    \fbox{ \begin{minipage}{4.75in}
313    {\em S/R GAD\_C4\_ADV\_X} ({\em gad\_c4\_adv\_x.F})
314    
315    $F_x$: {\bf uT} (argument)
316    
317    $U$: {\bf uTrans} (argument)
318    
319    $\tau$: {\bf tracer} (argument)
320    
321    {\em S/R GAD\_C4\_ADV\_Y} ({\em gad\_c4\_adv\_y.F})
322    
323    $F_y$: {\bf vT} (argument)
324    
325    $V$: {\bf vTrans} (argument)
326    
327    $\tau$: {\bf tracer} (argument)
328    
329    {\em S/R GAD\_C4\_ADV\_R} ({\em gad\_c4\_adv\_r.F})
330    
331    $F_r$: {\bf wT} (argument)
332    
333    $W$: {\bf rTrans} (argument)
334    
335    $\tau$: {\bf tracer} (argument)
336    
337    \end{minipage} }
338    
339    
340  \subsection{First order upwind advection}  \subsection{First order upwind advection}
341    
# Line 272  if the limiter is set to zero. Line 359  if the limiter is set to zero.
359    
360    
361  \section{Non-linear advection schemes}  \section{Non-linear advection schemes}
362    \begin{rawhtml}
363    <!-- CMIREDIR:non-linear_advection_schemes: -->
364    \end{rawhtml}
365    
366  Non-linear advection schemes invoke non-linear interpolation and are  Non-linear advection schemes invoke non-linear interpolation and are
367  widely used in computational fluid dynamics (non-linear does not refer  widely used in computational fluid dynamics (non-linear does not refer
# Line 310  r = \frac{ \tau_{i-1} - \tau_{i-2} }{ \t Line 400  r = \frac{ \tau_{i-1} - \tau_{i-2} }{ \t
400  r = \frac{ \tau_{i+1} - \tau_{i} }{ \tau_{i} - \tau_{i-1} } & \forall & u < 0  r = \frac{ \tau_{i+1} - \tau_{i} }{ \tau_{i} - \tau_{i-1} } & \forall & u < 0
401  \end{eqnarray}  \end{eqnarray}
402  as it's argument. There are many choices of limiter function but we  as it's argument. There are many choices of limiter function but we
403  only provide the Superbee limiter \cite{Roe85}:  only provide the Superbee limiter \cite{roe:85}:
404  \begin{equation}  \begin{equation}
405  \psi(r) = \max[0,\min[1,2r],\min[2,r]]  \psi(r) = \max[0,\min[1,2r],\min[2,r]]
406  \end{equation}  \end{equation}
407    
408    \fbox{ \begin{minipage}{4.75in}
409    {\em S/R GAD\_FLUXLIMIT\_ADV\_X} ({\em gad\_fluxlimit\_adv\_x.F})
410    
411    $F_x$: {\bf uT} (argument)
412    
413    $U$: {\bf uTrans} (argument)
414    
415    $\tau$: {\bf tracer} (argument)
416    
417    {\em S/R GAD\_FLUXLIMIT\_ADV\_Y} ({\em gad\_fluxlimit\_adv\_y.F})
418    
419    $F_y$: {\bf vT} (argument)
420    
421    $V$: {\bf vTrans} (argument)
422    
423    $\tau$: {\bf tracer} (argument)
424    
425    {\em S/R GAD\_FLUXLIMIT\_ADV\_R} ({\em gad\_fluxlimit\_adv\_r.F})
426    
427    $F_r$: {\bf wT} (argument)
428    
429    $W$: {\bf rTrans} (argument)
430    
431    $\tau$: {\bf tracer} (argument)
432    
433    \end{minipage} }
434    
435    
436  \subsection{Third order direct space time}  \subsection{Third order direct space time}
437    
438  The direct-space-time method deals with space and time discretization  The direct-space-time method deals with space and time discretization
439  together (other methods that treat space and time seperately are known  together (other methods that treat space and time separately are known
440  collectively as the ``Method of Lines''). The Lax-Wendroff scheme  collectively as the ``Method of Lines''). The Lax-Wendroff scheme
441  falls into this category; it adds sufficient diffusion to a second  falls into this category; it adds sufficient diffusion to a second
442  order flux that the forward-in-time method is stable. The upwind  order flux that the forward-in-time method is stable. The upwind
# Line 337  where Line 454  where
454  d_1 & = & \frac{1}{6} ( 2 - |c| ) ( 1 - |c| ) \\  d_1 & = & \frac{1}{6} ( 2 - |c| ) ( 1 - |c| ) \\
455  d_2 & = & \frac{1}{6} ( 1 - |c| ) ( 1 + |c| )  d_2 & = & \frac{1}{6} ( 1 - |c| ) ( 1 + |c| )
456  \end{eqnarray}  \end{eqnarray}
457  The coefficients $d_0$ and $d_1$ approach $1/3$ and $1/6$ repectively  The coefficients $d_0$ and $d_1$ approach $1/3$ and $1/6$ respectively
458  as the Courant number, $c$, vanishes. In this limit, the conventional  as the Courant number, $c$, vanishes. In this limit, the conventional
459  third order upwind method is recovered. For finite Courant number, the  third order upwind method is recovered. For finite Courant number, the
460  deviations from the linear method are analogous to the diffusion added  deviations from the linear method are analogous to the diffusion added
# Line 345  to centered second order advection in th Line 462  to centered second order advection in th
462    
463  The DST3 method described above must be used in a forward-in-time  The DST3 method described above must be used in a forward-in-time
464  manner and is stable for $0 \le |c| \le 1$. Although the scheme  manner and is stable for $0 \le |c| \le 1$. Although the scheme
465  appears to be forward-in-time, it is in fact second order in time and  appears to be forward-in-time, it is in fact third order in time and
466  the accuracy increases with the Courant number! For low Courant  the accuracy increases with the Courant number! For low Courant
467  number, DST3 produces very similar results (indistinguishable in  number, DST3 produces very similar results (indistinguishable in
468  Fig.~\ref{fig:advect-1d-lo}) to the linear third order method but for  Fig.~\ref{fig:advect-1d-lo}) to the linear third order method but for
# Line 353  large Courant number, where the linear u Line 470  large Courant number, where the linear u
470  unstable, the scheme is extremely accurate  unstable, the scheme is extremely accurate
471  (Fig.~\ref{fig:advect-1d-hi}) with only minor overshoots.  (Fig.~\ref{fig:advect-1d-hi}) with only minor overshoots.
472    
473    \fbox{ \begin{minipage}{4.75in}
474    {\em S/R GAD\_DST3\_ADV\_X} ({\em gad\_dst3\_adv\_x.F})
475    
476    $F_x$: {\bf uT} (argument)
477    
478    $U$: {\bf uTrans} (argument)
479    
480    $\tau$: {\bf tracer} (argument)
481    
482    {\em S/R GAD\_DST3\_ADV\_Y} ({\em gad\_dst3\_adv\_y.F})
483    
484    $F_y$: {\bf vT} (argument)
485    
486    $V$: {\bf vTrans} (argument)
487    
488    $\tau$: {\bf tracer} (argument)
489    
490    {\em S/R GAD\_DST3\_ADV\_R} ({\em gad\_dst3\_adv\_r.F})
491    
492    $F_r$: {\bf wT} (argument)
493    
494    $W$: {\bf rTrans} (argument)
495    
496    $\tau$: {\bf tracer} (argument)
497    
498    \end{minipage} }
499    
500    
501  \subsection{Third order direct space time with flux limiting}  \subsection{Third order direct space time with flux limiting}
502    
503  The overshoots in the DST3 method can be controlled with a flux limiter.  The overshoots in the DST3 method can be controlled with a flux limiter.
# Line 373  and the limiter is the Sweby limiter: Line 518  and the limiter is the Sweby limiter:
518  \psi(r) = \max[0, \min[\min(1,d_0+d_1r],\frac{1-c}{c}r ]]  \psi(r) = \max[0, \min[\min(1,d_0+d_1r],\frac{1-c}{c}r ]]
519  \end{equation}  \end{equation}
520    
521    \fbox{ \begin{minipage}{4.75in}
522    {\em S/R GAD\_DST3FL\_ADV\_X} ({\em gad\_dst3\_adv\_x.F})
523    
524    $F_x$: {\bf uT} (argument)
525    
526    $U$: {\bf uTrans} (argument)
527    
528    $\tau$: {\bf tracer} (argument)
529    
530    {\em S/R GAD\_DST3FL\_ADV\_Y} ({\em gad\_dst3\_adv\_y.F})
531    
532    $F_y$: {\bf vT} (argument)
533    
534    $V$: {\bf vTrans} (argument)
535    
536    $\tau$: {\bf tracer} (argument)
537    
538    {\em S/R GAD\_DST3FL\_ADV\_R} ({\em gad\_dst3\_adv\_r.F})
539    
540    $F_r$: {\bf wT} (argument)
541    
542    $W$: {\bf rTrans} (argument)
543    
544    $\tau$: {\bf tracer} (argument)
545    
546    \end{minipage} }
547    
548    
549  \subsection{Multi-dimensional advection}  \subsection{Multi-dimensional advection}
550    
551  In many of the aforementioned advection schemes the behaviour in  \begin{figure}
552    \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-lo-diag.eps}}
553    \caption{
554    Comparison of advection schemes in two dimensions; diagonal advection
555    of a resolved Gaussian feature. Courant number is 0.01 with
556    30$\times$30 points and solutions are shown for T=1/2. White lines
557    indicate zero crossing (ie. the presence of false minima).  The left
558    column shows the second order schemes; top) centered second order with
559    Adams-Bashforth, middle) Lax-Wendroff and bottom) Superbee flux
560    limited. The middle column shows the third order schemes; top) upwind
561    biased third order with Adams-Bashforth, middle) third order direct
562    space-time method and bottom) the same with flux limiting. The top
563    right panel shows the centered fourth order scheme with
564    Adams-Bashforth and right middle panel shows a fourth order variant on
565    the DST method. Bottom right panel shows the Superbee flux limiter
566    (second order) applied independently in each direction (method of
567    lines).
568    \label{fig:advect-2d-lo-diag}
569    }
570    \end{figure}
571    
572    \begin{figure}
573    \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-mid-diag.eps}}
574    \caption{
575    Comparison of advection schemes in two dimensions; diagonal advection
576    of a resolved Gaussian feature. Courant number is 0.27 with
577    30$\times$30 points and solutions are shown for T=1/2. White lines
578    indicate zero crossing (ie. the presence of false minima).  The left
579    column shows the second order schemes; top) centered second order with
580    Adams-Bashforth, middle) Lax-Wendroff and bottom) Superbee flux
581    limited. The middle column shows the third order schemes; top) upwind
582    biased third order with Adams-Bashforth, middle) third order direct
583    space-time method and bottom) the same with flux limiting. The top
584    right panel shows the centered fourth order scheme with
585    Adams-Bashforth and right middle panel shows a fourth order variant on
586    the DST method. Bottom right panel shows the Superbee flux limiter
587    (second order) applied independently in each direction (method of
588    lines).
589    \label{fig:advect-2d-mid-diag}
590    }
591    \end{figure}
592    
593    \begin{figure}
594    \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-hi-diag.eps}}
595    \caption{
596    Comparison of advection schemes in two dimensions; diagonal advection
597    of a resolved Gaussian feature. Courant number is 0.47 with
598    30$\times$30 points and solutions are shown for T=1/2. White lines
599    indicate zero crossings and initial maximum values (ie. the presence
600    of false extrema).  The left column shows the second order schemes;
601    top) centered second order with Adams-Bashforth, middle) Lax-Wendroff
602    and bottom) Superbee flux limited. The middle column shows the third
603    order schemes; top) upwind biased third order with Adams-Bashforth,
604    middle) third order direct space-time method and bottom) the same with
605    flux limiting. The top right panel shows the centered fourth order
606    scheme with Adams-Bashforth and right middle panel shows a fourth
607    order variant on the DST method. Bottom right panel shows the Superbee
608    flux limiter (second order) applied independently in each direction
609    (method of lines).
610    \label{fig:advect-2d-hi-diag}
611    }
612    \end{figure}
613    
614    
615    
616    In many of the aforementioned advection schemes the behavior in
617  multiple dimensions is not necessarily as good as the one dimensional  multiple dimensions is not necessarily as good as the one dimensional
618  behaviour. For instance, a shape preserving monotonic scheme in one  behavior. For instance, a shape preserving monotonic scheme in one
619  dimension can have severe shape distortion in two dimensions if the  dimension can have severe shape distortion in two dimensions if the
620  two components of horizontal fluxes are treated independently. There  two components of horizontal fluxes are treated independently. There
621  is a large body of literature on the subject dealing with this problem  is a large body of literature on the subject dealing with this problem
# Line 389  as if in one dimension: Line 627  as if in one dimension:
627  \tau^{n+1/3} & = & \tau^{n}  \tau^{n+1/3} & = & \tau^{n}
628  - \Delta t \left( \frac{1}{\Delta x} \delta_i F^x(\tau^{n})  - \Delta t \left( \frac{1}{\Delta x} \delta_i F^x(\tau^{n})
629             + \tau^{n} \frac{1}{\Delta x} \delta_i u \right) \\             + \tau^{n} \frac{1}{\Delta x} \delta_i u \right) \\
630  \tau^{n+2/3} & = & \tau^{n}  \tau^{n+2/3} & = & \tau^{n+1/3}
631  - \Delta t \left( \frac{1}{\Delta y} \delta_j F^y(\tau^{n+1/3})  - \Delta t \left( \frac{1}{\Delta y} \delta_j F^y(\tau^{n+1/3})
632             + \tau^{n} \frac{1}{\Delta y} \delta_i v \right) \\             + \tau^{n} \frac{1}{\Delta y} \delta_i v \right) \\
633  \tau^{n+3/3} & = & \tau^{n}  \tau^{n+3/3} & = & \tau^{n+2/3}
634  - \Delta t \left( \frac{1}{\Delta r} \delta_k F^x(\tau^{n+2/3})  - \Delta t \left( \frac{1}{\Delta r} \delta_k F^x(\tau^{n+2/3})
635             + \tau^{n} \frac{1}{\Delta r} \delta_i w \right)             + \tau^{n} \frac{1}{\Delta r} \delta_i w \right)
636  \end{eqnarray}  \end{eqnarray}
637    
638  In order to incorporate this method into the general model algorithm,  In order to incorporate this method into the general model algorithm,
639  we compute the effective tendancy rather than update the tracer so  we compute the effective tendency rather than update the tracer so
640  that other terms such as diffusion are using the $n$ time-level and  that other terms such as diffusion are using the $n$ time-level and
641  not the updated $n+3/3$ quantities:  not the updated $n+3/3$ quantities:
642  \begin{equation}  \begin{equation}
# Line 408  So that the over all time-stepping looks Line 646  So that the over all time-stepping looks
646  \begin{equation}  \begin{equation}
647  \tau^{n+1} = \tau^{n} + \Delta t \left( G^{n+1/2}_{adv} + G_{diff}(\tau^{n}) + G^{n}_{forcing} \right)  \tau^{n+1} = \tau^{n} + \Delta t \left( G^{n+1/2}_{adv} + G_{diff}(\tau^{n}) + G^{n}_{forcing} \right)
648  \end{equation}  \end{equation}
649    
650    \fbox{ \begin{minipage}{4.75in}
651    {\em S/R GAD\_ADVECTION} ({\em gad\_advection.F})
652    
653    $\tau$: {\bf Tracer} (argument)
654    
655    $G^{n+1/2}_{adv}$: {\bf Gtracer} (argument)
656    
657    $F_x, F_y, F_r$: {\bf af} (local)
658    
659    $U$: {\bf uTrans} (local)
660    
661    $V$: {\bf vTrans} (local)
662    
663    $W$: {\bf rTrans} (local)
664    
665    \end{minipage} }
666    
667    
668    \section{Comparison of advection schemes}
669    \label{sect:tracer_advection_schemes}
670    \begin{rawhtml}
671    <!-- CMIREDIR:comparison_of_advection_schemes: -->
672    \end{rawhtml}
673    
674    \begin{table}[htb]
675    \centering
676     \begin{tabular}[htb]{|l|c|c|c|c|l|}
677       \hline
678       Advection Scheme & code & use  & use Multi- & Stencil & comments \\
679                        &      & A.B. & dimension & (1 dim) & \\
680       \hline \hline
681       $1^{rst}$order upwind  & 1 &  No & Yes & 3 pts & linear/$\tau$, non-linear/v\\
682       \hline
683       centered $2^{nd}$order & 2 &  Yes & No & 3 pts & linear \\
684       \hline
685       $3^{rd}$order upwind   & 3 &  Yes & No & 5 pts & linear/$\tau$\\
686       \hline
687       centered $4^{th}$order & 4 &  Yes & No & 5 pts & linear \\
688       \hline \hline
689       $2^{nd}$order DST (Lax-Wendroff)  & 20 &
690                             No & Yes & 3 pts & linear/$\tau$, non-linear/v\\
691       \hline
692       $3^{rd}$order DST & 30 &  No & Yes & 5 pts & linear/$\tau$, non-linear/v\\
693       \hline \hline
694       $2^{nd}$order Flux Limiters & 77 &  No & Yes & 5 pts & non-linear \\
695       \hline
696       $3^{nd}$order DST Flux limiter & 33 &  No & Yes & 5 pts & non-linear \\
697       \hline
698     \end{tabular}
699     \caption{Summary of the different advection schemes available in MITgcm.
700              ``A.B.'' stands for Adams-Bashforth and ``DST'' for direct space time.
701              The code corresponds to the number used to select the corresponding
702              advection scheme in the parameter file (e.g., {\bf tempAdvScheme}=3 in
703              file {\em data} selects the $3^{rd}$ order upwind advection scheme
704              for temperature).
705       }
706     \label{tab:advectionShemes_summary}
707    \end{table}
708    
709    
710    Figs.~\ref{fig:advect-2d-lo-diag}, \ref{fig:advect-2d-mid-diag} and
711    \ref{fig:advect-2d-hi-diag} show solutions to a simple diagonal
712    advection problem using a selection of schemes for low, moderate and
713    high Courant numbers, respectively. The top row shows the linear
714    schemes, integrated with the Adams-Bashforth method. Theses schemes
715    are clearly unstable for the high Courant number and weakly unstable
716    for the moderate Courant number. The presence of false extrema is very
717    apparent for all Courant numbers. The middle row shows solutions
718    obtained with the unlimited but multi-dimensional schemes. These
719    solutions also exhibit false extrema though the pattern now shows
720    symmetry due to the multi-dimensional scheme. Also, the schemes are
721    stable at high Courant number where the linear schemes weren't. The
722    bottom row (left and middle) shows the limited schemes and most
723    obvious is the absence of false extrema. The accuracy and stability of
724    the unlimited non-linear schemes is retained at high Courant number
725    but at low Courant number the tendency is to loose amplitude in sharp
726    peaks due to diffusion. The one dimensional tests shown in
727    Figs.~\ref{fig:advect-1d-lo} and \ref{fig:advect-1d-hi} showed this
728    phenomenon.
729    
730    Finally, the bottom left and right panels use the same advection
731    scheme but the right does not use the multi-dimensional method. At low
732    Courant number this appears to not matter but for moderate Courant
733    number severe distortion of the feature is apparent. Moreover, the
734    stability of the multi-dimensional scheme is determined by the maximum
735    Courant number applied of each dimension while the stability of the
736    method of lines is determined by the sum. Hence, in the high Courant
737    number plot, the scheme is unstable.
738    
739    With many advection schemes implemented in the code two questions
740    arise: ``Which scheme is best?'' and ``Why don't you just offer the
741    best advection scheme?''. Unfortunately, no one advection scheme is
742    ``the best'' for all particular applications and for new applications
743    it is often a matter of trial to determine which is most
744    suitable. Here are some guidelines but these are not the rule;
745    \begin{itemize}
746    \item If you have a coarsely resolved model, using a
747    positive or upwind biased scheme will introduce significant diffusion
748    to the solution and using a centered higher order scheme will
749    introduce more noise. In this case, simplest may be best.
750    \item If you have a high resolution model, using a higher order
751    scheme will give a more accurate solution but scale-selective
752    diffusion might need to be employed. The flux limited methods
753    offer similar accuracy in this regime.
754    \item If your solution has shocks or propagating fronts then a
755    flux limited scheme is almost essential.
756    \item If your time-step is limited by advection, the multi-dimensional
757    non-linear schemes have the most stability (up to Courant number 1).
758    \item If you need to know how much diffusion/dissipation has occurred you
759    will have a lot of trouble figuring it out with a non-linear method.
760    \item The presence of false extrema is non-physical and this alone is the
761    strongest argument for using a positive scheme.
762    \end{itemize}

Legend:
Removed from v.1.3  
changed lines
  Added in v.1.22

  ViewVC Help
Powered by ViewVC 1.1.22