View source for Time-stepping Schemes Review

<div class="title">Time-stepping Schemes Review</div>
Numerical time stepping uses a discrete approximation to: 

<math display="block"> \frac{\partial \phi(t)}{\partial t} = {\cal F}(t) </math><!--\eqno{(1)}-->

where <math>\phi</math> represents one of <math>u</math>, <math>v</math>, <math>C</math>, or <math>\zeta</math>
and <math>{\cal F}(t)</math>
represents all the right-hand-side terms. In ROMS, the goal is to find time-stepping schemes which are accurate where they are valid and damping on unresolved signals ([[Bibliography#ShchepetkinAF_2005a | Shchepetkin and McWilliams (2005)]] and [[Bibliography#ShchepetkinAF_2008b | Shchepetkin and McWilliams (2009)]]). Also, the preference is for time-stepping schemes requiring only one set of the right-hand-side terms so that different time-stepping schemes can be used for different terms in the equations. Finally, as mentioned in Table [[Numerical_Solution_Technique#Table_timestep1 | Time Step]], not all versions of ROMS use the same time-stepping algorithm. We list some time-stepping schemes here which are used or have been used by the ROMS/SCRUM family of models, plus a few to help explain some of the more esoteric ones.

== Euler ==
The simplest approximation is the Euler time step:

<math display="block"> \frac{\phi(t + \Delta t) - \phi(t)}{\Delta t} = {\cal F}(t) </math><!--\eqno{(2)}-->

where you predict the next <math>\phi</math> value based only on the current fields.  This method is accurate to first order in <math>\Delta t</math>; however, it is unconditionally unstable with respect to advection.

== Leapfrog ==
The leapfrog time step is accurate to O(<math>\Delta t^2</math>):

<math display="block"> \frac{\phi(t + \Delta t) - \phi(t - \Delta t)}{2\Delta t} = {\cal F}(t). </math><!--\eqno{(3)}-->

This time step is more accurate, but it is unconditionally unstable with respect to diffusion.  Also, the even and odd time steps tend to diverge in a computational mode. This computational mode can be damped by taking correction steps.  SCRUM's time step on the depth-integrated equations was a leapfrog step with a trapezoidal correction (LF-TR) on every step, which uses a leapfrog step to obtain an initial guess of <math>\phi(t+\Delta t)</math>.  We will call the right-hand-side terms calculated from this initial guess <math>{\cal F}^*(t+\Delta t)</math>:

<math display="block"> \frac{\phi(t + \Delta t) - \phi(t)}{\Delta t} = \frac{1}{2}
  \left[ {\cal F}(t) + {\cal F}^*(t+\Delta t) \right]. </math><!--\eqno{(4)}-->
  
This leapfrog-trapezoidal time step is stable with respect to diffusion and it strongly damps the computational mode.  However, the right-hand-side terms are computed twice per time step.

== Third-order Adams-Bashforth (AB3) ==
The time step on SCRUM's full 3-D fields is done with a third-order Adams-Bashforth step.  It uses three time-levels of the right-hand-side terms:

<math display="block"> \frac{\phi(t + \Delta t) - \phi(t)}{\Delta t} =
  \alpha {\cal F}(t) +
  \beta  {\cal F}(t - \Delta t) +
  \gamma {\cal F}(t - 2 \Delta t) </math><!--\eqno{(5)}-->
  
where the coefficients <math>\alpha</math>, <math>\beta</math> and <math>\gamma</math> are chosen to obtain a third-order estimate of <math>\phi(t + \Delta t)</math>.  We use a Taylor series expansion:

<math display="block"> \frac{\phi(t + \Delta t) - \phi(t)}{\Delta t} =
   \phi^{\prime} + \frac{\Delta t}{2} \phi^{\prime \prime} +
   \frac{\Delta t^2}{6} \phi^{\prime \prime \prime} + \cdots )</math><!--\eqno{(6)}-->
   
where

<math display="block"> \begin{align}{\cal F}(t) & = \phi^{\prime} \\
   {\cal F}(t - \Delta t) & = \phi^{\prime}
   - \Delta t \phi^{\prime \prime}
   + \frac{\Delta t^2}{2} \phi^{\prime \prime \prime} + \cdots \\
   {\cal F}(t - 2\Delta t) & = \phi^{\prime}
   - 2\Delta t \phi^{\prime \prime}
   + 2\Delta t^2 \phi^{\prime \prime \prime} + \cdots \end{align}</math>
   
We find that the coefficients are:

<math display="block">
   \alpha = \frac{23}{12}, \qquad
   \beta  = - \frac{4}{3}, \qquad
   \gamma = \frac{5}{12} 
</math>

This requires one time level for the physical fields and three time levels of the right-hand-side information and requires special treatment on startup.

== Forward-Backward ==
In the discrete approximation equation <!--replace "discrete approximation equation" with equation (1)--> above, we assume that multiple equations for any number of variables are time stepped synchronously. For coupled equations, we can actually do better by time stepping asynchronously. Consider these equations:

<math display="block"> \begin{align}\frac{\partial \zeta}{\partial t} &= {\cal F}(u) \\
   \frac{\partial u}{\partial t} &= {\cal G}(\zeta) \end{align}</math><!--\eqno{(7)}-->

If we time step them alternately, we can always be using the newest information:

<math display="block"> \begin{align}\zeta^{n+1} &= \zeta^n + {\cal F}(u^n) \Delta t \\
   u^{n+1} &= u^n + {\cal G}(\zeta^{n+1}) \Delta t \end{align}</math><!--\eqno{(8)}-->
   
This scheme is second-order accurate and is stable for longer time steps than many other schemes. It is however unstable for advection terms.

== Forward-Backward Feedback (RK2-FB) ==
One option for solving the coupled discrete approximation equations <!--replace "coupled discrete approximation equations" with equation (7)--> is a predictor-corrector with predictor step:

<math display="block"> \begin{align}\zeta^{n+1,\star} &= \zeta^n + {\cal F}(u^n)\Delta t \\
   u^{n+1,\star} &= u^n + \left[\beta {\cal G}(\zeta^{n+1,\star}) +
   (1-\beta) {\cal G}(\zeta^n)\right] \Delta t \end{align}</math><!--\eqno{(9)}-->

and corrector step:

<math display="block"> \begin{align}\zeta^{n+1} &= \zeta^n + \frac{1}{2} \left[{\cal F}(u^{n+1,\star}) +
   {\cal F}(u^n) \right] \Delta t \\
   u^{n+1} &= u^n + \frac{1}{2} \left[\epsilon {\cal G}(\zeta^{n+1}) +
   (1-\epsilon){\cal G}(\zeta^{n+1,\star}) + {\cal G}(\zeta^n)
   \right] \Delta t \end{align}</math><!--\eqno{(9)}-->

Setting <math>\beta = \epsilon = 0</math> in the above, it becomes a standard second order Runge-Kutta scheme, which is unstable for a non-dissipative system. Adding the <math>\beta</math> and <math>\epsilon</math> terms adds Forward-Backward feedback to this algorithm, and allows us to improve both its accuracy and stability. The choice of <math>\beta = 1/3</math> and <math>\epsilon = 2/3</math> leads to a stable third-order scheme with <math>\alpha_{\max} = 2.14093</math> .

== LF-TR and LF-AM3 with FB Feedback==
Another possibility is a leapfrog predictor:

<math display="block"> \begin{align} \zeta^{n+1,\star} &= \zeta^{n-1} + 2{\cal F}(u^n) \Delta t \\
   u^{n+1,\star} &= u^{n-1} + 2\left\{ (1-2\beta){\cal G}(\zeta^n)
   + \beta \left[{\cal G}(\zeta^{n+1,\star}) + {\cal G}(\zeta^{n-1}) \right]
   \right\} \Delta t \end{align}</math><!--\eqno{(11)}-->

and either a two-time trapezoidal or a three-time Adams-Moulton corrector:

<math display="block"> \begin{align} \zeta^{n+1} &= \zeta^n + \left[ \left( \frac{1}{2}-\gamma \right)
   {\cal F}(u^{n+1,\star}) + \left( \frac{1}{2}+ 2\gamma \right) {\cal F}(u^n) -
   \gamma {\cal F}(u^{n-1}) \right] \Delta t \\
   u^{n+1} &= u^n + \left\{ \left( \frac{1}{2}-\gamma \right) \left[
   \epsilon {\cal G}(\zeta^{n+1}) +
   (1-\epsilon){\cal G}(\zeta^{n+1,\star}) \right] +
   \left( \frac{1}{2} + 2\gamma \right) {\cal G}(\zeta^n) - \gamma  {\cal
   G}(\zeta^{n-1}) \right\} \Delta t \end{align}</math><!--\eqno{(12)}-->

where the parameters <math>\beta</math> and <math>\epsilon</math> introduce FB-feedback during both stages, while <math>\gamma</math> controls the type of corrector scheme. Without FB-feedback, we have:

<math display="block">\beta = \epsilon = 0 \Rightarrow \begin{cases}
   \gamma = 0 &\Rightarrow \hbox{LF-TR,} &\alpha_{\max} = \sqrt{2} \\
   \gamma = 1/12 &\Rightarrow \hbox{LF-AM3,} &\alpha_{\max} = 1.5874 \\
   \gamma = 0.0804 &\Rightarrow \hbox{max stability,} &\alpha_{\max} = 1.5876 \end{cases}</math>

Keeping <math>\gamma = 1/12</math> maintains third-order accuracy. The most accurate choices for <math>\beta</math> and <math>\epsilon</math> are <math>\beta = 17/120</math> and <math>\epsilon = 11/20</math>, which yields a fourth-order scheme with low numerical dispersion and diffusion and <math>\alpha_{\max} = 1.851640</math>.

== Generalized FB with an AB3-AM4 Step ==
One drawback of the previous two schemes is the need to evaluate the right-hand-side (r.h.s.) terms twice per time step. The original forward-backward scheme manages to achieve <math>\alpha_{\max} = 2</math> while only evaluating each r.h.s. term once. It is possible to contruct a robust scheme using a combination of a three-time AB3-like step for <math>\zeta</math> with a four-time AM4-like step for <math>u</math>:

<math display="block"> \begin{align}\zeta^{n+1} &= \zeta^n + \left[ \left( \frac{3}{2}+\beta \right)
   {\cal F}(u^n) - \left( \frac{1}{2}+ 2\beta \right) {\cal F}(u^{n-1}) +
   \beta {\cal F}(u^{n-2}) \right] \Delta t \\
   u^{n+1} &= u^n + \left[ \left( \frac{1}{2}+\gamma + 2\epsilon \right)
   {\cal G}(\zeta^{n+1}) +
   \left( \frac{1}{2}-2\gamma-3\epsilon \right){\cal G}(\zeta^n) +
   \gamma {\cal G}(\zeta^{n-1}) + \epsilon  {\cal G}(\zeta^{n-2}) \right]
   \Delta t \end{align}</math><!--\eqno{(13)}-->

This is second-order accurate for any set of values for <math>\beta</math>, <math>\gamma</math>, and <math>\epsilon</math>. It is third-order accurate if <math>\beta = 5/12</math>. However, it has a wider stability range for <math>\beta = 0.281105</math>. [[Bibliography#ShchepetkinAF_2008b | Shchepetkin and McWilliams (2009)]] choose to use this scheme for the barotropic mode in their model with <math>\beta= 0.281105</math>, <math>\gamma = 0.0880</math>, and <math>\epsilon = 0.013</math>, to obtain <math>\alpha_{\max} = 1.7802</math>, close to the value of 2 for pure FB, but with better stability properties for the combination of waves, advection, and Coriolis terms.