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

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

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

revision 1.19 by edhill, Sat Oct 16 03:40:12 2004 UTC revision 1.21 by jmc, Tue Apr 4 20:16:39 2006 UTC
# Line 225  The rigid-lid approximation can be easil Line 225  The rigid-lid approximation can be easil
225  of the free-surface equation which can be written:  of the free-surface equation which can be written:
226  \begin{equation}  \begin{equation}
227  \partial_t \eta + \partial_x H \widehat{u} + \partial_y H \widehat{v} = P-E+R  \partial_t \eta + \partial_x H \widehat{u} + \partial_y H \widehat{v} = P-E+R
228  \label{eq:linear-free-surface=P-E+R}  \label{eq:linear-free-surface=P-E}
229  \end{equation}  \end{equation}
230  which differs from the depth integrated continuity equation with  which differs from the depth integrated continuity equation with
231  rigid-lid (\ref{eq:rigid-lid-continuity}) by the time-dependent term  rigid-lid (\ref{eq:rigid-lid-continuity}) by the time-dependent term
# Line 233  and fresh-water source term. Line 233  and fresh-water source term.
233    
234  Equation \ref{eq:discrete-time-cont-rigid-lid} in the rigid-lid  Equation \ref{eq:discrete-time-cont-rigid-lid} in the rigid-lid
235  pressure method is then replaced by the time discretization of  pressure method is then replaced by the time discretization of
236  \ref{eq:linear-free-surface=P-E+R} which is:  \ref{eq:linear-free-surface=P-E} which is:
237  \begin{equation}  \begin{equation}
238  \eta^{n+1}  \eta^{n+1}
239  + \Delta t \partial_x H \widehat{u^{n+1}}  + \Delta t \partial_x H \widehat{u^{n+1}}
240  + \Delta t \partial_y H \widehat{v^{n+1}}  + \Delta t \partial_y H \widehat{v^{n+1}}
241  =  =
242  \eta^{n}  \eta^{n}
243  + \Delta t ( P - E + R )  + \Delta t ( P - E )
244  \label{eq:discrete-time-backward-free-surface}  \label{eq:discrete-time-backward-free-surface}
245  \end{equation}  \end{equation}
246  where the use of flow at time level $n+1$ makes the method implicit  where the use of flow at time level $n+1$ makes the method implicit
# Line 255  As for the rigid-lid pressure method, eq Line 255  As for the rigid-lid pressure method, eq
255  \begin{eqnarray}  \begin{eqnarray}
256  u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-backward-free-surface} \\  u^{*} & = & u^{n} + \Delta t G_u^{(n+1/2)} \label{eq:ustar-backward-free-surface} \\
257  v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-backward-free-surface} \\  v^{*} & = & v^{n} + \Delta t G_v^{(n+1/2)} \label{eq:vstar-backward-free-surface} \\
258  \eta^* & = & \epsilon_{fs} \left( \eta^{n} +P-E+R \right)- \Delta t  \eta^* & = & \epsilon_{fs} \left( \eta^{n} + \Delta t (P-E) \right)- \Delta t
259    \partial_x H \widehat{u^{*}}    \partial_x H \widehat{u^{*}}
260  + \partial_y H \widehat{v^{*}}  + \partial_y H \widehat{v^{*}}
261  \\  \\
262    \partial_x g H \partial_x \eta^{n+1}    \partial_x g H \partial_x \eta^{n+1}
263  + \partial_y g H \partial_y \eta^{n+1}  & + & \partial_y g H \partial_y \eta^{n+1}
264  - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}   - \frac{\epsilon_{fs} \eta^{n+1}}{\Delta t^2}
265  & = &   =
266  - \frac{\eta^*}{\Delta t^2}  - \frac{\eta^*}{\Delta t^2}
267  \label{eq:elliptic-backward-free-surface}  \label{eq:elliptic-backward-free-surface}
268  \\  \\
# Line 442  FORWARD\_STEP \\ Line 442  FORWARD\_STEP \\
442  \>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ (\ref{eq:Gt-n-sync})\\  \>\>\> GAD\_CALC\_RHS \` $G_\theta^n = G_\theta( u, \theta^n )$ (\ref{eq:Gt-n-sync})\\
443  \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\  \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\
444  \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-sync}) \\  \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-sync}) \\
445  \>\> TIMESTEP\_TRACER \` $\tau^*$ (\ref{eq:tstar-sync}) \\  \>\> TIMESTEP\_TRACER \` $\theta^*$ (\ref{eq:tstar-sync}) \\
446  \>\> IMPLDIFF \` $\tau^{(n+1)}$ (\ref{eq:t-n+1-sync}) \\  \>\> IMPLDIFF \` $\theta^{(n+1)}$ (\ref{eq:t-n+1-sync}) \\
447  \> DYNAMICS \\  \> DYNAMICS \\
448  \>\> CALC\_PHI\_HYD \` $\phi_{hyd}^n$ (\ref{eq:phi-hyd-sync}) \\  \>\> CALC\_PHI\_HYD \` $\phi_{hyd}^n$ (\ref{eq:phi-hyd-sync}) \\
449  \>\> MOM\_FLUXFORM or MOM\_VECINV \` $G_{\vec{\bf v}}^n$ (\ref{eq:Gv-n-sync})\\  \>\> MOM\_FLUXFORM or MOM\_VECINV \` $G_{\vec{\bf v}}^n$ (\ref{eq:Gv-n-sync})\\
# Line 465  FORWARD\_STEP \\ Line 465  FORWARD\_STEP \\
465  Calling tree for the overall synchronous algorithm using  Calling tree for the overall synchronous algorithm using
466  Adams-Bashforth time-stepping.  Adams-Bashforth time-stepping.
467  The place where the model geometry  The place where the model geometry
468  ({\em hFac} factors) is updated is added here but is only relevant  ({\bf hFac} factors) is updated is added here but is only relevant
469  for the non-linear free-surface algorithm.  for the non-linear free-surface algorithm.
470  For completeness, the external forcing,  For completeness, the external forcing,
471  ocean and atmospheric physics have been added, although they are mainly  ocean and atmospheric physics have been added, although they are mainly
# Line 547  extrapolate tendencies to $n$ (dashed ar Line 547  extrapolate tendencies to $n$ (dashed ar
547  The hydrostatic pressure/geo-potential $\phi_{hyd}$ is evaluated directly  The hydrostatic pressure/geo-potential $\phi_{hyd}$ is evaluated directly
548  at time level $n$ (vertical arrows) and used with the extrapolated tendencies  at time level $n$ (vertical arrows) and used with the extrapolated tendencies
549  to step forward the flow variables from $n-1/2$ to $n+1/2$ (solid arc-arrow).  to step forward the flow variables from $n-1/2$ to $n+1/2$ (solid arc-arrow).
550  The implicit-in-time operator ${\cal L_{u,v}}$ (vertical arrows) is  The implicit-in-time operator ${\cal L}_{\bf u,v}$ (vertical arrows) is
551  then applied to the previous estimation of the the flow field ($*$-variables)  then applied to the previous estimation of the the flow field ($*$-variables)
552  and yields to the two velocity components $u,v$ at time level $n+1/2$.  and yields to the two velocity components $u,v$ at time level $n+1/2$.
553  These are then used to calculate the advection term (dashed arc-arrow)  These are then used to calculate the advection term (dashed arc-arrow)
# Line 575  thermodynamics solver is delayed from ha Line 575  thermodynamics solver is delayed from ha
575  allowing the use of the most recent velocities to compute  allowing the use of the most recent velocities to compute
576  the advection terms. Once the thermodynamics fields are  the advection terms. Once the thermodynamics fields are
577  updated, the hydrostatic pressure is computed  updated, the hydrostatic pressure is computed
578  to step frowrad the dynamics  to step forwrad the dynamics.
579  Note that the pressure gradient must also be taken out of the  Note that the pressure gradient must also be taken out of the
580  Adams-Bashforth extrapolation. Also, retaining the integer time-levels,  Adams-Bashforth extrapolation. Also, retaining the integer time-levels,
581  $n$ and $n+1$, does not give a user the sense of where variables are  $n$ and $n+1$, does not give a user the sense of where variables are
# Line 583  located in time.  Instead, we re-write t Line 583  located in time.  Instead, we re-write t
583  \ref{eq:Gt-n-sync} to \ref{eq:v-n+1-sync}, annotating the  \ref{eq:Gt-n-sync} to \ref{eq:v-n+1-sync}, annotating the
584  position in time of variables appropriately:  position in time of variables appropriately:
585  \begin{eqnarray}  \begin{eqnarray}
586    \phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n}) dr
587    \label{eq:phi-hyd-staggered} \\
588  \vec{\bf G}_{\vec{\bf v}}^{n-1/2} & = & \vec{\bf G}_{\vec{\bf v}} ( \vec{\bf v}^{n-1/2} )  \vec{\bf G}_{\vec{\bf v}}^{n-1/2} & = & \vec{\bf G}_{\vec{\bf v}} ( \vec{\bf v}^{n-1/2} )
589  \label{eq:Gv-n-staggered} \\  \label{eq:Gv-n-staggered} \\
590  \vec{\bf G}_{\vec{\bf v}}^{(n)} & = & (3/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-1/2} - (1/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-3/2}  \vec{\bf G}_{\vec{\bf v}}^{(n)} & = & (3/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-1/2} - (1/2 + \epsilon_{AB} ) \vec{\bf G}_{\vec{\bf v}}^{n-3/2}
591  \label{eq:Gv-n+5-staggered} \\  \label{eq:Gv-n+5-staggered} \\
 \phi^{n}_{hyd} & = & \int b(\theta^{n},S^{n}) dr  
 \label{eq:phi-hyd-staggered} \\  
592  \vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \Delta t \left( \vec{\bf G}_{\vec{\bf v}}^{(n)} - \nabla \phi_{hyd}^{n} \right)  \vec{\bf v}^{*} & = & \vec{\bf v}^{n-1/2} + \Delta t \left( \vec{\bf G}_{\vec{\bf v}}^{(n)} - \nabla \phi_{hyd}^{n} \right)
593  \label{eq:vstar-staggered} \\  \label{eq:vstar-staggered} \\
594  \vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* )  \vec{\bf v}^{**} & = & {\cal L}_{\vec{\bf v}}^{-1} ( \vec{\bf v}^* )
# Line 613  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil Line 613  G_{\theta,S}^{(n+1/2)} & = & (3/2+\epsil
613  The corresponding calling tree is given in  The corresponding calling tree is given in
614  \ref{fig:call-tree-adams-bashforth-staggered}.  \ref{fig:call-tree-adams-bashforth-staggered}.
615  The staggered algorithm is activated with the run-time flag  The staggered algorithm is activated with the run-time flag
616  {\bf staggerTimeStep=.TRUE.} in parameter file {\em data},  {\bf staggerTimeStep}{\em=.TRUE.} in parameter file {\em data},
617  namelist {\em PARM01}.  namelist {\em PARM01}.
618    
619  \begin{figure}  \begin{figure}
# Line 643  FORWARD\_STEP \\ Line 643  FORWARD\_STEP \\
643       (\ref{eq:Gt-n-staggered})\\       (\ref{eq:Gt-n-staggered})\\
644  \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\  \>\>\> EXTERNAL\_FORCING \` $G_\theta^n = G_\theta^n + {\cal Q}$ \\
645  \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-staggered}) \\  \>\>\> ADAMS\_BASHFORTH2 \` $G_\theta^{(n+1/2)}$ (\ref{eq:Gt-n+5-staggered}) \\
646  \>\> TIMESTEP\_TRACER \` $\tau^*$ (\ref{eq:tstar-staggered}) \\  \>\> TIMESTEP\_TRACER \` $\theta^*$ (\ref{eq:tstar-staggered}) \\
647  \>\> IMPLDIFF \` $\tau^{(n+1)}$ (\ref{eq:t-n+1-staggered}) \\  \>\> IMPLDIFF \` $\theta^{(n+1)}$ (\ref{eq:t-n+1-staggered}) \\
648  \> TRACERS\_CORRECTION\_STEP  \\  \> TRACERS\_CORRECTION\_STEP  \\
649  \>\> CYCLE\_TRACER \` $\theta^{n+1}$ \\  \>\> CYCLE\_TRACER \` $\theta^{n+1}$ \\
650  \>\> FILTER \` Shapiro Filter, Zonal Filter (FFT) \\  \>\> FILTER \` Shapiro Filter, Zonal Filter (FFT) \\
# Line 654  FORWARD\_STEP \\ Line 654  FORWARD\_STEP \\
654  Calling tree for the overall staggered algorithm using  Calling tree for the overall staggered algorithm using
655  Adams-Bashforth time-stepping.  Adams-Bashforth time-stepping.
656  The place where the model geometry  The place where the model geometry
657  ({\em hFac} factors) is updated is added here but is only relevant  ({\bf hFac} factors) is updated is added here but is only relevant
658  for the non-linear free-surface algorithm.  for the non-linear free-surface algorithm.
659  }  }
660  \label{fig:call-tree-adams-bashforth-staggered}  \label{fig:call-tree-adams-bashforth-staggered}
# Line 775  $w^{n+1}$. Line 775  $w^{n+1}$.
775    
776    
777  \section{Variants on the Free Surface}  \section{Variants on the Free Surface}
778    \label{sect:free-surface}
779    
780  We now describe the various formulations of the free-surface that  We now describe the various formulations of the free-surface that
781  include non-linear forms, implicit in time using Crank-Nicholson,  include non-linear forms, implicit in time using Crank-Nicholson,
# Line 813  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h) Line 814  $\eta^{n+1}$: {\bf etaN} (\em DYNVARS.h)
814    
815    
816  Once ${\eta}^{n+1}$ has been found, substituting into  Once ${\eta}^{n+1}$ has been found, substituting into
817  \ref{eq:discrete-time-u,eq:discrete-time-v} yields $\vec{\bf v}^{n+1}$ if the model is  \ref{eq:discrete-time-u}, \ref{eq:discrete-time-v} yields $\vec{\bf v}^{n+1}$
818  hydrostatic ($\epsilon_{nh}=0$):  if the model is hydrostatic ($\epsilon_{nh}=0$):
819  $$  $$
820  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
821  - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}  - \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
# Line 887  at the same point in the code. Line 888  at the same point in the code.
888    
889    
890  \subsection{Crank-Nickelson barotropic time stepping}  \subsection{Crank-Nickelson barotropic time stepping}
891    \label{sect:freesurf-CrankNick}
892    
893  The full implicit time stepping described previously is  The full implicit time stepping described previously is
894  unconditionally stable but damps the fast gravity waves, resulting in  unconditionally stable but damps the fast gravity waves, resulting in
# Line 901  stable, Crank-Nickelson scheme; $(\beta, Line 903  stable, Crank-Nickelson scheme; $(\beta,
903  corresponds to the forward - backward scheme that conserves energy but is  corresponds to the forward - backward scheme that conserves energy but is
904  only stable for small time steps.\\  only stable for small time steps.\\
905  In the code, $\beta,\gamma$ are defined as parameters, respectively  In the code, $\beta,\gamma$ are defined as parameters, respectively
906  {\it implicSurfPress}, {\it implicDiv2DFlow}. They are read from  {\bf implicSurfPress}, {\bf implicDiv2DFlow}. They are read from
907  the main data file "{\it data}" and are set by default to 1,1.  the main parameter file "{\em data}" and are set by default to 1,1.
908    
909  Equations \ref{eq:ustar-backward-free-surface} --  Equations \ref{eq:ustar-backward-free-surface} --
910  \ref{eq:vn+1-backward-free-surface} are modified as follows:  \ref{eq:vn+1-backward-free-surface} are modified as follows:
911  $$  \begin{eqnarray*}
912  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }  \frac{ \vec{\bf v}^{n+1} }{ \Delta t }
913  + {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]  + {\bf \nabla}_h b_s [ \beta {\eta}^{n+1} + (1-\beta) {\eta}^{n} ]
914  + \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1}  + \epsilon_{nh} {\bf \nabla}_h {\phi'_{nh}}^{n+1}
915   = \frac{ \vec{\bf v}^* }{ \Delta t }   = \frac{ \vec{\bf v}^* }{ \Delta t }
916  $$  \end{eqnarray*}
917  $$  \begin{eqnarray}
918  \epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}  \epsilon_{fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
919  + {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}  + {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o}
920  [ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr  [ \gamma \vec{\bf v}^{n+1} + (1-\gamma) \vec{\bf v}^{n}] dr
921  = \epsilon_{fw} (P-E)  = \epsilon_{fw} (P-E)
922  $$  \label{eq:eta-n+1-CrankNick}
923    \end{eqnarray}
924  where:  where:
925  \begin{eqnarray*}  \begin{eqnarray*}
926  \vec{\bf v}^* & = &  \vec{\bf v}^* & = &
# Line 939  $$ Line 942  $$
942  {\bf \nabla}_h {\eta}^{n+1}  {\bf \nabla}_h {\eta}^{n+1}
943  = {\eta}^*  = {\eta}^*
944  $$  $$
945  and then to compute (correction step):  and then to compute ({\em CORRECTION\_STEP}):
946  $$  $$
947  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}  \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
948  - \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}  - \beta \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
949  $$  $$
950    
951  The non-hydrostatic part is solved as described previously.  %The non-hydrostatic part is solved as described previously.
952    
953  Note that:  \noindent
954    Notes:
955  \begin{enumerate}  \begin{enumerate}
956    \item The RHS term of equation \ref{eq:eta-n+1-CrankNick}
957    corresponds the contribution of fresh water flux (P-E)
958    to the free-surface variations ($\epsilon_{fw}=1$,
959    {\bf useRealFreshWater}{\em=TRUE} in parameter file {\em data}).
960    In order to remain consistent with the tracer equation, specially in
961    the non-linear free-surface formulation, this term is also
962    affected by the Crank-Nickelson time stepping. The RHS reads:
963    $\epsilon_{fw} ( \gamma (P-E)^{n+1/2} + (1-\gamma) (P-E)^{n-1/2} )$
964  \item The non-hydrostatic part of the code has not yet been  \item The non-hydrostatic part of the code has not yet been
965  updated, so that this option cannot be used with $(\beta,\gamma) \neq (1,1)$.  updated, and therefore cannot be used with $(\beta,\gamma) \neq (1,1)$.
966  \item The stability criteria with Crank-Nickelson time stepping  \item The stability criteria with Crank-Nickelson time stepping
967  for the pure linear gravity wave problem in cartesian coordinates is:  for the pure linear gravity wave problem in cartesian coordinates is:
968  \begin{itemize}  \begin{itemize}

Legend:
Removed from v.1.19  
changed lines
  Added in v.1.21

  ViewVC Help
Powered by ViewVC 1.1.22