doc/algorithm/crank-nicol.rst

4f2617d475 Jeff*0001 .. _crank-nicolson_baro:
                0002
                0003 Crank-Nicolson barotropic time stepping
                0004 ---------------------------------------
                0005
                0006 The full implicit time stepping described previously is
                0007 unconditionally stable but damps the fast gravity waves, resulting in
                0008 a loss of potential energy. The modification presented now allows one
7c73d7406b timo*0009 to combine an implicit part (:math:`\gamma,\beta`) and an explicit
                0010 part (:math:`1-\gamma,1-\beta`) for the surface pressure gradient
                0011 (:math:`\gamma`) and for the barotropic flow divergence
                0012 (:math:`\beta`). For instance, :math:`\gamma=\beta=1` is the previous fully implicit
                0013 scheme; :math:`\gamma=\beta=1/2` is the non damping (energy
4f2617d475 Jeff*0014 conserving), unconditionally stable, Crank-Nicolson scheme;
7c73d7406b timo*0015 :math:`(\gamma,\beta)=(1,0)` or :math:`=(0,1)` corresponds to the
4f2617d475 Jeff*0016 forward - backward scheme that conserves energy but is only stable for
7c73d7406b timo*0017 small time steps. In the code, :math:`\gamma,\beta` are defined as parameters,
4f2617d475 Jeff*0018 respectively :varlink:`implicSurfPress`, :varlink:`implicDiv2DFlow`. They are read
                0019 from the main parameter file ``data`` (namelist ``PARM01``) and are set
                0020 by default to 1,1.
                0021
                ** Warning **Wide character in print at /usr/local/share/lxr/source line 1030, <$git> line 23.


0022 Equations :eq:`ustar-backward-free-surface` –
                0023 :eq:`vn+1-backward-free-surface` are modified as follows:
                0024
                0025 .. math::
                0026    \frac{ \vec{\bf v}^{n+1} }{ \Delta t }
7c73d7406b timo*0027    + {\bf \nabla}_h b_s [ \gamma {\eta}^{n+1} + (1-\gamma) {\eta}^{n} ]
0bad585a21 Navi*0028    + \epsilon_{\rm nh} {\bf \nabla}_h {\phi'_{\rm nh}}^{n+1}
4f2617d475 Jeff*0029    = \frac{ \vec{\bf v}^{n} }{ \Delta t }
                0030    + \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
0bad585a21 Navi*0031    + {\bf \nabla}_h {\phi'_{\rm hyd}}^{(n+1/2)}
4f2617d475 Jeff*0032
                0033 .. math::
0bad585a21 Navi*0034    \epsilon_{\rm fs} \frac{ {\eta}^{n+1} - {\eta}^{n} }{ \Delta t}
                0035    + {\bf \nabla}_h \cdot \int_{R_{\rm fixed}}^{R_o}
7c73d7406b timo*0036    [ \beta \vec{\bf v}^{n+1} + (1-\beta) \vec{\bf v}^{n}] dr
0bad585a21 Navi*0037    = \epsilon_{\rm fw} ({\mathcal{P-E}})
4f2617d475 Jeff*0038    :label: eta-n+1-CrankNick
                0039
                0040 We set
                0041
                0042 .. math::
                0043      \begin{aligned}
                0044      \vec{\bf v}^* & = &
                0045      \vec{\bf v} ^{n} + \Delta t \vec{\bf G}_{\vec{\bf v}} ^{(n+1/2)}
7c73d7406b timo*0046      + (\gamma-1) \Delta t {\bf \nabla}_h b_s {\eta}^{n}
0bad585a21 Navi*0047      + \Delta t {\bf \nabla}_h {\phi'_{\rm hyd}}^{(n+1/2)}
4f2617d475 Jeff*0048      \\
                0049      {\eta}^* & = &
0bad585a21 Navi*0050      \epsilon_{fs} {\eta}^{n} + \epsilon_{\rm fw} \Delta t ({\mathcal{P-E}})
                0051      - \Delta t {\bf \nabla}_h \cdot \int_{R_{\rm fixed}}^{R_o}
7c73d7406b timo*0052      [ \beta \vec{\bf v}^* + (1-\beta) \vec{\bf v}^{n}] dr\end{aligned}
4f2617d475 Jeff*0053
0bad585a21 Navi*0054 In the hydrostatic case :math:`\epsilon_{\rm nh}=0`, allowing us to find
4f2617d475 Jeff*0055 :math:`{\eta}^{n+1}`, thus:
                0056
                0057 .. math::
0bad585a21 Navi*0058      \epsilon_{\rm fs} {\eta}^{n+1} -
                0059      {\bf \nabla}_h \cdot \gamma\beta \Delta t^2 b_s (R_o - R_{\rm fixed})
4f2617d475 Jeff*0060      {\bf \nabla}_h {\eta}^{n+1}
                0061      = {\eta}^*
                0062
                0063 and then to compute (:filelink:`CORRECTION_STEP <model/src/correction_step.F>`):
                0064
                0065 .. math::
                0066      \vec{\bf v}^{n+1} = \vec{\bf v}^{*}
7c73d7406b timo*0067      - \gamma \Delta t {\bf \nabla}_h b_s {\eta}^{n+1}
4f2617d475 Jeff*0068
                0069 Notes:
                0070
                0071 #. The RHS term of equation :eq:`eta-n+1-CrankNick` corresponds the
94151a9b18 Jeff*0072    contribution of fresh water flux ({\mathcal{P-E}}) to the free-surface variations
0bad585a21 Navi*0073    (:math:`\epsilon_{fw}=1`, :varlink:`useRealFreshWaterFlux` ``=.TRUE.`` in parameter
4f2617d475 Jeff*0074    file ``data``). In order to remain consistent with the tracer equation,
                0075    specially in the non-linear free-surface formulation, this term is
                0076    also affected by the Crank-Nicolson time stepping. The RHS reads:
0bad585a21 Navi*0077    :math:`\epsilon_{\rm fw} ( \beta ({\mathcal{P-E}})^{n+1/2} + (1-\beta) ({\mathcal{P-E}})^{n-1/2} )`
4f2617d475 Jeff*0078
                0079
                0080 #. The stability criteria with Crank-Nicolson time stepping for the pure
                0081    linear gravity wave problem in cartesian coordinates is:
                0082
7c73d7406b timo*0083    -  :math:`\gamma + \beta < 1` : unstable
4f2617d475 Jeff*0084
7c73d7406b timo*0085    -  :math:`\gamma \geq 1/2` and :math:`\beta \geq 1/2` : stable
4f2617d475 Jeff*0086
0bad585a21 Navi*0087    -  :math:`\gamma + \beta \geq 1` : stable if :math:`c_{\rm max}^2 (\gamma - 1/2)(\beta - 1/2) + 1 \geq 0`
                0088       with :math:`c_{\rm max} = 2 \Delta t \sqrt{gH} \sqrt{ \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} }`
4f2617d475 Jeff*0089
                0090
                0091 #. A similar mixed forward/backward time-stepping is also available for
0bad585a21 Navi*0092    the non-hydrostatic algorithm, with a fraction :math:`\gamma_{\rm nh}`
                0093    (:math:`0 < \gamma_{\rm nh} \leq 1`) of the non-hydrostatic pressure
4f2617d475 Jeff*0094    gradient being evaluated at time step :math:`n+1` (backward in time)
0bad585a21 Navi*0095    and the remaining part (:math:`1 - \gamma_{\rm nh}`) being evaluated at
4f2617d475 Jeff*0096    time step :math:`n` (forward in time). The run-time parameter
                0097    :varlink:`implicitNHPress` corresponding to the implicit fraction
0bad585a21 Navi*0098    :math:`\gamma_{\rm nh}` of the non-hydrostatic pressure is set by default
7c73d7406b timo*0099    to the implicit fraction :math:`\gamma` of surface pressure
4f2617d475 Jeff*0100    (:varlink:`implicSurfPress`), but can also be specified independently (in
                0101    main parameter file ``data``, namelist ``PARM01``).