Changeset 2211
- Timestamp:
- 2010-10-12T12:08:06+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_STP.tex
r2164 r2211 22 22 Having defined the continuous equations in Chap.~\ref{PE}, we need now to choose 23 23 a time discretization. In the present chapter, we provide a general description of the \NEMO 24 time stepping strategy and itsconsequences for the order in which the equations are24 time stepping strategy and the consequences for the order in which the equations are 25 25 solved. 26 26 … … 40 40 corresponding time evolution equation; $\rdt$ is the time step; and the 41 41 superscripts indicate the time at which a quantity is evaluated. Each term of the 42 RHS is evaluated at a specific time step (s)depending on the physics with which42 RHS is evaluated at a specific time step depending on the physics with which 43 43 it is associated. 44 44 45 45 The choice of the time step used for this evaluation is discussed below as 46 well as the implications in terms ofstarting or restarting a model47 simulation. Note that the time stepping is generally performed in a one step46 well as the implications for starting or restarting a model 47 simulation. Note that the time stepping calculation is generally performed in a single 48 48 operation. With such a complex and nonlinear system of equations it would be 49 49 dangerous to let a prognostic variable evolve in time for each term separately. 50 50 51 The three level scheme requires three arrays for each prognostic variable s.52 For each variable $x$ there is $x_b$ (before) and $x_n$ (now). The third array,51 The three level scheme requires three arrays for each prognostic variable. 52 For each variable $x$ there is $x_b$ (before), $x_n$ (now) and $x_a$. The third array, 53 53 although referred to as $x_a$ (after) in the code, is usually not the variable at 54 54 the after time step; but rather it is used to store the time derivative (RHS in 55 55 \eqref{Eq_STP}) prior to time-stepping the equation. Generally, the time 56 stepping is performed once at each time step in \mdl{tranxt} and \mdl{dynnxt}57 modules, except for implicit vertical diffusion or sea surface height when58 time-splitting options are used.56 stepping is performed once at each time step in the \mdl{tranxt} and \mdl{dynnxt} 57 modules, except when using implicit vertical diffusion or calculating sea surface height 58 in which case time-splitting options are used. 59 59 60 60 % ------------------------------------------------------------------------------------------------------------- … … 64 64 \label{STP_leap_frog} 65 65 66 The time stepping used for non-diffusive processes is the well-known leapfrog 67 scheme \citep{Mesinger_Arakawa_Bk76}. It is a time centred scheme, $i.e.$ 66 The time stepping used for processes other than diffusion is the well-known leapfrog 67 scheme \citep{Mesinger_Arakawa_Bk76}. This scheme is widely used for advection 68 processes in low-viscosity fluids. It is a time centred scheme, $i.e.$ 68 69 the RHS in \eqref{Eq_STP} is evaluated at time step $t$, the now time step. 69 It is only used for non-diffusive terms, that ismomentum and tracer advection,70 pressure gradient, and Coriolis terms . This scheme is widely used for advective71 processes in low-viscosity fluids.It is an efficient method that achieves70 It may be used for momentum and tracer advection, 71 pressure gradient, and Coriolis terms, but not for diffusion terms. 72 It is an efficient method that achieves 72 73 second-order accuracy with just one right hand side evaluation per time step. 73 74 Moreover, it does not artificially damp linear oscillatory motion nor does it produce 74 75 instability by amplifying the oscillations. These advantages are somewhat diminished 75 76 by the large phase-speed error of the leapfrog scheme, and the unsuitability 76 of leapfrog differencing for the representation of diffusi veand Rayleigh77 of leapfrog differencing for the representation of diffusion and Rayleigh 77 78 damping processes. However, the scheme allows the coexistence of a numerical 78 79 and a physical mode due to its leading third order dispersive error. In other words a 79 80 divergence of odd and even time steps may occur. To prevent it, the leapfrog scheme 80 is often used in association with a Robert-Asselin time filter ( thereafter the LF-RA scheme).81 is often used in association with a Robert-Asselin time filter (hereafter the LF-RA scheme). 81 82 This filter, first designed by \citet{Robert_JMSJ66} and more comprehensively studied 82 83 by \citet{Asselin_MWR72}, is a kind of laplacian diffusion in time that mixes odd and … … 87 88 where the subscript $F$ denotes filtered values and $\gamma$ is the Asselin 88 89 coefficient. $\gamma$ is initialized as \np{rn\_atfp} (namelist parameter). 89 Its default value is \ textit{now} (see \S~\ref{STP_mLF}) \np{rn\_atfp}=$10^{-3}$,90 Its default value is \np{rn\_atfp}=$10^{-3}$ (see \S~\ref{STP_mLF}), 90 91 causing only a weak dissipation of high frequency motions (\citep{Farge1987}). 91 The addition of a time filter d oes, nevertheless, degradethe accuracy of the92 The addition of a time filter degrades the accuracy of the 92 93 calculation from second to first order. However, the second order truncation 93 94 error is proportional to $\gamma$, which is small compared to 1. Therefore, 94 95 the LF-RA is a quasi second order accurate scheme. The LF-RA scheme 95 has beenpreferred to other time differencing schemes such as96 is preferred to other time differencing schemes such as 96 97 predictor corrector or trapezoidal schemes, because the user has an explicit 97 98 and simple control of the magnitude of the time diffusion of the scheme. 98 In association with the 2nd order centred space discretisation of the 99 advective terms in the momentum and tracer equations, it avoids implicit 100 numerical diffusion in both the time and space discretisations of the 101 advective term: they are both set explicitly by the user through the Robert-Asselin 102 filter parameter and the viscous and diffusive coefficients. 99 When used with the 2nd order space centred discretisation of the 100 advection terms in the momentum and tracer equations, LF-RA avoids implicit 101 numerical diffusion: diffusion is set explicitly by the user through the Robert-Asselin 102 filter parameter and the viscosity and diffusion coefficients. 103 103 104 104 % ------------------------------------------------------------------------------------------------------------- … … 109 109 110 110 The leapfrog differencing scheme is unsuitable for the representation of 111 diffusi veand damping processes. For a tendancy $D_x$, representing a112 diffusi veterm or a restoring term to a tracer climatology111 diffusion and damping processes. For a tendancy $D_x$, representing a 112 diffusion term or a restoring term to a tracer climatology 113 113 (when present, see \S~\ref{TRA_dmp}), a forward time differencing scheme 114 114 is used : … … 117 117 \end{equation} 118 118 119 This is diffusive in time and conditionally stable. For example, the119 This is diffusive in time and conditionally stable. The 120 120 conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies_Bk04}: 121 121 \begin{equation} \label{Eq_STP_euler_stability} … … 134 134 135 135 For the vertical diffusion terms, a forward time differencing scheme can be 136 used, but usually the numerical stability condition imp lies a strong136 used, but usually the numerical stability condition imposes a strong 137 137 constraint on the time step. Two solutions are available in \NEMO to overcome 138 138 the stability constraint: $(a)$ a forward time differencing scheme using a 139 139 time splitting technique (\np{ln\_zdfexp}=.true.) or $(b)$ a backward (or implicit) 140 time differencing scheme by\np{ln\_zdfexp}=.false.). In $(a)$, the master140 time differencing scheme (\np{ln\_zdfexp}=.false.). In $(a)$, the master 141 141 time step $\Delta $t is cut into $N$ fractional time steps so that the 142 stability criterion is reduced by a factor of $N$. The computation is doneas142 stability criterion is reduced by a factor of $N$. The computation is performed as 143 143 follows: 144 144 \begin{equation} \label{Eq_STP_ts} … … 152 152 \end{equation} 153 153 with DF a vertical diffusion term. The number of fractional time steps, $N$, is given 154 by setting \np{n \_zdfexp}, (namelist parameter). The scheme $(b)$ is unconditionally154 by setting \np{nn\_zdfexp}, (namelist parameter). The scheme $(b)$ is unconditionally 155 155 stable but diffusive. It can be written as follows: 156 156 \begin{equation} \label{Eq_STP_imp} … … 159 159 160 160 This scheme is rather time consuming since it requires a matrix inversion, 161 but it becomes attractive since a splitting factor of 3 or more is needed162 forthe forward time differencing scheme. For example, the finite difference161 but it becomes attractive since a value of 3 or more is needed for N in 162 the forward time differencing scheme. For example, the finite difference 163 163 approximation of the temperature equation is: 164 164 \begin{equation} \label{Eq_STP_imp_zdf} … … 168 168 \end{equation} 169 169 where RHS is the right hand side of the equation except for the vertical diffusion term. 170 \sgacomment{why change from T to u in the following equation?} 170 171 We rewrite \eqref{Eq_STP_imp} as: 171 172 \begin{equation} \label{Eq_STP_imp_mat} … … 179 180 \end{align*} 180 181 181 \eqref{Eq_STP_imp_mat} is a linear system of equations w hichassociated182 matrix is tridiagonal. Moreover, $c(k)$ and $d(k)$ are positive and the diagonal182 \eqref{Eq_STP_imp_mat} is a linear system of equations with an associated 183 matrix which is tridiagonal. Moreover, $c(k)$ and $d(k)$ are positive and the diagonal 183 184 term is greater than the sum of the two extra-diagonal terms, therefore a special 184 185 adaptation of the Gauss elimination procedure is used to find the solution … … 198 199 \includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_TimeStepping_flowchart.pdf} 199 200 \caption{Sketch of the leapfrog time stepping sequence in \NEMO from \citet{Leclair_Madec_OM09}. 200 The use of semi-implicit computation of the hydrostatic pressure gradient requires201 the tracer equation to be step forward prior to the momentum equation.202 The need of theknowledge of the vertical scale factor (here denoted as $h$)203 requires the sea surface height and the continuity equation to be step forward201 The use of a semi-implicit computation of the hydrostatic pressure gradient requires 202 the tracer equation to be stepped forward prior to the momentum equation. 203 The need for knowledge of the vertical scale factor (here denoted as $h$) 204 requires the sea surface height and the continuity equation to be stepped forward 204 205 prior to the computation of the tracer equation. 205 Note that how the surface pressure gradient $\nabla p_s$ is evaluated is not represented here206 Note that the method for the evaluation of the surface pressure gradient $\nabla p_s$ is not presented here 206 207 (see \S~\ref{DYN_spg}). } 207 208 \end{center} \end{figure} … … 212 213 by introducing a semi-implicit computation of the hydrostatic pressure gradient term 213 214 \citep{Brown_Campana_MWR78}. Instead of evaluating the pressure at $t$, a linear 214 combination of $t-\rdt$, $t$ and $t+\rdt$ is used (see \S~\ref{DYN_hpg_imp}).215 This technique, controlled by \np{nn\_dynhpg\_rst} namelist parameter, does not216 require a significant additional computational coast when tracers and thus density217 is time stepped before the dynamics , a choice madein \NEMO215 combination of values at $t-\rdt$, $t$ and $t+\rdt$ is used (see \S~\ref{DYN_hpg_imp}). 216 This technique, controlled by the \np{nn\_dynhpg\_rst} namelist parameter, does not 217 introduce a significant additional computational cost when tracers and thus density 218 is time stepped before the dynamics. This time step ordering is used in \NEMO 218 219 (Fig.\ref{Fig_TimeStep_flowchart}). 219 220 … … 225 226 as the Forward-Backward scheme used in MOM \citep{Griffies_al_OS05} and more 226 227 efficient than the LF-AM3 scheme (leapfrog time stepping combined with a third order 227 Adams-Moult hon interpolation for the predictor phase) used in ROMS228 Adams-Moulton interpolation for the predictor phase) used in ROMS 228 229 \citep{Shchepetkin_McWilliams_OM05}. 229 230 230 In fact, this technique is efficient when the physical phenomen umthat231 controls the time-step is internal gravity waves (IGWs). Indeed, it is231 In fact, this technique is efficient when the physical phenomenon that 232 limits the time-step is internal gravity waves (IGWs). Indeed, it is 232 233 equivalent to applying a time filter to the pressure gradient to eliminate high 233 234 frequency IGWs. Obviously, the doubling of the time-step is achievable only 234 235 if no other factors control the time-step, such as the stability limits associated 235 236 with advection, diffusion or Coriolis terms. For example, it is useless in low resolution 236 global ocean configurations as inertial oscillations invicinity of the North Pole237 are the limiting factor ofthe time step. It is also often useless in very high238 resolution configurations where thestrong currents and small grid cells exert237 global ocean configurations, since inertial oscillations in the vicinity of the North Pole 238 are the limiting factor for the time step. It is also often useless in very high 239 resolution configurations where strong currents and small grid cells exert 239 240 the strongest constraint on the time step. 241 \sgacomment{ not sure "useless" is the right word here. "valueless", "inefficient"?} 240 242 241 243 … … 247 249 248 250 Significant changes have been introduced by \cite{Leclair_Madec_OM09} in the 249 LF-RA scheme in order to ensure t he tracer conservation and to allow the use of250 a much smaller value of the Asselin filter parameter. The modifications concern251 both the forcing and filtering treatment in the LF-RA scheme.251 LF-RA scheme in order to ensure tracer conservation and to allow the use of 252 a much smaller value of the Asselin filter parameter. The modifications affect 253 both the forcing and filtering treatments in the LF-RA scheme. 252 254 253 255 In a classical LF-RA environment, the forcing term is centred in time, $i.e.$ 254 256 it is time-stepped over a $2\rdt$ period: $x^t = x^t + 2\rdt Q^t $ where $Q$ 255 is the f orcing applied on$x$, and the filter is given by \eqref{Eq_STP_asselin}.257 is the filtered forcing applied to $x$, and the filter is given by \eqref{Eq_STP_asselin}. 256 258 In the modified LF-RA environment, these two formulations have been replaced by: 257 259 \begin{align} … … 261 263 - \gamma\,\rdt \, \left[ Q^{t+\rdt/2} - Q^{t-\rdt/2} \right] \label{Eq_STP_RA} 262 264 \end{align} 265 \sgacomment{Q(t)=f(x(t-dt),x(t),x(t+dt)), Q(t-dt/2)?} 263 266 264 267 The change in the forcing formulation given by \eqref{Eq_STP_forcing} 265 (see Fig.\ref{Fig_MLF_forcing}) has a paramount effect: the forcing term no268 (see Fig.\ref{Fig_MLF_forcing}) has a significant effect: the forcing term no 266 269 longer excites the divergence of odd and even time steps \citep{Leclair_Madec_OM09}. 267 This property allows to improve the LF-RA scheme in two aspects.270 This property improves the LF-RA scheme in two respects. 268 271 First, the LF-RA becomes a truly quasi-second order scheme. Indeed, 269 \eqref{Eq_STP_forcing} used in combination with a careful treatment of thestatic270 instabilit ies(\S\ref{ZDF_evd}) and of the TKE physics (\S\ref{ZDF_tke_ene}),271 the two main other sources of time step divergence, allows a reduction by two order272 \eqref{Eq_STP_forcing} used in combination with a careful treatment of static 273 instability (\S\ref{ZDF_evd}) and of the TKE physics (\S\ref{ZDF_tke_ene}), 274 the two other main sources of time step divergence, allows a reduction by two orders 272 275 of magnitude of the Asselin filter parameter. 273 276 Second, the LF-RA can now ensure the local and global conservation of tracers. 274 Indeed, the time filtering is no more required on the forcing part. The influence of 277 Indeed, time filtering is no longer required on the forcing part. 278 \sgacomment{ but Q is described above as the forcing part!} The influence of 275 279 the forcing in the Asselin filter can be removed by adding a new term in the filter 276 280 (last term in \eqref{Eq_STP_RA} compared to \eqref{Eq_STP_asselin}). Since … … 282 286 in the time filter, \eqref{Eq_STP_RA}, allows an exact evaluation of the 283 287 contribution due to the forcing term between any two time steps, 284 even if separated by only $\rdt$ as the time filter is no more applied onthe288 even if separated by only $\rdt$ since the time filter is no longer applied to the 285 289 forcing term. 286 290 … … 294 298 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 295 299 300 \sgacomment{two methods in caption sound the same} 296 301 297 302 % ------------------------------------------------------------------------------------------------------------- … … 309 314 x^1 = x^0 + \rdt \ \text{RHS}^0 310 315 \end{equation} 311 This is done simply by keeping the leapfrog environment but setting a t the first time step312 only both a half $\rdt$ and the equality of all \textit{before} and \textit{now} fields.316 This is done simply by keeping the leapfrog environment but setting all $x^0$ (\textit{before}) 317 and $x^{1/2}$ (\textit{now}) fields equal at the first time step. 313 318 314 319 It is also possible to restart from a previous computation, by using a … … 322 327 Note that when a semi-implicit scheme is used to evaluate the hydrostatic pressure 323 328 gradient (see \S\ref{DYN_hpg_imp}), an extra three-dimensional field has to be 324 added in the restart file to ensure an exact restartability. This is done onlyoptionally325 via the namelist parameter \np{nn\_dynhpg\_rst}, so that a reduction of the size of326 restart file can be obtained when therestartability is not a key issue (operational327 oceanography or ensemble simulation for seasonal forcast).329 added to the restart file to ensure an exact restartability. This is done optionally 330 via the namelist parameter \np{nn\_dynhpg\_rst}, so that the size of the 331 restart file can be reduced when restartability is not a key issue (operational 332 oceanography or in ensemble simulations for seasonal forecasting). 328 333 329 334 Note the size of the time step used, $\rdt$, is also saved in the restart file. 330 335 When restarting, if the the time step has been changed, a restart using an Euler time 331 stepping is imposed.336 stepping scheme is imposed. 332 337 %%% 333 338 \gmcomment{
Note: See TracChangeset
for help on using the changeset viewer.
