# source:branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/DOC/TexFiles/Chapters/Chap_MISC.tex@5600 Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File size: 35.6 KB
Line
1% ================================================================
2% Chapter � Miscellaneous Topics
3% ================================================================
4\chapter{Miscellaneous Topics}
5\label{MISC}
6\minitoc
7
8\newpage
9$\$\newline    % force a new ligne
10
11% ================================================================
12% Representation of Unresolved Straits
13% ================================================================
14\section{Representation of Unresolved Straits}
15\label{MISC_strait}
16
17In climate modeling, it often occurs that a crucial connections between water masses
18is broken as the grid mesh is too coarse to resolve narrow straits. For example, coarse
19grid spacing typically closes off the Mediterranean from the Atlantic at the Strait of
20Gibraltar. In this case, it is important for climate models to include the effects of salty
21water entering the Atlantic from the Mediterranean. Likewise, it is important for the
22Mediterranean to replenish its supply of water from the Atlantic to balance the net
23evaporation occurring over the Mediterranean region. This problem occurs even in
24eddy permitting simulations. For example, in ORCA 1/4\deg several straits of the Indonesian
25archipelago (Ombai, Lombok...) are much narrow than even a single ocean grid-point.
26
27We describe briefly here the three methods that can be used in \NEMO to handle
28such improperly resolved straits. The first two consist of opening the strait by hand
29while ensuring that the mass exchanges through the strait are not too large by
30either artificially reducing the surface of the strait grid-cells or, locally increasing
31the lateral friction. In the third one, the strait is closed but exchanges of mass,
32heat and salt across the land are allowed.
33Note that such modifications are so specific to a given configuration that no attempt
34has been made to set them in a generic way. However, examples of how
35they can be set up is given in the ORCA 2\deg and 0.5\deg configurations. For example,
36for details of implementation in ORCA2, search:
37\vspace{-10pt}
38\begin{alltt}
39\tiny
40\begin{verbatim}
41IF( cp_cfg == "orca" .AND. jp_cfg == 2 )
42\end{verbatim}
43\end{alltt}
44
45% -------------------------------------------------------------------------------------------------------------
47% -------------------------------------------------------------------------------------------------------------
49\label{MISC_strait_hand}
50
51$\bullet$ reduced scale factor in the cross-strait direction to a value in better agreement
52with the true mean width of the strait. (Fig.~\ref{Fig_MISC_strait_hand}).
53This technique is sometime called "partially open face" or "partially closed cells".
54The key issue here is only to reduce the faces of $T$-cell ($i.e.$ change the value
55of the horizontal scale factors at $u$- or $v$-point) but not the volume of the $T$-cell.
56Indeed, reducing the volume of strait $T$-cell can easily produce a numerical
57instability at that grid point that would require a reduction of the model time step.
58The changes associated with strait management are done in \mdl{domhgr},
59just after the definition or reading of the horizontal scale factors.
60
61$\bullet$ increase of the viscous boundary layer thickness by local increase of the
62fmask value at the coast (Fig.~\ref{Fig_MISC_strait_hand}). This is done in
63\mdl{dommsk} together with the setting of the coastal value of fmask
64(see Section \ref{LBC_coast})
65
66%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
67\begin{figure}[!tbp]     \begin{center}
68\includegraphics[width=0.80\textwidth]{./TexFiles/Figures/Fig_Gibraltar.pdf}
69\includegraphics[width=0.80\textwidth]{./TexFiles/Figures/Fig_Gibraltar2.pdf}
70\caption{   \label{Fig_MISC_strait_hand}
71Example of the Gibraltar strait defined in a $1\deg \times 1\deg$ mesh.
72\textit{Top}: using partially open cells. The meridional scale factor at $v$-point
73is reduced on both sides of the strait to account for the real width of the strait
74(about 20 km). Note that the scale factors of the strait $T$-point remains unchanged.
75\textit{Bottom}: using viscous boundary layers. The four fmask parameters
76along the strait coastlines are set to a value larger than 4, $i.e.$ "strong" no-slip
77case (see Fig.\ref{Fig_LBC_shlat}) creating a large viscous boundary layer
78that allows a reduced transport through the strait.}
79\end{center}   \end{figure}
80%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
81
82% -------------------------------------------------------------------------------------------------------------
84% -------------------------------------------------------------------------------------------------------------
86\label{MISC_strait_cla}
87%--------------------------------------------namcla--------------------------------------------------------
88\namdisplay{namcla}
89%--------------------------------------------------------------------------------------------------------------
90
91\colorbox{yellow}{Add a short description of CLA staff here or in lateral boundary condition chapter?}
92Options are defined through the  \ngn{namcla} namelist variables.
93
94%The problem is resolved here by allowing the mixing of tracers and mass/volume between non-adjacent water columns at nominated regions within the model. Momentum is not mixed. The scheme conserves total tracer content, and total volume (the latter in $z*$- or $s*$-coordinate), and maintains compatibility between the tracer and mass/volume budgets.
95
96% ================================================================
97% Closed seas
98% ================================================================
99\section{Closed seas (\mdl{closea})}
100\label{MISC_closea}
101
102\colorbox{yellow}{Add here a short description of the way closed seas are managed}
103
104
105% ================================================================
106% Sub-Domain Functionality (\textit{nizoom, njzoom}, namelist parameters)
107% ================================================================
108\section{Sub-Domain Functionality (\np{jpizoom}, \np{jpjzoom})}
109\label{MISC_zoom}
110
111The sub-domain functionality, also improperly called the zoom option
112(improperly because it is not associated with a change in model resolution)
113is a quite simple function that allows a simulation over a sub-domain of an
114already defined configuration ($i.e.$ without defining a new mesh, initial
115state and forcings). This option can be useful for testing the user settings
116of surface boundary conditions, or the initial ocean state of a huge ocean
117model configuration while having a small computer memory requirement.
118It can also be used to easily test specific physics in a sub-domain (for example,
119see \citep{Madec_al_JPO96} for a test of the coupling used in the global ocean
120version of OPA between sea-ice and ocean model over the Arctic or Antarctic
121ocean, using a sub-domain). In the standard model, this option does not
122include any specific treatment for the ocean boundaries of the sub-domain:
123they are considered as artificial vertical walls. Nevertheless, it is quite easy
124to add a restoring term toward a climatology in the vicinity of such boundaries
125(see \S\ref{TRA_dmp}).
126
127In order to easily define a sub-domain over which the computation can be
128performed, the dimension of all input arrays (ocean mesh, bathymetry,
129forcing, initial state, ...) are defined as \np{jpidta}, \np{jpjdta} and \np{jpkdta}
130( in \ngn{namcfg} namelist), while the computational domain is defined through
131\np{jpiglo}, \np{jpjglo} and \jp{jpk} (\ngn{namcfg} namelist). When running the
132model over the whole domain, the user sets \np{jpiglo}=\np{jpidta} \np{jpjglo}=\np{jpjdta}
133and \jp{jpk}=\jp{jpkdta}. When running the model over a sub-domain, the user
134has to provide the size of the sub-domain, (\np{jpiglo}, \np{jpjglo}, \np{jpkglo}),
135and the indices of the south western corner as \np{jpizoom} and \np{jpjzoom} in
136the  \ngn{namcfg} namelist (Fig.~\ref{Fig_LBC_zoom}).
137
138Note that a third set of dimensions exist, \jp{jpi}, \jp{jpj} and \jp{jpk} which is
139actually used to perform the computation. It is set by default to \jp{jpi}=\np{jpjglo}
140and \jp{jpj}=\np{jpjglo}, except for massively parallel computing where the
141computational domain is laid out on local processor memories following a 2D
142horizontal splitting. % (see {\S}IV.2-c) ref to the section to be updated
143
144\subsection{Simple subsetting of input files via netCDF attributes}
145
146The extended grids for use with the under-shelf ice cavities will result in redundant rows
147around Antarctica if the ice cavities are not active. A simple mechanism for subsetting
148input files associated with the extended domains has been implemented to avoid the need to
149maintain different sets of input fields for use with or without active ice cavities. The
150existing 'zoom' options are overly complex for this task and marked for deletion anyway.
151This alternative subsetting operates for the j-direction only and works by optionally
152looking for and using a global file attribute (named: \np{open\_ocean\_jstart}) to
153determine the starting j-row for input. The use of this option is best explained with an
154example: Consider an ORCA1 configuration using the extended grid bathymetry and coordinate
155files:
156\vspace{-10pt}
157\begin{alltt}
158\tiny
159\begin{verbatim}
160eORCA1_bathymetry_v2.nc
161eORCA1_coordinates.nc
162\end{verbatim}
163\end{alltt}
164\noindent These files define a horizontal domain of 362x332. Assuming the first row with
165open ocean wet points in the non-isf bathymetry for this set is row 42 (Fortran indexing)
166then the formally correct setting for \np{open\_ocean\_jstart} is 41. Using this value as the
167first row to be read will result in a 362x292 domain which is the same size as the original
168ORCA1 domain. Thus the extended coordinates and bathymetry files can be used with all the
169original input files for ORCA1 if the ice cavities are not active (\np{ln\_isfcav =
170.false.}). Full instructions for achieving this are:
171
172\noindent Add the new attribute to any input files requiring a j-row offset, i.e:
173\vspace{-10pt}
174\begin{alltt}
175\tiny
176\begin{verbatim}
177ncatted  -a open_ocean_jstart,global,a,d,41 eORCA1_coordinates.nc
178ncatted  -a open_ocean_jstart,global,a,d,41 eORCA1_bathymetry_v2.nc
179\end{verbatim}
180\end{alltt}
181
182\noindent Add the logical switch to \ngn{namcfg} in the configuration namelist and set true:
183%--------------------------------------------namcfg--------------------------------------------------------
184\namdisplay{namcfg_orca1}
185%--------------------------------------------------------------------------------------------------------------
186
187\noindent Note the j-size of the global domain is the (extended j-size minus
188\np{open\_ocean\_jstart} + 1 ) and this must match the size of all datasets other than
189bathymetry and coordinates currently. However the option can be extended to any global, 2D
190and 3D, netcdf, input field by adding the:
191\vspace{-10pt}
192\begin{alltt}
193\tiny
194\begin{verbatim}
195lrowattr=ln_use_jattr
196\end{verbatim}
197\end{alltt}
198optional argument to the appropriate \np{iom\_get} call and the \np{open\_ocean\_jstart} attribute to the corresponding input files. It remains the users responsibility to set \np{jpjdta} and \np{jpjglo} values in the \np{namelist\_cfg} file according to their needs.
199
200%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
201\begin{figure}[!ht]    \begin{center}
202\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_LBC_zoom.pdf}
203\caption{   \label{Fig_LBC_zoom}
204Position of a model domain compared to the data input domain when the zoom functionality is used.}
205\end{center}   \end{figure}
206%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
207
208
209% ================================================================
210% Accelerating the Convergence
211% ================================================================
212\section{Accelerating the Convergence (\np{nn\_acc} = 1)}
213\label{MISC_acc}
214%--------------------------------------------namdom-------------------------------------------------------
215\namdisplay{namdom}
216%--------------------------------------------------------------------------------------------------------------
217
218Searching an equilibrium state with an global ocean model requires a very long time
219integration period (a few thousand years for a global model). Due to the size of
220the time step required for numerical stability (less than a few hours),
221this usually requires a large elapsed time. In order to overcome this problem,
222\citet{Bryan1984} introduces a technique that is intended to accelerate
223the spin up to equilibrium. It uses a larger time step in
224the tracer evolution equations than in the momentum evolution
225equations. It does not affect the equilibrium solution but modifies the
226trajectory to reach it.
227
228Options are defined through the  \ngn{namdom} namelist variables.
229The acceleration of convergence option is used when \np{nn\_acc}=1. In that case,
230$\rdt=rn\_rdt$ is the time step of dynamics while $\widetilde{\rdt}=rdttra$ is the
231tracer time-step. the former is set from the \np{rn\_rdt} namelist parameter while the latter
232is computed using a hyperbolic tangent profile and the following namelist parameters :
233\np{rn\_rdtmin}, \np{rn\_rdtmax} and \np{rn\_rdth}. Those three parameters correspond
234to the surface value the deep ocean value and the depth at which the transition occurs, respectively.
235The set of prognostic equations to solve becomes:
236\begin{equation} \label{Eq_acc}
237\begin{split}
238\frac{\partial \textbf{U}_h }{\partial t}
239   &\equiv \frac{\textbf{U}_h ^{t+1}-\textbf{U}_h^{t-1} }{2\rdt} = \ldots \\
240\frac{\partial T}{\partial t} &\equiv \frac{T^{t+1}-T^{t-1}}{2 \widetilde{\rdt}} = \ldots \\
241\frac{\partial S}{\partial t} &\equiv \frac{S^{t+1} -S^{t-1}}{2 \widetilde{\rdt}} = \ldots \\
242\end{split}
243\end{equation}
244
245\citet{Bryan1984} has examined the consequences of this distorted physics.
246Free waves have a slower phase speed, their meridional structure is slightly
247modified, and the growth rate of baroclinically unstable waves is reduced
248but with a wider range of instability. This technique is efficient for
249searching for an equilibrium state in coarse resolution models. However its
250application is not suitable for many oceanic problems: it cannot be used for
251transient or time evolving problems (in particular, it is very questionable
252to use this technique when there is a seasonal cycle in the forcing fields),
253and it cannot be used in high-resolution models where baroclinically
254unstable processes are important. Moreover, the vertical variation of
255$\widetilde{ \rdt}$ implies that the heat and salt contents are no longer
256conserved due to the vertical coupling of the ocean level through both
257advection and diffusion. Therefore \np{rn\_rdtmin} = \np{rn\_rdtmax} should be
258a more clever choice.
259
260
261% ================================================================
262% Accuracy and Reproducibility
263% ================================================================
264\section{Accuracy and Reproducibility (\mdl{lib\_fortran})}
265\label{MISC_fortran}
266
267\subsection{Issues with intrinsinc SIGN function (\key{nosignedzero})}
268\label{MISC_sign}
269
270The SIGN(A, B) is the \textsc {Fortran} intrinsic function delivers the magnitude
271of A with the sign of B. For example, SIGN(-3.0,2.0) has the value 3.0.
272The problematic case is when the second argument is zero, because, on platforms
273that support IEEE arithmetic, zero is actually a signed number.
274There is a positive zero and a negative zero.
275
276In \textsc{Fortran}~90, the processor was required always to deliver a positive result for SIGN(A, B)
277if B was zero. Nevertheless, in \textsc{Fortran}~95, the processor is allowed to do the correct thing
278and deliver ABS(A) when B is a positive zero and -ABS(A) when B is a negative zero.
279This change in the specification becomes apparent only when B is of type real, and is zero,
280and the processor is capable of distinguishing between positive and negative zero,
281and B is negative real zero. Then SIGN delivers a negative result where, under \textsc{Fortran}~90
282rules,  it used to return a positive result.
283This change may be especially sensitive for the ice model, so we overwrite the intrinsinc
284function with our own function simply performing :   \\
285\verb?   IF( B >= 0.e0 ) THEN   ;   SIGN(A,B) = ABS(A)  ?    \\
286\verb?   ELSE                   ;   SIGN(A,B) =-ABS(A)     ?  \\
287\verb?   ENDIF    ? \\
288This feature can be found in \mdl{lib\_fortran} module and is effective when \key{nosignedzero}
289is defined. We use a CPP key as the overwritting of a intrinsic function can present
290performance issues with some computers/compilers.
291
292
293\subsection{MPP reproducibility}
294\label{MISC_glosum}
295
296The numerical reproducibility of simulations on distributed memory parallel computers
297is a critical issue. In particular, within NEMO global summation of distributed arrays
298is most susceptible to rounding errors, and their propagation and accumulation cause
299uncertainty in final simulation reproducibility on different numbers of processors.
300To avoid so, based on \citet{He_Ding_JSC01} review of different technics,
301we use a so called self-compensated summation method. The idea is to estimate
302the roundoff error, store it in a buffer, and then add it back in the next addition.
303
304Suppose we need to calculate $b = a_1 + a_2 + a_3$. The following algorithm
305will allow to split the sum in two ($sum_1 = a_{1} + a_{2}$ and $b = sum_2 = sum_1 + a_3$)
306with exactly the same rounding errors as the sum performed all at once.
307\begin{align*}
308   sum_1 \ \  &= a_1 + a_2 \\
309   error_1     &= a_2 + ( a_1 - sum_1 ) \\
310   sum_2 \ \  &= sum_1 + a_3 + error_1 \\
311   error_2     &= a_3 + error_1 + ( sum_1 - sum_2 ) \\
312   b \qquad \ &= sum_2 \\
313\end{align*}
314This feature can be found in \mdl{lib\_fortran} module and is effective when \key{mpp\_rep}.
315In that case, all calls to glob\_sum function (summation over the entire basin excluding
316duplicated rows and columns due to cyclic or north fold boundary condition as well as
317overlap MPP areas).
318Note this implementation may be sensitive to the optimization level.
319
320\subsection{MPP scalability}
321\label{MISC_mppsca}
322
323The default method of communicating values across the north-fold in distributed memory applications
324(\key{mpp\_mpi}) uses a \textsc{MPI\_ALLGATHER} function to exchange values from each processing
325region in the northern row with every other processing region in the northern row. This enables a
326global width array containing the top 4 rows to be collated on every northern row processor and then
327folded with a simple algorithm. Although conceptually simple, this "All to All" communication will
328hamper performance scalability for large numbers of northern row processors. From version 3.4
329onwards an alternative method is available which only performs direct "Peer to Peer" communications
330between each processor and its immediate "neighbours" across the fold line. This is achieved by
331using the default \textsc{MPI\_ALLGATHER} method during initialisation to help identify the "active"
332neighbours. Stored lists of these neighbours are then used in all subsequent north-fold exchanges to
333restrict exchanges to those between associated regions. The collated global width array for each
334region is thus only partially filled but is guaranteed to be set at all the locations actually
335required by each individual for the fold operation. This alternative method should give identical
336results to the default \textsc{ALLGATHER} method and is recommended for large values of \np{jpni}.
337The new method is activated by setting \np{ln\_nnogather} to be true ({\bf nammpp}). The
338reproducibility of results using the two methods should be confirmed for each new, non-reference
339configuration.
340
341% ================================================================
342% Model optimisation, Control Print and Benchmark
343% ================================================================
344\section{Model Optimisation, Control Print and Benchmark}
345\label{MISC_opt}
346%--------------------------------------------namctl-------------------------------------------------------
347\namdisplay{namctl}
348%--------------------------------------------------------------------------------------------------------------
349
350 \gmcomment{why not make these bullets into subsections?}
351Options are defined through the  \ngn{namctl} namelist variables.
352
353$\bullet$ Vector optimisation:
354
355\key{vectopt\_loop} enables the internal loops to collapse. This is very
356a very efficient way to increase the length of vector calculations and thus
357to speed up the model on vector computers.
358
359% Add here also one word on NPROMA technique that has been found useless, since compiler have made significant progress during the last decade.
360
361% Add also one word on NEC specific optimisation (Novercheck option for example)
362
363$\bullet$ Control print %: describe here 4 things:
364
3651- \np{ln\_ctl} : compute and print the trends averaged over the interior domain
366in all TRA, DYN, LDF and ZDF modules. This option is very helpful when
367diagnosing the origin of an undesired change in model results.
368
3692- also \np{ln\_ctl} but using the nictl and njctl namelist parameters to check
370the source of differences between mono and multi processor runs.
371
3723- \key{esopa} (to be rename key\_nemo) : which is another option for model
373management. When defined, this key forces the activation of all options and
374CPP keys. For example, all tracer and momentum advection schemes are called!
375Therefore the model results have no physical meaning.
376However, this option forces both the compiler and the model to run through
377all the \textsc{Fortran} lines of the model. This allows the user to check for obvious
378compilation or execution errors with all CPP options, and errors in namelist options.
379
3804- last digit comparison (\np{nn\_bit\_cmp}). In an MPP simulation, the computation of
381a sum over the whole domain is performed as the summation over all processors of
382each of their sums over their interior domains. This double sum never gives exactly
383the same result as a single sum over the whole domain, due to truncation differences.
384The "bit comparison" option has been introduced in order to be able to check that
385mono-processor and multi-processor runs give exactly the same results.
386%THIS is to be updated with the mpp_sum_glo  introduced in v3.3
387% nn_bit_cmp  today only check that the nn_cla = 0 (no cross land advection)
388
389$\bullet$  Benchmark (\np{nn\_bench}). This option defines a benchmark run based on
390a GYRE configuration (see \S\ref{CFG_gyre}) in which the resolution remains the same
391whatever the domain size. This allows a very large model domain to be used, just by
392changing the domain size (\jp{jpiglo}, \jp{jpjglo}) and without adjusting either the time-step
393or the physical parameterisations.
394
395
396% ================================================================
397% Elliptic solvers (SOL)
398% ================================================================
399\section{Elliptic solvers (SOL)}
400\label{MISC_sol}
401%--------------------------------------------namdom-------------------------------------------------------
402\namdisplay{namsol}
403%--------------------------------------------------------------------------------------------------------------
404
405When the filtered sea surface height option is used, the surface pressure gradient is
406computed in \mdl{dynspg\_flt}. The force added in the momentum equation is solved implicitely.
407It is thus solution of an elliptic equation \eqref{Eq_PE_flt} for which two solvers are available:
408a Successive-Over-Relaxation scheme (SOR) and a preconditioned conjugate gradient
410the value of \np{nn\_solv}   \ngn{namsol} namelist variable.
411
412The PCG is a very efficient method for solving elliptic equations on vector computers.
413It is a fast and rather easy method to use; which are attractive features for a large
414number of ocean situations (variable bottom topography, complex coastal geometry,
415variable grid spacing, open or cyclic boundaries, etc ...). It does not require
416a search for an optimal parameter as in the SOR method. However, the SOR has
417been retained because it is a linear solver, which is a very useful property when
418using the adjoint model of \NEMO.
419
420At each time step, the time derivative of the sea surface height at time step $t+1$
421(or equivalently the divergence of the \textit{after} barotropic transport) that appears
422in the filtering forced is the solution of the elliptic equation obtained from the horizontal
423divergence of the vertical summation of \eqref{Eq_PE_flt}.
424Introducing the following coefficients:
425\begin{equation}  \label{Eq_sol_matrix}
426\begin{aligned}
427&c_{i,j}^{NS}  &&= {2 \rdt }^2 \; \frac{H_v (i,j) \; e_{1v} (i,j)}{e_{2v}(i,j)}              \\
428&c_{i,j}^{EW} &&= {2 \rdt }^2 \; \frac{H_u (i,j) \; e_{2u} (i,j)}{e_{1u}(i,j)}            \\
429&b_{i,j} &&= \delta_i \left[ e_{2u}M_u \right] - \delta_j \left[ e_{1v}M_v \right]\ ,   \\
430\end{aligned}
431\end{equation}
432the resulting five-point finite difference equation is given by:
433\begin{equation}  \label{Eq_solmat}
434\begin{split}
435       c_{i+1,j}^{NS} D_{i+1,j}  + \;  c_{i,j+1}^{EW} D_{i,j+1}
436  +   c_{i,j}    ^{NS} D_{i-1,j}   + \;  c_{i,j}    ^{EW} D_{i,j-1}                                          &    \\
437  -    \left(c_{i+1,j}^{NS} + c_{i,j+1}^{EW} + c_{i,j}^{NS} + c_{i,j}^{EW} \right)   D_{i,j}  &=  b_{i,j}
438\end{split}
439\end{equation}
440\eqref{Eq_solmat} is a linear symmetric system of equations. All the elements of
441the corresponding matrix \textbf{A} vanish except those of five diagonals. With
442the natural ordering of the grid points (i.e. from west to east and from
443south to north), the structure of \textbf{A} is block-tridiagonal with
444tridiagonal or diagonal blocks. \textbf{A} is a positive-definite symmetric
445matrix of size $(jpi \cdot jpj)^2$, and \textbf{B}, the right hand side of
446\eqref{Eq_solmat}, is a vector.
447
448Note that in the linear free surface case, the depth that appears in \eqref{Eq_sol_matrix}
449does not vary with time, and thus the matrix can be computed once for all. In non-linear free surface
450(\key{vvl} defined) the matrix have to be updated at each time step.
451
452% -------------------------------------------------------------------------------------------------------------
453%       Successive Over Relaxation
454% -------------------------------------------------------------------------------------------------------------
455\subsection{Successive Over Relaxation (\np{nn\_solv}=2, \mdl{solsor})}
456\label{MISC_solsor}
457
458Let us introduce the four cardinal coefficients:
459\begin{align*}
460a_{i,j}^S &= c_{i,j    }^{NS}/d_{i,j}     &\qquad  a_{i,j}^W &= c_{i,j}^{EW}/d_{i,j}       \\
461a_{i,j}^E &= c_{i,j+1}^{EW}/d_{i,j}    &\qquad   a_{i,j}^N &= c_{i+1,j}^{NS}/d_{i,j}
462\end{align*}
463where $d_{i,j} = c_{i,j}^{NS}+ c_{i+1,j}^{NS} + c_{i,j}^{EW} + c_{i,j+1}^{EW}$
464(i.e. the diagonal of the matrix). \eqref{Eq_solmat} can be rewritten as:
465\begin{equation}  \label{Eq_solmat_p}
466\begin{split}
467a_{i,j}^{N}  D_{i+1,j} +\,a_{i,j}^{E}  D_{i,j+1} +\, a_{i,j}^{S}  D_{i-1,j} +\,a_{i,j}^{W} D_{i,j-1}  -  D_{i,j} = \tilde{b}_{i,j}
468\end{split}
469\end{equation}
470with $\tilde b_{i,j} = b_{i,j}/d_{i,j}$. \eqref{Eq_solmat_p} is the equation actually solved
471with the SOR method. This method used is an iterative one. Its algorithm can be
472summarised as follows (see \citet{Haltiner1980} for a further discussion):
473
474initialisation (evaluate a first guess from previous time step computations)
475\begin{equation}
476D_{i,j}^0 = 2 \, D_{i,j}^t - D_{i,j}^{t-1}
477\end{equation}
478iteration $n$, from $n=0$ until convergence, do :
479\begin{equation} \label{Eq_sor_algo}
480\begin{split}
481R_{i,j}^n  = &a_{i,j}^{N} D_{i+1,j}^n       +\,a_{i,j}^{E}  D_{i,j+1} ^n
482         +\, a_{i,j}^{S}  D_{i-1,j} ^{n+1}+\,a_{i,j}^{W} D_{i,j-1} ^{n+1}
483                 -  D_{i,j}^n - \tilde{b}_{i,j}                                           \\
484D_{i,j} ^{n+1}  = &D_{i,j} ^{n}   + \omega \;R_{i,j}^n
485\end{split}
486\end{equation}
487where \textit{$\omega$ }satisfies $1\leq \omega \leq 2$. An optimal value exists for
488\textit{$\omega$} which significantly accelerates the convergence, but it has to be
489adjusted empirically for each model domain (except for a uniform grid where an
490analytical expression for \textit{$\omega$} can be found \citep{Richtmyer1967}).
491The value of $\omega$ is set using \np{rn\_sor}, a \textbf{namelist} parameter.
492The convergence test is of the form:
493\begin{equation}
494\delta = \frac{\sum\limits_{i,j}{R_{i,j}^n}{R_{i,j}^n}}
495                    {\sum\limits_{i,j}{ \tilde{b}_{i,j}^n}{\tilde{b}_{i,j}^n}} \leq \epsilon
496\end{equation}
497where $\epsilon$ is the absolute precision that is required. It is recommended
498that a value smaller or equal to $10^{-6}$ is used for $\epsilon$ since larger
499values may lead to numerically induced basin scale barotropic oscillations.
500The precision is specified by setting \np{rn\_eps} (\textbf{namelist} parameter).
501In addition, two other tests are used to halt the iterative algorithm. They involve
502the number of iterations and the modulus of the right hand side. If the former
503exceeds a specified value, \np{nn\_max} (\textbf{namelist} parameter),
504or the latter is greater than $10^{15}$, the whole model computation is stopped
505and the last computed time step fields are saved in a abort.nc NetCDF file.
506In both cases, this usually indicates that there is something wrong in the model
507configuration (an error in the mesh, the initial state, the input forcing,
508or the magnitude of the time step or of the mixing coefficients). A typical value of
509$nn\_max$ is a few hundred when $\epsilon = 10^{-6}$, increasing to a few
510thousand when $\epsilon = 10^{-12}$.
511The vectorization of the SOR algorithm is not straightforward. The scheme
512contains two linear recurrences on $i$ and $j$. This inhibits the vectorisation.
513\eqref{Eq_sor_algo} can be been rewritten as:
514\begin{equation}
515\begin{split}
516R_{i,j}^n
517= &a_{i,j}^{N}  D_{i+1,j}^n +\,a_{i,j}^{E}  D_{i,j+1} ^n
518 +\,a_{i,j}^{S}  D_{i-1,j} ^{n}+\,_{i,j}^{W} D_{i,j-1} ^{n} -  D_{i,j}^n - \tilde{b}_{i,j}      \\
519R_{i,j}^n = &R_{i,j}^n - \omega \;a_{i,j}^{S}\; R_{i,j-1}^n                                             \\
520R_{i,j}^n = &R_{i,j}^n - \omega \;a_{i,j}^{W}\; R_{i-1,j}^n
521\end{split}
522\end{equation}
523This technique slightly increases the number of iteration required to reach the convergence,
524but this is largely compensated by the gain obtained by the suppression of the recurrences.
525
526Another technique have been chosen, the so-called red-black SOR. It consist in solving successively
527\eqref{Eq_sor_algo} for odd and even grid points. It also slightly reduced the convergence rate
528but allows the vectorisation. In addition, and this is the reason why it has been chosen, it is able to handle the north fold boundary condition used in ORCA configuration ($i.e.$ tri-polar global ocean mesh).
529
530The SOR method is very flexible and can be used under a wide range of conditions,
531including irregular boundaries, interior boundary points, etc. Proofs of convergence, etc.
532may be found in the standard numerical methods texts for partial differential equations.
533
534% -------------------------------------------------------------------------------------------------------------
536% -------------------------------------------------------------------------------------------------------------
537\subsection{Preconditioned Conjugate Gradient  (\np{nn\_solv}=1, \mdl{solpcg}) }
538\label{MISC_solpcg}
539
540\textbf{A} is a definite positive symmetric matrix, thus solving the linear
541system \eqref{Eq_solmat} is equivalent to the minimisation of a quadratic
542functional:
543\begin{equation*}
544\textbf{Ax} = \textbf{b} \leftrightarrow \textbf{x} =\text{inf}_{y} \,\phi (\textbf{y})
546\phi (\textbf{y}) = 1/2 \langle \textbf{Ay},\textbf{y}\rangle - \langle \textbf{b},\textbf{y} \rangle
547\end{equation*}
548where $\langle , \rangle$ is the canonical dot product. The idea of the
549conjugate gradient method is to search for the solution in the following
550iterative way: assuming that $\textbf{x}^n$ has been obtained, $\textbf{x}^{n+1}$
551is found from $\textbf {x}^{n+1}={\textbf {x}}^n+\alpha^n{\textbf {d}}^n$ which satisfies:
552\begin{equation*}
553{\textbf{ x}}^{n+1}=\text{inf} _{{\textbf{ y}}={\textbf{ x}}^n+\alpha^n \,{\textbf{ d}}^n} \,\phi ({\textbf{ y}})\;\;\Leftrightarrow \;\;\frac{d\phi }{d\alpha}=0
554\end{equation*}
555and expressing $\phi (\textbf{y})$ as a function of \textit{$\alpha$}, we obtain the
556value that minimises the functional:
557\begin{equation*}
558\alpha ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle  / \langle {\textbf{ A d}^n, \textbf{d}^n} \rangle
559\end{equation*}
560where $\textbf{r}^n = \textbf{b}-\textbf{A x}^n = \textbf{A} (\textbf{x}-\textbf{x}^n)$
561is the error at rank $n$. The descent vector $\textbf{d}^n$ s chosen to be dependent
562on the error: $\textbf{d}^n = \textbf{r}^n + \beta^n \,\textbf{d}^{n-1}$. $\beta ^n$
563is searched such that the descent vectors form an orthogonal basis for the dot
564product linked to \textbf{A}. Expressing the condition
565$\langle \textbf{A d}^n, \textbf{d}^{n-1} \rangle = 0$ the value of $\beta ^n$ is found:
566 $\beta ^n = \langle{ \textbf{r}^n , \textbf{r}^n} \rangle / \langle {\textbf{r}^{n-1}, \textbf{r}^{n-1}} \rangle$.
567 As a result, the errors $\textbf{r}^n$ form an orthogonal
568base for the canonic dot product while the descent vectors $\textbf{d}^n$ form
569an orthogonal base for the dot product linked to \textbf{A}. The resulting
570algorithm is thus the following one:
571
572initialisation :
573\begin{equation*}
574\begin{split}
575\textbf{x}^0 &= D_{i,j}^0   = 2 D_{i,j}^t - D_{i,j}^{t-1}       \quad, \text{the initial guess }     \\
576\textbf{r}^0 &= \textbf{d}^0 = \textbf{b} - \textbf{A x}^0       \\
577\gamma_0 &= \langle{ \textbf{r}^0 , \textbf{r}^0} \rangle
578\end{split}
579\end{equation*}
580
581iteration $n,$ from $n=0$ until convergence, do :
582\begin{equation}
583\begin{split}
584\text{z}^n& = \textbf{A d}^n \\
585\alpha_n &= \gamma_n /  \langle{ \textbf{z}^n , \textbf{d}^n} \rangle \\
586\textbf{x}^{n+1} &= \textbf{x}^n + \alpha_n \,\textbf{d}^n \\
587\textbf{r}^{n+1} &= \textbf{r}^n - \alpha_n \,\textbf{z}^n \\
588\gamma_{n+1} &= \langle{ \textbf{r}^{n+1} , \textbf{r}^{n+1}} \rangle \\
589\beta_{n+1} &= \gamma_{n+1}/\gamma_{n}  \\
590\textbf{d}^{n+1} &= \textbf{r}^{n+1} + \beta_{n+1}\; \textbf{d}^{n}\\
591\end{split}
592\end{equation}
593
594
595The convergence test is:
596\begin{equation}
597\delta = \gamma_{n}\; / \langle{ \textbf{b} , \textbf{b}} \rangle \leq \epsilon
598\end{equation}
599where $\epsilon$ is the absolute precision that is required. As for the SOR algorithm,
600the whole model computation is stopped when the number of iterations, \np{nn\_max}, or
601the modulus of the right hand side of the convergence equation exceeds a
602specified value (see \S\ref{MISC_solsor} for a further discussion). The required
603precision and the maximum number of iterations allowed are specified by setting
604\np{rn\_eps} and \np{nn\_max} (\textbf{namelist} parameters).
605
606It can be demonstrated that the above algorithm is optimal, provides the exact
607solution in a number of iterations equal to the size of the matrix, and that
608the convergence rate is faster as the matrix is closer to the identity matrix,
609$i.e.$ its eigenvalues are closer to 1. Therefore, it is more efficient to solve
610a better conditioned system which has the same solution. For that purpose,
611we introduce a preconditioning matrix \textbf{Q} which is an approximation
612of \textbf{A} but much easier to invert than \textbf{A}, and solve the system:
613\begin{equation} \label{Eq_pmat}
614\textbf{Q}^{-1} \textbf{A x} = \textbf{Q}^{-1} \textbf{b}
615\end{equation}
616
617The same algorithm can be used to solve \eqref{Eq_pmat} if instead of the
618canonical dot product the following one is used:
619${\langle{ \textbf{a} , \textbf{b}} \rangle}_Q = \langle{ \textbf{a} , \textbf{Q b}} \rangle$, and
620if $\textbf{\~{b}} = \textbf{Q}^{-1}\;\textbf{b}$ and $\textbf{\~{A}} = \textbf{Q}^{-1}\;\textbf{A}$
621are substituted to \textbf{b} and \textbf{A} \citep{Madec_al_OM88}.
622In \NEMO, \textbf{Q} is chosen as the diagonal of \textbf{ A}, i.e. the simplest form for
623\textbf{Q} so that it can be easily inverted. In this case, the discrete formulation of
624\eqref{Eq_pmat} is in fact given by \eqref{Eq_solmat_p} and thus the matrix and
625right hand side are computed independently from the solver used.
626
627% ================================================================
628
629
630
631
632
633
634
635
636
637
638
639
Note: See TracBrowser for help on using the repository browser.