# Changeset 11388 for NEMO/trunk

Ignore:
Timestamp:
2019-08-01T18:09:07+02:00 (21 months ago)
Message:

Add adaptive-implicit vertical advection section to ZDF chapter. May still need a few tweaks next week and refers to an updated traadv_fct.F90 which is not yet submitted. Finals tests are in progress

File:
1 edited

### Legend:

Unmodified
 r11318 \label{subsec:ZDF_aimp} This refers to \citep{shchepetkin_OM15}. TBC ... The adaptive-implicit vertical advection option in NEMO is based on the work of \citep{shchepetkin_OM15}.  In common with most ocean models, the timestep used with NEMO needs to satisfy multiple criteria associated with different physical processes in order to maintain numerical stability. \citep{shchepetkin_OM15} pointed out that the vertical CFL criterion is commonly the most limiting. \citep{lemarie.debreu.ea_OM15} examined the constraints for a range of time and space discretizations and provide the CFL stability criteria for a range of advection schemes. The values for the Leap-Frog with Robert asselin filter time-stepping (as used in NEMO) are reproduced in \autoref{tab:zad_Aimp_CFLcrit}. Treating the vertical advection implicitly can avoid these restrictions but at the cost of large dispersive errors and, possibly, large numerical viscosity. The adaptive-implicit vertical advection option provides a targetted use of the implicit scheme only when and where potential breaches of the vertical CFL condition occur. In many practical applications these events may occur remote from the main area of interest or due to short-lived conditions such that the extra numerical diffusion or viscosity does not greatly affect the overall solution. With such applications, setting: \forcode{ln_zad_Aimp = .true.} should allow much longer model timesteps to be used whilst retaining the accuracy of the high order explicit schemes over most of the domain. \begin{table}[htbp] \begin{center} % \begin{tabular}{cp{70pt}cp{70pt}cp{70pt}cp{70pt}} \begin{tabular}{r|ccc} \hline spatial discretization &   2nd order centered   & 3rd order upwind & 4th order compact  \\ advective CFL criterion     & 0.904 &   0.472  &   0.522    \\ \hline \end{tabular} \caption{ \protect\label{tab:zad_Aimp_CFLcrit} The advective CFL criteria for a range of spatial discretizations for the Leap-Frog with Robert Asselin filter time-stepping ($\nu=0.1$) as given in \citep{lemarie.debreu.ea_OM15}. } \end{center} \end{table} In particular, the advection scheme remains explicit everywhere except where and when local vertical velocities exceed a threshold set just below the explicit stability limit. Once the threshold is reached a tapered transition towards an implicit scheme is used by partitioning the vertical velocity into a part that can be treated explicitly and any excess that must be treated implicitly. The partitioning is achieved via a Courant-number dependent weighting algorithm as described in \citep{shchepetkin_OM15}. The local cell Courant number ($Cu$) used for this partitioning is: \label{eq:Eqn_zad_Aimp_Courant} \begin{split} Cu &= {2 \rdt \over e^n_{3t_{ijk}}} \bigg (\big [ \texttt{Max}(w^n_{ijk},0.0) - \texttt{Min}(w^n_{ijk+1},0.0) \big ]    \\ &\phantom{=} +\big [ \texttt{Max}(e_{{2_u}ij}e^n_{{3_{u}}ijk}u^n_{ijk},0.0) - \texttt{Min}(e_{{2_u}i-1j}e^n_{{3_{u}}i-1jk}u^n_{i-1jk},0.0) \big ] \big / e_{{1_t}ij}e_{{2_t}ij}            \\ &\phantom{=} +\big [ \texttt{Max}(e_{{1_v}ij}e^n_{{3_{v}}ijk}v^n_{ijk},0.0) - \texttt{Min}(e_{{1_v}ij-1}e^n_{{3_{v}}ij-1k}v^n_{ij-1k},0.0) \big ] \big / e_{{1_t}ij}e_{{2_t}ij} \bigg )    \\ \end{split} \noindent and the tapering algorithm follows \citep{shchepetkin_OM15} as: \begin{align} \label{eq:Eqn_zad_Aimp_partition} Cu_{min} &= 0.15 \nonumber \\ Cu_{max} &= 0.27 \nonumber \\ Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\ Fcu    &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\ Cf &= \begin{cases} 0.0                                                        &\text{if $Cu \leq Cu_{min}$} \\ (Cu - Cu_{min})^2 / (Fcu +  (Cu - Cu_{min})^2)             &\text{else if $Cu < Cu_{cut}$} \\ (Cu - Cu_{max}) / Cu                                       &\text{else} \end{cases} \end{align} \noindent With these settings the coefficient ($Cf$) is shown in \autoref{fig:zad_Aimp_coeff} \begin{figure}[!t] \begin{center} \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_coeff} \caption{ \protect\label{fig:zad_Aimp_coeff} The value of the partitioning coefficient ($Cf$) used to partition vertical velocities into parts to be treated implicitly and explicitly for a range of typical Courant numbers (\forcode{ln_zad_Aimp=.true.}) } \end{center} \end{figure} \noindent The partitioning coefficient is used to determine the part of the vertical velocity that must be handled implicitly ($w_i$) and to subtract this from the total vertical velocity ($w_n$) to leave that which can continue to be handled explicitly: \begin{align} \label{eq:Eqn_zad_Aimp_partition2} w_{i_{ijk}} &= Cf_{ijk} w_{n_{ijk}}     \nonumber \\ w_{n_{ijk}} &= (1-Cf_{ijk}) w_{n_{ijk}} \end{align} \noindent Note that the coefficient is such that the treatment is never fully implicit; the three cases from \autoref{eq:Eqn_zad_Aimp_partition} can be considered as: fully-explicit; mixed explicit/implicit and mostly-implicit.  The $w_i$ component is added to the implicit solvers for the vertical mixing in \mdl{dynzdf} and \mdl{trazdf} in a similar way to \citep{shchepetkin_OM15}.  This is sufficient for the flux-limited advection scheme (\forcode{ln_traadv_mus}) but further intervention is required when using the flux-corrected scheme (\forcode{ln_traadv_fct}). For these schemes the implicit upstream fluxes must be added to both the monotonic guess and to the higher order solution when calculating the antidiffusive fluxes. The implicit vertical fluxes are then removed since they are added by the implicit solver later on. The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be used in a wide range of simulations. The following test simulation, however, does illustrate the potential benefits and will hopefully encourage further testing and feedback from users: \subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case} The \href{https://forge.ipsl.jussieu.fr/nemo/chrome/site/doc/NEMO/guide/html/test\_cases.html\#overflow}{OVERFLOW test case} provides a simple illustration of the adaptive-implicit advection in action. The example here differs from the basic test case by only a few extra physics choices namely: \begin{verbatim} ln_dynldf_OFF = .false. ln_dynldf_lap = .true. ln_dynldf_hor = .true. ln_zdfnpc     = .true. ln_traadv_fct = .true. nn_fct_h   =  4 nn_fct_v   =  4 \end{verbatim} \noindent which were chosen to provide a slightly more stable and less noisy solution. The result when using the default value of \forcode{nn_rdt = 10.} without adaptive-implicit vertical velocity is illustrated in \autoref{fig:zad_Aimp_overflow_frames}. The mass of cold water, initially sitting on the shelf, moves down the slope and forms a bottom-trapped, dense plume. Even with these extra physics choices the model is close to stability limits and attempts with \forcode{nn_rdt = 30.} will fail after about 5.5 hours with excessively high horizontal velocities. This time-scale corresponds with the time the plume reaches the steepest part of the topography and, although detected as a horizontal CFL breach, the instability originates from a breach of the vertical CFL limit. This is a good candidate, therefore, for use of the adaptive-implicit vertical advection scheme. The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps are shown in \autoref{fig:zad_Aimp_overflow_all_rdt} (together with the equivalent frames from the base run).  In this simple example the use of the adaptive-implicit vertcal advection scheme has enabled a 12x increase in the model timestep without significantly altering the solution (although at this extreme the plume is more diffuse and has not travelled so far). Notably, the solution with and without the scheme is slightly different even with \forcode{nn_rdt = 10.}; suggesting that the base run was close enough to instability to trigger the scheme despite completing successfully.  To assist in diagnosing how active the scheme is, in both location and time, the 3D implicit and explicit components of the vertical velocity are available via XIOS as \texttt{wimp} and \texttt{wexp} respectively. Likewise, the partitioning coefficient ($Cf$) is also available as \texttt{wi\_cff}. \begin{figure}[!t] \begin{center} \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_overflow_frames} \caption{ \protect\label{fig:zad_Aimp_overflow_frames} A time-series of temperature vertical cross-sections for the OVERFLOW test case. These results are for the default settings with \forcode{nn_rdt=10.0} and without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}). } \end{center} \end{figure} \begin{figure}[!t] \begin{center} \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_overflow_all_rdt} \caption{ \protect\label{fig:zad_Aimp_overflow_all_rdt} Sample temperature vertical cross-sections from mid- and end-run using different values for \forcode{nn_rdt} and with or without adaptive implicit vertical advection. Without the adaptive implicit vertical advection only the run with the shortest timestep is able to run to completion. Note also that the colour-scale has been chosen to confirm that temperatures remain within the original range of 10$^o$ to 20$^o$. } \end{center} \end{figure}