/[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.6 by cnh, Wed Oct 24 14:24:06 2001 UTC revision 1.14 by afe, Tue Mar 23 15:29:40 2004 UTC
# Line 2  Line 2 
2  % $Name$  % $Name$
3    
4  \section{Tracer equations}  \section{Tracer equations}
5  \label{sec:tracer_equations}  \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
# Line 15  part of the tracer equations and the var Line 15  part of the tracer equations and the var
15  described here.  described here.
16    
17  \subsection{Time-stepping of tracers: ABII}  \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  The default advection scheme is the centered second order method which
21  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 42  only affects the surface layer since the Line 43  only affects the surface layer since the
43  everywhere else. This term is therefore referred to as the surface  everywhere else. This term is therefore referred to as the surface
44  correction term. Global conservation is not possible using the  correction term. Global conservation is not possible using the
45  flux-form (as here) and a linearized free-surface  flux-form (as here) and a linearized free-surface
46  (\cite{Griffies00,Campin02}).  (\cite{griffies:00,campin:02}).
47    
48  The continuity equation can be recovered by setting  The continuity equation can be recovered by setting
49  $G_{diff}=G_{forc}=0$ and $\tau=1$.  $G_{diff}=G_{forc}=0$ and $\tau=1$.
50    
51  The driver routine that calls the routines to calculate tendancies are  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  {\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  (moisture), respectively. These in turn call a generic advection
54  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
55  flow field and relevent tracer as arguments and returns the collective  flow field and relevant tracer as arguments and returns the collective
56  tendancy due to advection and diffusion. Forcing is add subsequently  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 tendancy  in {\em S/R CALC\_GT} or {\em S/R CALC\_GS} to the same tendency
58  array.  array.
59    
60  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
# Line 67  $F_r$: {\bf fVerT} (argument) Line 68  $F_r$: {\bf fVerT} (argument)
68    
69  \end{minipage} }  \end{minipage} }
70    
71  The space and time discretizations are treated seperately (method of  The space and time discretization are treated separately (method of
72  lines). Tendancies are calculated at time levels $n$ and $n-1$ and  lines). Tendencies are calculated at time levels $n$ and $n-1$ and
73  extrapolated to $n+1/2$ using the Adams-Bashforth method:  extrapolated to $n+1/2$ using the Adams-Bashforth method:
74  \marginpar{$\epsilon$: {\bf AB\_eps}}  \marginpar{$\epsilon$: {\bf AB\_eps}}
75  \begin{equation}  \begin{equation}
# Line 76  G^{(n+1/2)} = Line 77  G^{(n+1/2)} =
77  (\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)}
78  \end{equation}  \end{equation}
79  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
80  step $n$. The tendancy at $n-1$ is not re-calculated but rather the  step $n$. The tendency at $n-1$ is not re-calculated but rather the
81  tendancy at $n$ is stored in a global array for later re-use.  tendency at $n$ is stored in a global array for later re-use.
82    
83  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
84  {\em S/R ADAMS\_BASHFORTH2} ({\em model/src/adams\_bashforth2.F})  {\em S/R ADAMS\_BASHFORTH2} ({\em model/src/adams\_bashforth2.F})
# Line 92  $\epsilon$: {\bf ABeps} (PARAMS.h) Line 93  $\epsilon$: {\bf ABeps} (PARAMS.h)
93    
94  \end{minipage} }  \end{minipage} }
95    
96  The tracers are stepped forward in time using the extrapolated tendancy:  The tracers are stepped forward in time using the extrapolated tendency:
97  \begin{equation}  \begin{equation}
98  \tau^{(n+1)} = \tau^{(n)} + \Delta t G^{(n+1/2)}  \tau^{(n+1)} = \tau^{(n)} + \Delta t G^{(n+1/2)}
99  \end{equation}  \end{equation}
# Line 112  $\Delta t$: {\bf deltaTtracer} (PARAMS.h Line 113  $\Delta t$: {\bf deltaTtracer} (PARAMS.h
113  \end{minipage} }  \end{minipage} }
114    
115  Strictly speaking the ABII scheme should be applied only to the  Strictly speaking the ABII scheme should be applied only to the
116  advection terms. However, this scheme is only used in conjuction with  advection terms. However, this scheme is only used in conjunction with
117  the standard second, third and fourth order advection  the standard second, third and fourth order advection
118  schemes. Selection of any other advection scheme disables  schemes. Selection of any other advection scheme disables
119  Adams-Bashforth for tracers so that explicit diffusion and forcing use  Adams-Bashforth for tracers so that explicit diffusion and forcing use
# Line 122  the forward method. Line 123  the forward method.
123    
124    
125  \section{Linear advection schemes}  \section{Linear advection schemes}
126    \label{sect:tracer-advection}
127    \begin{rawhtml}
128    <!-- CMIREDIR:linear_advection_schemes -->
129    \end{rawhtml}
130    
131  \begin{figure}  \begin{figure}
132  \resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-lo.eps}}  \resizebox{5.5in}{!}{\includegraphics{part2/advect-1d-lo.eps}}
# Line 170  commonly used in the field and most fami Line 175  commonly used in the field and most fami
175  \subsection{Centered second order advection-diffusion}  \subsection{Centered second order advection-diffusion}
176    
177  The basic discretization, centered second order, is the default. It is  The basic discretization, centered second order, is the default. It is
178  designed to be consistant with the continuity equation to facilitate  designed to be consistent with the continuity equation to facilitate
179  conservation properties analogous to the continuum. However, centered  conservation properties analogous to the continuum. However, centered
180  second order advection is notoriously noisey and must be used in  second order advection is notoriously noisy and must be used in
181  conjuction with some finite amount of diffusion to produce a sensible  conjunction with some finite amount of diffusion to produce a sensible
182  solution.  solution.
183    
184  The advection operator is discretized:  The advection operator is discretized:
# Line 199  W & = & {\cal A}_c w Line 204  W & = & {\cal A}_c w
204    
205  For non-divergent flow, this discretization can be shown to conserve  For non-divergent flow, this discretization can be shown to conserve
206  the tracer both locally and globally and to globally conserve tracer  the tracer both locally and globally and to globally conserve tracer
207  variance, $\tau^2$. The proof is given in \cite{Adcroft95,Adcroft97}.  variance, $\tau^2$. The proof is given in \cite{adcroft:95,adcroft:97}.
208    
209  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
210  {\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 247  F_r & = & W \overline{\tau - \frac{1}{6} Line 252  F_r & = & W \overline{\tau - \frac{1}{6}
252  \end{eqnarray}  \end{eqnarray}
253    
254  At boundaries, $\delta_{\hat{n}} \tau$ is set to zero allowing  At boundaries, $\delta_{\hat{n}} \tau$ is set to zero allowing
255  $\delta_{nn}$ to be evaluated. We are currently examing the accuracy  $\delta_{nn}$ to be evaluated. We are currently examine the accuracy
256  of this boundary condition and the effect on the solution.  of this boundary condition and the effect on the solution.
257    
258  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
# Line 281  $\tau$: {\bf tracer} (argument) Line 286  $\tau$: {\bf tracer} (argument)
286    
287  Centered fourth order advection is formally the most accurate scheme  Centered fourth order advection is formally the most accurate scheme
288  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
289  simultation where dynamical scales are well resolved. However, the  simulation where dynamical scales are well resolved. However, the
290  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
291  used with some finite amount of diffusion. Bi-harmonic is recommended  used with some finite amount of diffusion. Bi-harmonic is recommended
292  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
293  well resolved gradient the fourth order scheme worked so hard to  well resolved gradient the fourth order scheme worked so hard to
# Line 296  F_r & = & W \overline{\tau - \frac{1}{6} Line 301  F_r & = & W \overline{\tau - \frac{1}{6}
301  \end{eqnarray}  \end{eqnarray}
302    
303  As for the third order scheme, the best discretization near boundaries  As for the third order scheme, the best discretization near boundaries
304  is under investigation but currenlty $\delta_i \tau=0$ on a boundary.  is under investigation but currently $\delta_i \tau=0$ on a boundary.
305    
306  \fbox{ \begin{minipage}{4.75in}  \fbox{ \begin{minipage}{4.75in}
307  {\em S/R GAD\_C4\_ADV\_X} ({\em gad\_c4\_adv\_x.F})  {\em S/R GAD\_C4\_ADV\_X} ({\em gad\_c4\_adv\_x.F})
# Line 386  r = \frac{ \tau_{i-1} - \tau_{i-2} }{ \t Line 391  r = \frac{ \tau_{i-1} - \tau_{i-2} }{ \t
391  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
392  \end{eqnarray}  \end{eqnarray}
393  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
394  only provide the Superbee limiter \cite{Roe85}:  only provide the Superbee limiter \cite{roe:85}:
395  \begin{equation}  \begin{equation}
396  \psi(r) = \max[0,\min[1,2r],\min[2,r]]  \psi(r) = \max[0,\min[1,2r],\min[2,r]]
397  \end{equation}  \end{equation}
# Line 422  $\tau$: {\bf tracer} (argument) Line 427  $\tau$: {\bf tracer} (argument)
427  \subsection{Third order direct space time}  \subsection{Third order direct space time}
428    
429  The direct-space-time method deals with space and time discretization  The direct-space-time method deals with space and time discretization
430  together (other methods that treat space and time seperately are known  together (other methods that treat space and time separately are known
431  collectively as the ``Method of Lines''). The Lax-Wendroff scheme  collectively as the ``Method of Lines''). The Lax-Wendroff scheme
432  falls into this category; it adds sufficient diffusion to a second  falls into this category; it adds sufficient diffusion to a second
433  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 440  where Line 445  where
445  d_1 & = & \frac{1}{6} ( 2 - |c| ) ( 1 - |c| ) \\  d_1 & = & \frac{1}{6} ( 2 - |c| ) ( 1 - |c| ) \\
446  d_2 & = & \frac{1}{6} ( 1 - |c| ) ( 1 + |c| )  d_2 & = & \frac{1}{6} ( 1 - |c| ) ( 1 + |c| )
447  \end{eqnarray}  \end{eqnarray}
448  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
449  as the Courant number, $c$, vanishes. In this limit, the conventional  as the Courant number, $c$, vanishes. In this limit, the conventional
450  third order upwind method is recovered. For finite Courant number, the  third order upwind method is recovered. For finite Courant number, the
451  deviations from the linear method are analogous to the diffusion added  deviations from the linear method are analogous to the diffusion added
# Line 448  to centered second order advection in th Line 453  to centered second order advection in th
453    
454  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
455  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
456  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
457  the accuracy increases with the Courant number! For low Courant  the accuracy increases with the Courant number! For low Courant
458  number, DST3 produces very similar results (indistinguishable in  number, DST3 produces very similar results (indistinguishable in
459  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 538  $\tau$: {\bf tracer} (argument) Line 543  $\tau$: {\bf tracer} (argument)
543  \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-lo-diag.eps}}  \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-lo-diag.eps}}
544  \caption{  \caption{
545  Comparison of advection schemes in two dimensions; diagonal advection  Comparison of advection schemes in two dimensions; diagonal advection
546  of a resolved Guassian feature. Courant number is 0.01 with  of a resolved Gaussian feature. Courant number is 0.01 with
547  30$\times$30 points and solutions are shown for T=1/2. White lines  30$\times$30 points and solutions are shown for T=1/2. White lines
548  indicate zero crossing (ie. the presence of false minima).  The left  indicate zero crossing (ie. the presence of false minima).  The left
549  column shows the second order schemes; top) centered second order with  column shows the second order schemes; top) centered second order with
# Line 549  space-time method and bottom) the same w Line 554  space-time method and bottom) the same w
554  right panel shows the centered fourth order scheme with  right panel shows the centered fourth order scheme with
555  Adams-Bashforth and right middle panel shows a fourth order variant on  Adams-Bashforth and right middle panel shows a fourth order variant on
556  the DST method. Bottom right panel shows the Superbee flux limiter  the DST method. Bottom right panel shows the Superbee flux limiter
557  (second order) applied independantly in each direction (method of  (second order) applied independently in each direction (method of
558  lines).  lines).
559  \label{fig:advect-2d-lo-diag}  \label{fig:advect-2d-lo-diag}
560  }  }
# Line 559  lines). Line 564  lines).
564  \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-mid-diag.eps}}  \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-mid-diag.eps}}
565  \caption{  \caption{
566  Comparison of advection schemes in two dimensions; diagonal advection  Comparison of advection schemes in two dimensions; diagonal advection
567  of a resolved Guassian feature. Courant number is 0.27 with  of a resolved Gaussian feature. Courant number is 0.27 with
568  30$\times$30 points and solutions are shown for T=1/2. White lines  30$\times$30 points and solutions are shown for T=1/2. White lines
569  indicate zero crossing (ie. the presence of false minima).  The left  indicate zero crossing (ie. the presence of false minima).  The left
570  column shows the second order schemes; top) centered second order with  column shows the second order schemes; top) centered second order with
# Line 570  space-time method and bottom) the same w Line 575  space-time method and bottom) the same w
575  right panel shows the centered fourth order scheme with  right panel shows the centered fourth order scheme with
576  Adams-Bashforth and right middle panel shows a fourth order variant on  Adams-Bashforth and right middle panel shows a fourth order variant on
577  the DST method. Bottom right panel shows the Superbee flux limiter  the DST method. Bottom right panel shows the Superbee flux limiter
578  (second order) applied independantly in each direction (method of  (second order) applied independently in each direction (method of
579  lines).  lines).
580  \label{fig:advect-2d-mid-diag}  \label{fig:advect-2d-mid-diag}
581  }  }
# Line 580  lines). Line 585  lines).
585  \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-hi-diag.eps}}  \resizebox{5.5in}{!}{\includegraphics{part2/advect-2d-hi-diag.eps}}
586  \caption{  \caption{
587  Comparison of advection schemes in two dimensions; diagonal advection  Comparison of advection schemes in two dimensions; diagonal advection
588  of a resolved Guassian feature. Courant number is 0.47 with  of a resolved Gaussian feature. Courant number is 0.47 with
589  30$\times$30 points and solutions are shown for T=1/2. White lines  30$\times$30 points and solutions are shown for T=1/2. White lines
590  indicate zero crossings and initial maximum values (ie. the presence  indicate zero crossings and initial maximum values (ie. the presence
591  of false extrema).  The left column shows the second order schemes;  of false extrema).  The left column shows the second order schemes;
# Line 591  middle) third order direct space-time me Line 596  middle) third order direct space-time me
596  flux limiting. The top right panel shows the centered fourth order  flux limiting. The top right panel shows the centered fourth order
597  scheme with Adams-Bashforth and right middle panel shows a fourth  scheme with Adams-Bashforth and right middle panel shows a fourth
598  order variant on the DST method. Bottom right panel shows the Superbee  order variant on the DST method. Bottom right panel shows the Superbee
599  flux limiter (second order) applied independantly in each direction  flux limiter (second order) applied independently in each direction
600  (method of lines).  (method of lines).
601  \label{fig:advect-2d-hi-diag}  \label{fig:advect-2d-hi-diag}
602  }  }
# Line 599  flux limiter (second order) applied inde Line 604  flux limiter (second order) applied inde
604    
605    
606    
607  In many of the aforementioned advection schemes the behaviour in  In many of the aforementioned advection schemes the behavior in
608  multiple dimensions is not necessarily as good as the one dimensional  multiple dimensions is not necessarily as good as the one dimensional
609  behaviour. For instance, a shape preserving monotonic scheme in one  behavior. For instance, a shape preserving monotonic scheme in one
610  dimension can have severe shape distortion in two dimensions if the  dimension can have severe shape distortion in two dimensions if the
611  two components of horizontal fluxes are treated independently. There  two components of horizontal fluxes are treated independently. There
612  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 622  as if in one dimension: Line 627  as if in one dimension:
627  \end{eqnarray}  \end{eqnarray}
628    
629  In order to incorporate this method into the general model algorithm,  In order to incorporate this method into the general model algorithm,
630  we compute the effective tendancy rather than update the tracer so  we compute the effective tendency rather than update the tracer so
631  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
632  not the updated $n+3/3$ quantities:  not the updated $n+3/3$ quantities:
633  \begin{equation}  \begin{equation}
# Line 668  stable at high Courant number where the Line 673  stable at high Courant number where the
673  bottom row (left and middle) shows the limited schemes and most  bottom row (left and middle) shows the limited schemes and most
674  obvious is the absence of false extrema. The accuracy and stability of  obvious is the absence of false extrema. The accuracy and stability of
675  the unlimited non-linear schemes is retained at high Courant number  the unlimited non-linear schemes is retained at high Courant number
676  but at low Courant number the tendancy is to loose amplitude in sharp  but at low Courant number the tendency is to loose amplitude in sharp
677  peaks due to diffusion. The one dimensional tests shown in  peaks due to diffusion. The one dimensional tests shown in
678  Figs.~\ref{fig:advect-1d-lo} and \ref{fig:advect-1d-hi} showed this  Figs.~\ref{fig:advect-1d-lo} and \ref{fig:advect-1d-hi} showed this
679  phenomenum.  phenomenon.
680    
681  Finally, the bottom left and right panels use the same advection  Finally, the bottom left and right panels use the same advection
682  scheme but the right does not use the mutli-dimensional method. At low  scheme but the right does not use the multi-dimensional method. At low
683  Courant number this appears to not matter but for moderate Courant  Courant number this appears to not matter but for moderate Courant
684  number severe distortion of the feature is apparent. Moreoever, the  number severe distortion of the feature is apparent. Moreover, the
685  stability of the multi-dimensional scheme is determined by the maximum  stability of the multi-dimensional scheme is determined by the maximum
686  Courant number applied of each dimension while the stability of the  Courant number applied of each dimension while the stability of the
687  method of lines is determined by the sum. Hence, in the high Courant  method of lines is determined by the sum. Hence, in the high Courant
# Line 697  introduce more noise. In this case, simp Line 702  introduce more noise. In this case, simp
702  scheme will give a more accurate solution but scale-selective  scheme will give a more accurate solution but scale-selective
703  diffusion might need to be employed. The flux limited methods  diffusion might need to be employed. The flux limited methods
704  offer similar accuracy in this regime.  offer similar accuracy in this regime.
705  \item If your solution has shocks or propagatin fronts then a  \item If your solution has shocks or propagating fronts then a
706  flux limited scheme is almost essential.  flux limited scheme is almost essential.
707  \item If your time-step is limited by advection, the multi-dimensional  \item If your time-step is limited by advection, the multi-dimensional
708  non-linear schemes have the most stablility (upto Courant number 1).  non-linear schemes have the most stability (up to Courant number 1).
709  \item If you need to know how much diffusion/dissipation has occured you  \item If you need to know how much diffusion/dissipation has occurred you
710  will have a lot of trouble figuring it out with a non-linear method.  will have a lot of trouble figuring it out with a non-linear method.
711  \item The presence of false extrema is unphysical and this alone is the  \item The presence of false extrema is non-physical and this alone is the
712  strongest argument for using a positive scheme.  strongest argument for using a positive scheme.
713  \end{itemize}  \end{itemize}

Legend:
Removed from v.1.6  
changed lines
  Added in v.1.14

  ViewVC Help
Powered by ViewVC 1.1.22