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

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.24

  ViewVC Help
Powered by ViewVC 1.1.22