/[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.5 by cnh, Wed Oct 24 14:19:34 2001 UTC
# Line 2  Line 2 
2  % $Name$  % $Name$
3    
4  \section{Tracer equations}  \section{Tracer equations}
5    \label{sec:tracer_eqautions}
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
# Line 66  $F_r$: {\bf fVerT} (argument) Line 67  $F_r$: {\bf fVerT} (argument)
67    
68  \end{minipage} }  \end{minipage} }
69    
   
70  The space and time discretizations are treated seperately (method of  The space and time discretizations are treated seperately (method of
71  lines). The Adams-Bashforth time discretization reads:  lines). Tendancies are calculated at time levels $n$ and $n-1$ and
72    extrapolated to $n+1/2$ using the Adams-Bashforth method:
73  \marginpar{$\epsilon$: {\bf AB\_eps}}  \marginpar{$\epsilon$: {\bf AB\_eps}}
 \marginpar{$\Delta t$: {\bf deltaTtracer}}  
74  \begin{equation}  \begin{equation}
75  \tau^{(n+1)} = \tau^{(n)} + \Delta t \left(  G^{(n+1/2)} =
76  (\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)  
77  \end{equation}  \end{equation}
78  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
79  step $n$.  step $n$. The tendancy at $n-1$ is not re-calculated but rather the
80    tendancy at $n$ is stored in a global array for later re-use.
81    
82  Strictly speaking the ABII scheme should be applied only to the  \fbox{ \begin{minipage}{4.75in}
83  advection terms. However, this scheme is only used in conjuction with  {\em S/R ADAMS\_BASHFORTH2} ({\em model/src/adams\_bashforth2.F})
84  the standard second, third and fourth order advection  
85  schemes. Selection of any other advection scheme disables  $G^{(n+1/2)}$: {\bf gTracer} (argument on exit)
86  Adams-Bashforth for tracers so that explicit diffusion and forcing use  
87  the forward method.  $G^{(n)}$: {\bf gTracer} (argument on entry)
88    
89    $G^{(n-1)}$: {\bf gTrNm1} (argument)
90    
91    $\epsilon$: {\bf ABeps} (PARAMS.h)
92    
93    \end{minipage} }
94    
95    The tracers are stepped forward in time using the extrapolated tendancy:
96    \begin{equation}
97    \tau^{(n+1)} = \tau^{(n)} + \Delta t G^{(n+1/2)}
98    \end{equation}
99    \marginpar{$\Delta t$: {\bf deltaTtracer}}
100    
101  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
102  {\em S/R TIMESTEP\_TRACER} ({\em model/src/timestep\_tracer.F})  {\em S/R TIMESTEP\_TRACER} ({\em model/src/timestep\_tracer.F})
103    
104  $\tau$: {\bf tracer} (argument)  $\tau^{(n+1)}$: {\bf gTracer} (argument on exit)
105    
106  $G^{(n)}$: {\bf gTracer} (argument)  $\tau^{(n)}$: {\bf tracer} (argument on entry)
107    
108  $G^{(n-1)}$: {\bf gTrNm1} (argument)  $G^{(n+1/2)}$: {\bf gTracer} (argument)
109    
110  $\Delta t$: {\bf deltaTtracer} (PARAMS.h)  $\Delta t$: {\bf deltaTtracer} (PARAMS.h)
111    
112  \end{minipage} }  \end{minipage} }
113    
114    Strictly speaking the ABII scheme should be applied only to the
115    advection terms. However, this scheme is only used in conjuction with
116    the standard second, third and fourth order advection
117    schemes. Selection of any other advection scheme disables
118    Adams-Bashforth for tracers so that explicit diffusion and forcing use
119    the forward method.
120    
121    
122    
123    
124    \section{Linear advection schemes}
125    
126  \begin{figure}  \begin{figure}
127  \resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-lo.eps}}  \resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-lo.eps}}
128  \caption{  \caption{
# Line 136  $\mu=c/(1-c)$. Line 160  $\mu=c/(1-c)$.
160  }  }
161  \end{figure}  \end{figure}
162    
 \section{Linear advection schemes}  
   
163  The advection schemes known as centered second order, centered fourth  The advection schemes known as centered second order, centered fourth
164  order, first order upwind and upwind biased third order are known as  order, first order upwind and upwind biased third order are known as
165  linear advection schemes because the coefficient for interpolation of  linear advection schemes because the coefficient for interpolation of
# Line 228  At boundaries, $\delta_{\hat{n}} \tau$ i Line 250  At boundaries, $\delta_{\hat{n}} \tau$ i
250  $\delta_{nn}$ to be evaluated. We are currently examing the accuracy  $\delta_{nn}$ to be evaluated. We are currently examing the accuracy
251  of this boundary condition and the effect on the solution.  of this boundary condition and the effect on the solution.
252    
253    \fbox{ \begin{minipage}{4.75in}
254    {\em S/R GAD\_U3\_ADV\_X} ({\em gad\_u3\_adv\_x.F})
255    
256    $F_x$: {\bf uT} (argument)
257    
258    $U$: {\bf uTrans} (argument)
259    
260    $\tau$: {\bf tracer} (argument)
261    
262    {\em S/R GAD\_U3\_ADV\_Y} ({\em gad\_u3\_adv\_y.F})
263    
264    $F_y$: {\bf vT} (argument)
265    
266    $V$: {\bf vTrans} (argument)
267    
268    $\tau$: {\bf tracer} (argument)
269    
270    {\em S/R GAD\_U3\_ADV\_R} ({\em gad\_u3\_adv\_r.F})
271    
272    $F_r$: {\bf wT} (argument)
273    
274    $W$: {\bf rTrans} (argument)
275    
276    $\tau$: {\bf tracer} (argument)
277    
278    \end{minipage} }
279    
280  \subsection{Centered fourth order advection}  \subsection{Centered fourth order advection}
281    
# Line 250  F_r & = & W \overline{\tau - \frac{1}{6} Line 298  F_r & = & W \overline{\tau - \frac{1}{6}
298  As for the third order scheme, the best discretization near boundaries  As for the third order scheme, the best discretization near boundaries
299  is under investigation but currenlty $\delta_i \tau=0$ on a boundary.  is under investigation but currenlty $\delta_i \tau=0$ on a boundary.
300    
301    \fbox{ \begin{minipage}{4.75in}
302    {\em S/R GAD\_C4\_ADV\_X} ({\em gad\_c4\_adv\_x.F})
303    
304    $F_x$: {\bf uT} (argument)
305    
306    $U$: {\bf uTrans} (argument)
307    
308    $\tau$: {\bf tracer} (argument)
309    
310    {\em S/R GAD\_C4\_ADV\_Y} ({\em gad\_c4\_adv\_y.F})
311    
312    $F_y$: {\bf vT} (argument)
313    
314    $V$: {\bf vTrans} (argument)
315    
316    $\tau$: {\bf tracer} (argument)
317    
318    {\em S/R GAD\_C4\_ADV\_R} ({\em gad\_c4\_adv\_r.F})
319    
320    $F_r$: {\bf wT} (argument)
321    
322    $W$: {\bf rTrans} (argument)
323    
324    $\tau$: {\bf tracer} (argument)
325    
326    \end{minipage} }
327    
328    
329  \subsection{First order upwind advection}  \subsection{First order upwind advection}
330    
331  Although the upwind scheme is the underlying scheme for the robust or  Although the upwind scheme is the underlying scheme for the robust or
# Line 315  only provide the Superbee limiter \cite{ Line 391  only provide the Superbee limiter \cite{
391  \psi(r) = \max[0,\min[1,2r],\min[2,r]]  \psi(r) = \max[0,\min[1,2r],\min[2,r]]
392  \end{equation}  \end{equation}
393    
394    \fbox{ \begin{minipage}{4.75in}
395    {\em S/R GAD\_FLUXLIMIT\_ADV\_X} ({\em gad\_fluxlimit\_adv\_x.F})
396    
397    $F_x$: {\bf uT} (argument)
398    
399    $U$: {\bf uTrans} (argument)
400    
401    $\tau$: {\bf tracer} (argument)
402    
403    {\em S/R GAD\_FLUXLIMIT\_ADV\_Y} ({\em gad\_fluxlimit\_adv\_y.F})
404    
405    $F_y$: {\bf vT} (argument)
406    
407    $V$: {\bf vTrans} (argument)
408    
409    $\tau$: {\bf tracer} (argument)
410    
411    {\em S/R GAD\_FLUXLIMIT\_ADV\_R} ({\em gad\_fluxlimit\_adv\_r.F})
412    
413    $F_r$: {\bf wT} (argument)
414    
415    $W$: {\bf rTrans} (argument)
416    
417    $\tau$: {\bf tracer} (argument)
418    
419    \end{minipage} }
420    
421    
422  \subsection{Third order direct space time}  \subsection{Third order direct space time}
423    
# Line 353  large Courant number, where the linear u Line 456  large Courant number, where the linear u
456  unstable, the scheme is extremely accurate  unstable, the scheme is extremely accurate
457  (Fig.~\ref{fig:advect-1d-hi}) with only minor overshoots.  (Fig.~\ref{fig:advect-1d-hi}) with only minor overshoots.
458    
459    \fbox{ \begin{minipage}{4.75in}
460    {\em S/R GAD\_DST3\_ADV\_X} ({\em gad\_dst3\_adv\_x.F})
461    
462    $F_x$: {\bf uT} (argument)
463    
464    $U$: {\bf uTrans} (argument)
465    
466    $\tau$: {\bf tracer} (argument)
467    
468    {\em S/R GAD\_DST3\_ADV\_Y} ({\em gad\_dst3\_adv\_y.F})
469    
470    $F_y$: {\bf vT} (argument)
471    
472    $V$: {\bf vTrans} (argument)
473    
474    $\tau$: {\bf tracer} (argument)
475    
476    {\em S/R GAD\_DST3\_ADV\_R} ({\em gad\_dst3\_adv\_r.F})
477    
478    $F_r$: {\bf wT} (argument)
479    
480    $W$: {\bf rTrans} (argument)
481    
482    $\tau$: {\bf tracer} (argument)
483    
484    \end{minipage} }
485    
486    
487  \subsection{Third order direct space time with flux limiting}  \subsection{Third order direct space time with flux limiting}
488    
489  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 504  and the limiter is the Sweby limiter:
504  \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 ]]
505  \end{equation}  \end{equation}
506    
507    \fbox{ \begin{minipage}{4.75in}
508    {\em S/R GAD\_DST3FL\_ADV\_X} ({\em gad\_dst3\_adv\_x.F})
509    
510    $F_x$: {\bf uT} (argument)
511    
512    $U$: {\bf uTrans} (argument)
513    
514    $\tau$: {\bf tracer} (argument)
515    
516    {\em S/R GAD\_DST3FL\_ADV\_Y} ({\em gad\_dst3\_adv\_y.F})
517    
518    $F_y$: {\bf vT} (argument)
519    
520    $V$: {\bf vTrans} (argument)
521    
522    $\tau$: {\bf tracer} (argument)
523    
524    {\em S/R GAD\_DST3FL\_ADV\_R} ({\em gad\_dst3\_adv\_r.F})
525    
526    $F_r$: {\bf wT} (argument)
527    
528    $W$: {\bf rTrans} (argument)
529    
530    $\tau$: {\bf tracer} (argument)
531    
532    \end{minipage} }
533    
534    
535  \subsection{Multi-dimensional advection}  \subsection{Multi-dimensional advection}
536    
537    \begin{figure}
538    \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-lo-diag.eps}}
539    \caption{
540    Comparison of advection schemes in two dimensions; diagonal advection
541    of a resolved Guassian feature. Courant number is 0.01 with
542    30$\times$30 points and solutions are shown for T=1/2. White lines
543    indicate zero crossing (ie. the presence of false minima).  The left
544    column shows the second order schemes; top) centered second order with
545    Adams-Bashforth, middle) Lax-Wendroff and bottom) Superbee flux
546    limited. The middle column shows the third order schemes; top) upwind
547    biased third order with Adams-Bashforth, middle) third order direct
548    space-time method and bottom) the same with flux limiting. The top
549    right panel shows the centered fourth order scheme with
550    Adams-Bashforth and right middle panel shows a fourth order variant on
551    the DST method. Bottom right panel shows the Superbee flux limiter
552    (second order) applied independantly in each direction (method of
553    lines).
554    \label{fig:advect-2d-lo-diag}
555    }
556    \end{figure}
557    
558    \begin{figure}
559    \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-mid-diag.eps}}
560    \caption{
561    Comparison of advection schemes in two dimensions; diagonal advection
562    of a resolved Guassian feature. Courant number is 0.27 with
563    30$\times$30 points and solutions are shown for T=1/2. White lines
564    indicate zero crossing (ie. the presence of false minima).  The left
565    column shows the second order schemes; top) centered second order with
566    Adams-Bashforth, middle) Lax-Wendroff and bottom) Superbee flux
567    limited. The middle column shows the third order schemes; top) upwind
568    biased third order with Adams-Bashforth, middle) third order direct
569    space-time method and bottom) the same with flux limiting. The top
570    right panel shows the centered fourth order scheme with
571    Adams-Bashforth and right middle panel shows a fourth order variant on
572    the DST method. Bottom right panel shows the Superbee flux limiter
573    (second order) applied independantly in each direction (method of
574    lines).
575    \label{fig:advect-2d-mid-diag}
576    }
577    \end{figure}
578    
579    \begin{figure}
580    \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-hi-diag.eps}}
581    \caption{
582    Comparison of advection schemes in two dimensions; diagonal advection
583    of a resolved Guassian feature. Courant number is 0.47 with
584    30$\times$30 points and solutions are shown for T=1/2. White lines
585    indicate zero crossings and initial maximum values (ie. the presence
586    of false extrema).  The left column shows the second order schemes;
587    top) centered second order with Adams-Bashforth, middle) Lax-Wendroff
588    and bottom) Superbee flux limited. The middle column shows the third
589    order schemes; top) upwind biased third order with Adams-Bashforth,
590    middle) third order direct space-time method and bottom) the same with
591    flux limiting. The top right panel shows the centered fourth order
592    scheme with Adams-Bashforth and right middle panel shows a fourth
593    order variant on the DST method. Bottom right panel shows the Superbee
594    flux limiter (second order) applied independantly in each direction
595    (method of lines).
596    \label{fig:advect-2d-hi-diag}
597    }
598    \end{figure}
599    
600    
601    
602  In many of the aforementioned advection schemes the behaviour in  In many of the aforementioned advection schemes the behaviour in
603  multiple dimensions is not necessarily as good as the one dimensional  multiple dimensions is not necessarily as good as the one dimensional
604  behaviour. For instance, a shape preserving monotonic scheme in one  behaviour. For instance, a shape preserving monotonic scheme in one
# Line 408  So that the over all time-stepping looks Line 632  So that the over all time-stepping looks
632  \begin{equation}  \begin{equation}
633  \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)
634  \end{equation}  \end{equation}
635    
636    \fbox{ \begin{minipage}{4.75in}
637    {\em S/R GAD\_ADVECTION} ({\em gad\_advection.F})
638    
639    $\tau$: {\bf Tracer} (argument)
640    
641    $G^{n+1/2}_{adv}$: {\bf Gtracer} (argument)
642    
643    $F_x, F_y, F_r$: {\bf af} (local)
644    
645    $U$: {\bf uTrans} (local)
646    
647    $V$: {\bf vTrans} (local)
648    
649    $W$: {\bf rTrans} (local)
650    
651    \end{minipage} }
652    
653    
654    \section{Comparison of advection schemes}
655    
656    Figs.~\ref{fig:advect-2d-lo-diag}, \ref{fig:advect-2d-mid-diag} and
657    \ref{fig:advect-2d-hi-diag} show solutions to a simple diagonal
658    advection problem using a selection of schemes for low, moderate and
659    high Courant numbers, respectively. The top row shows the linear
660    schemes, integrated with the Adams-Bashforth method. Theses schemes
661    are clearly unstable for the high Courant number and weakly unstable
662    for the moderate Courant number. The presence of false extrema is very
663    apparent for all Courant numbers. The middle row shows solutions
664    obtained with the unlimited but multi-dimensional schemes. These
665    solutions also exhibit false extrema though the pattern now shows
666    symmetry due to the multi-dimensional scheme. Also, the schemes are
667    stable at high Courant number where the linear schemes weren't. The
668    bottom row (left and middle) shows the limited schemes and most
669    obvious is the absence of false extrema. The accuracy and stability of
670    the unlimited non-linear schemes is retained at high Courant number
671    but at low Courant number the tendancy is to loose amplitude in sharp
672    peaks due to diffusion. The one dimensional tests shown in
673    Figs.~\ref{fig:advect-1d-lo} and \ref{fig:advect-1d-hi} showed this
674    phenomenum.
675    
676    Finally, the bottom left and right panels use the same advection
677    scheme but the right does not use the mutli-dimensional method. At low
678    Courant number this appears to not matter but for moderate Courant
679    number severe distortion of the feature is apparent. Moreoever, the
680    stability of the multi-dimensional scheme is determined by the maximum
681    Courant number applied of each dimension while the stability of the
682    method of lines is determined by the sum. Hence, in the high Courant
683    number plot, the scheme is unstable.
684    
685    With many advection schemes implemented in the code two questions
686    arise: ``Which scheme is best?'' and ``Why don't you just offer the
687    best advection scheme?''. Unfortunately, no one advection scheme is
688    ``the best'' for all particular applications and for new applications
689    it is often a matter of trial to determine which is most
690    suitable. Here are some guidelines but these are not the rule;
691    \begin{itemize}
692    \item If you have a coarsely resolved model, using a
693    positive or upwind biased scheme will introduce significant diffusion
694    to the solution and using a centered higher order scheme will
695    introduce more noise. In this case, simplest may be best.
696    \item If you have a high resolution model, using a higher order
697    scheme will give a more accurate solution but scale-selective
698    diffusion might need to be employed. The flux limited methods
699    offer similar accuracy in this regime.
700    \item If your solution has shocks or propagatin fronts then a
701    flux limited scheme is almost essential.
702    \item If your time-step is limited by advection, the multi-dimensional
703    non-linear schemes have the most stablility (upto Courant number 1).
704    \item If you need to know how much diffusion/dissipation has occured you
705    will have a lot of trouble figuring it out with a non-linear method.
706    \item The presence of false extrema is unphysical and this alone is the
707    strongest argument for using a positive scheme.
708    \end{itemize}

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

  ViewVC Help
Powered by ViewVC 1.1.22