Changeset 6289 for trunk/DOC/TexFiles/Chapters/Chap_MISC.tex
- Timestamp:
- 2016-02-05T00:47:05+01:00 (8 years ago)
- File:
-
- 1 edited
-
trunk/DOC/TexFiles/Chapters/Chap_MISC.tex (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/DOC/TexFiles/Chapters/Chap_MISC.tex
r5118 r6289 1 1 % ================================================================ 2 % Chapter �Miscellaneous Topics2 % Chapter ——— Miscellaneous Topics 3 3 % ================================================================ 4 4 \chapter{Miscellaneous Topics} … … 34 34 has been made to set them in a generic way. However, examples of how 35 35 they can be set up is given in the ORCA 2\deg and 0.5\deg configurations. For example, 36 for details of implementation in ORCA2, search: 37 \vspace{-10pt} 38 \begin{alltt} 39 \tiny 40 \begin{verbatim} 41 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) 42 \end{verbatim} 43 \end{alltt} 36 for details of implementation in ORCA2, search: 37 \texttt{ IF( cp\_cfg == "orca" .AND. jp\_cfg == 2 ) } 44 38 45 39 % ------------------------------------------------------------------------------------------------------------- … … 80 74 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 81 75 82 % -------------------------------------------------------------------------------------------------------------83 % Cross Land Advection84 % -------------------------------------------------------------------------------------------------------------85 \subsection{Cross Land Advection (\mdl{tracla})}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?}92 Options 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 76 96 77 % ================================================================ … … 208 189 209 190 % ================================================================ 210 % Accelerating the Convergence211 % ================================================================212 \section{Accelerating the Convergence (\np{nn\_acc} = 1)}213 \label{MISC_acc}214 %--------------------------------------------namdom-------------------------------------------------------215 \namdisplay{namdom}216 %--------------------------------------------------------------------------------------------------------------217 218 Searching an equilibrium state with an global ocean model requires a very long time219 integration period (a few thousand years for a global model). Due to the size of220 the time step required for numerical stability (less than a few hours),221 this usually requires a large elapsed time. In order to overcome this problem,222 \citet{Bryan1984} introduces a technique that is intended to accelerate223 the spin up to equilibrium. It uses a larger time step in224 the tracer evolution equations than in the momentum evolution225 equations. It does not affect the equilibrium solution but modifies the226 trajectory to reach it.227 228 Options are defined through the \ngn{namdom} namelist variables.229 The 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 the231 tracer time-step. the former is set from the \np{rn\_rdt} namelist parameter while the latter232 is 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 correspond234 to the surface value the deep ocean value and the depth at which the transition occurs, respectively.235 The 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.246 Free waves have a slower phase speed, their meridional structure is slightly247 modified, and the growth rate of baroclinically unstable waves is reduced248 but with a wider range of instability. This technique is efficient for249 searching for an equilibrium state in coarse resolution models. However its250 application is not suitable for many oceanic problems: it cannot be used for251 transient or time evolving problems (in particular, it is very questionable252 to use this technique when there is a seasonal cycle in the forcing fields),253 and it cannot be used in high-resolution models where baroclinically254 unstable processes are important. Moreover, the vertical variation of255 $\widetilde{ \rdt}$ implies that the heat and salt contents are no longer256 conserved due to the vertical coupling of the ocean level through both257 advection and diffusion. Therefore \np{rn\_rdtmin} = \np{rn\_rdtmax} should be258 a more clever choice.259 260 261 % ================================================================262 191 % Accuracy and Reproducibility 263 192 % ================================================================ … … 370 299 the source of differences between mono and multi processor runs. 371 300 372 3- \key{esopa} (to be rename key\_nemo) : which is another option for model 373 management. When defined, this key forces the activation of all options and 374 CPP keys. For example, all tracer and momentum advection schemes are called! 375 Therefore the model results have no physical meaning. 376 However, this option forces both the compiler and the model to run through 377 all the \textsc{Fortran} lines of the model. This allows the user to check for obvious 378 compilation or execution errors with all CPP options, and errors in namelist options. 379 380 4- last digit comparison (\np{nn\_bit\_cmp}). In an MPP simulation, the computation of 301 %%gm to be removed both here and in the code 302 3- last digit comparison (\np{nn\_bit\_cmp}). In an MPP simulation, the computation of 381 303 a sum over the whole domain is performed as the summation over all processors of 382 304 each of their sums over their interior domains. This double sum never gives exactly … … 386 308 %THIS is to be updated with the mpp_sum_glo introduced in v3.3 387 309 % nn_bit_cmp today only check that the nn_cla = 0 (no cross land advection) 310 %%gm end 388 311 389 312 $\bullet$ Benchmark (\np{nn\_bench}). This option defines a benchmark run based on … … 393 316 or the physical parameterisations. 394 317 395 396 % ================================================================ 397 % Elliptic solvers (SOL) 398 % ================================================================ 399 \section{Elliptic solvers (SOL)} 400 \label{MISC_sol} 401 %--------------------------------------------namdom------------------------------------------------------- 402 \namdisplay{namsol} 403 %-------------------------------------------------------------------------------------------------------------- 404 405 When the filtered sea surface height option is used, the surface pressure gradient is 406 computed in \mdl{dynspg\_flt}. The force added in the momentum equation is solved implicitely. 407 It is thus solution of an elliptic equation \eqref{Eq_PE_flt} for which two solvers are available: 408 a Successive-Over-Relaxation scheme (SOR) and a preconditioned conjugate gradient 409 scheme(PCG) \citep{Madec_al_OM88, Madec_PhD90}. The solver is selected trough the 410 the value of \np{nn\_solv} \ngn{namsol} namelist variable. 411 412 The PCG is a very efficient method for solving elliptic equations on vector computers. 413 It is a fast and rather easy method to use; which are attractive features for a large 414 number of ocean situations (variable bottom topography, complex coastal geometry, 415 variable grid spacing, open or cyclic boundaries, etc ...). It does not require 416 a search for an optimal parameter as in the SOR method. However, the SOR has 417 been retained because it is a linear solver, which is a very useful property when 418 using the adjoint model of \NEMO. 419 420 At 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 422 in the filtering forced is the solution of the elliptic equation obtained from the horizontal 423 divergence of the vertical summation of \eqref{Eq_PE_flt}. 424 Introducing 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} 432 the 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 441 the corresponding matrix \textbf{A} vanish except those of five diagonals. With 442 the natural ordering of the grid points (i.e. from west to east and from 443 south to north), the structure of \textbf{A} is block-tridiagonal with 444 tridiagonal or diagonal blocks. \textbf{A} is a positive-definite symmetric 445 matrix of size $(jpi \cdot jpj)^2$, and \textbf{B}, the right hand side of 446 \eqref{Eq_solmat}, is a vector. 447 448 Note that in the linear free surface case, the depth that appears in \eqref{Eq_sol_matrix} 449 does 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 458 Let us introduce the four cardinal coefficients: 459 \begin{align*} 460 a_{i,j}^S &= c_{i,j }^{NS}/d_{i,j} &\qquad a_{i,j}^W &= c_{i,j}^{EW}/d_{i,j} \\ 461 a_{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*} 463 where $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} 467 a_{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} 470 with $\tilde b_{i,j} = b_{i,j}/d_{i,j}$. \eqref{Eq_solmat_p} is the equation actually solved 471 with the SOR method. This method used is an iterative one. Its algorithm can be 472 summarised as follows (see \citet{Haltiner1980} for a further discussion): 473 474 initialisation (evaluate a first guess from previous time step computations) 475 \begin{equation} 476 D_{i,j}^0 = 2 \, D_{i,j}^t - D_{i,j}^{t-1} 477 \end{equation} 478 iteration $n$, from $n=0$ until convergence, do : 479 \begin{equation} \label{Eq_sor_algo} 480 \begin{split} 481 R_{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} \\ 484 D_{i,j} ^{n+1} = &D_{i,j} ^{n} + \omega \;R_{i,j}^n 485 \end{split} 486 \end{equation} 487 where \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 489 adjusted empirically for each model domain (except for a uniform grid where an 490 analytical expression for \textit{$\omega$} can be found \citep{Richtmyer1967}). 491 The value of $\omega$ is set using \np{rn\_sor}, a \textbf{namelist} parameter. 492 The 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} 497 where $\epsilon$ is the absolute precision that is required. It is recommended 498 that a value smaller or equal to $10^{-6}$ is used for $\epsilon$ since larger 499 values may lead to numerically induced basin scale barotropic oscillations. 500 The precision is specified by setting \np{rn\_eps} (\textbf{namelist} parameter). 501 In addition, two other tests are used to halt the iterative algorithm. They involve 502 the number of iterations and the modulus of the right hand side. If the former 503 exceeds a specified value, \np{nn\_max} (\textbf{namelist} parameter), 504 or the latter is greater than $10^{15}$, the whole model computation is stopped 505 and the last computed time step fields are saved in a abort.nc NetCDF file. 506 In both cases, this usually indicates that there is something wrong in the model 507 configuration (an error in the mesh, the initial state, the input forcing, 508 or 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 510 thousand when $\epsilon = 10^{-12}$. 511 The vectorization of the SOR algorithm is not straightforward. The scheme 512 contains 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} 516 R_{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} \\ 519 R_{i,j}^n = &R_{i,j}^n - \omega \;a_{i,j}^{S}\; R_{i,j-1}^n \\ 520 R_{i,j}^n = &R_{i,j}^n - \omega \;a_{i,j}^{W}\; R_{i-1,j}^n 521 \end{split} 522 \end{equation} 523 This technique slightly increases the number of iteration required to reach the convergence, 524 but this is largely compensated by the gain obtained by the suppression of the recurrences. 525 526 Another 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 528 but 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 530 The SOR method is very flexible and can be used under a wide range of conditions, 531 including irregular boundaries, interior boundary points, etc. Proofs of convergence, etc. 532 may be found in the standard numerical methods texts for partial differential equations. 533 534 % ------------------------------------------------------------------------------------------------------------- 535 % Preconditioned Conjugate Gradient 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 541 system \eqref{Eq_solmat} is equivalent to the minimisation of a quadratic 542 functional: 543 \begin{equation*} 544 \textbf{Ax} = \textbf{b} \leftrightarrow \textbf{x} =\text{inf}_{y} \,\phi (\textbf{y}) 545 \quad , \qquad 546 \phi (\textbf{y}) = 1/2 \langle \textbf{Ay},\textbf{y}\rangle - \langle \textbf{b},\textbf{y} \rangle 547 \end{equation*} 548 where $\langle , \rangle$ is the canonical dot product. The idea of the 549 conjugate gradient method is to search for the solution in the following 550 iterative way: assuming that $\textbf{x}^n$ has been obtained, $\textbf{x}^{n+1}$ 551 is 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*} 555 and expressing $\phi (\textbf{y})$ as a function of \textit{$\alpha $}, we obtain the 556 value 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*} 560 where $\textbf{r}^n = \textbf{b}-\textbf{A x}^n = \textbf{A} (\textbf{x}-\textbf{x}^n)$ 561 is the error at rank $n$. The descent vector $\textbf{d}^n$ s chosen to be dependent 562 on the error: $\textbf{d}^n = \textbf{r}^n + \beta^n \,\textbf{d}^{n-1}$. $\beta ^n$ 563 is searched such that the descent vectors form an orthogonal basis for the dot 564 product 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 568 base for the canonic dot product while the descent vectors $\textbf{d}^n$ form 569 an orthogonal base for the dot product linked to \textbf{A}. The resulting 570 algorithm is thus the following one: 571 572 initialisation : 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 581 iteration $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 595 The convergence test is: 596 \begin{equation} 597 \delta = \gamma_{n}\; / \langle{ \textbf{b} , \textbf{b}} \rangle \leq \epsilon 598 \end{equation} 599 where $\epsilon $ is the absolute precision that is required. As for the SOR algorithm, 600 the whole model computation is stopped when the number of iterations, \np{nn\_max}, or 601 the modulus of the right hand side of the convergence equation exceeds a 602 specified value (see \S\ref{MISC_solsor} for a further discussion). The required 603 precision and the maximum number of iterations allowed are specified by setting 604 \np{rn\_eps} and \np{nn\_max} (\textbf{namelist} parameters). 605 606 It can be demonstrated that the above algorithm is optimal, provides the exact 607 solution in a number of iterations equal to the size of the matrix, and that 608 the 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 610 a better conditioned system which has the same solution. For that purpose, 611 we introduce a preconditioning matrix \textbf{Q} which is an approximation 612 of \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 617 The same algorithm can be used to solve \eqref{Eq_pmat} if instead of the 618 canonical dot product the following one is used: 619 ${\langle{ \textbf{a} , \textbf{b}} \rangle}_Q = \langle{ \textbf{a} , \textbf{Q b}} \rangle$, and 620 if $\textbf{\~{b}} = \textbf{Q}^{-1}\;\textbf{b}$ and $\textbf{\~{A}} = \textbf{Q}^{-1}\;\textbf{A}$ 621 are substituted to \textbf{b} and \textbf{A} \citep{Madec_al_OM88}. 622 In \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 625 right hand side are computed independently from the solver used. 626 627 % ================================================================ 628 629 630 631 632 633 634 635 636 637 638 639 318 % ================================================================ 319 320 321 322 323
Note: See TracChangeset
for help on using the changeset viewer.
