/[MITgcm]/manual/s_phys_pkgs/text/seaice.tex
ViewVC logotype

Diff of /manual/s_phys_pkgs/text/seaice.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.19 by mlosch, Wed Dec 14 11:22:42 2011 UTC revision 1.20 by mlosch, Mon Mar 31 11:30:21 2014 UTC
# Line 61  no additional CPP options are required. Line 61  no additional CPP options are required.
61  (see Section \ref{sec:buildingCode}).  (see Section \ref{sec:buildingCode}).
62    
63  Parts of the SEAICE code can be enabled or disabled at compile time  Parts of the SEAICE code can be enabled or disabled at compile time
64  via CPP preprocessor flags. These options are set in either  via CPP preprocessor flags. These options are set in
65  \code{SEAICE\_OPTIONS.h} or in \code{ECCO\_CPPOPTIONS.h}.  \code{SEAICE\_OPTIONS.h}.
66  Table \ref{tab:pkg:seaice:cpp} summarizes these options.  Table \ref{tab:pkg:seaice:cpp} summarizes the most important ones.
67    
68  \begin{table}[!ht]  \begin{table}[!ht]
69  \centering  \centering
# Line 80  Table \ref{tab:pkg:seaice:cpp} summarize Line 80  Table \ref{tab:pkg:seaice:cpp} summarize
80          \code{SEAICE\_CGRID} &          \code{SEAICE\_CGRID} &
81            LSR solver on C-grid (rather than original B-grid) \\            LSR solver on C-grid (rather than original B-grid) \\
82          \code{SEAICE\_ALLOW\_EVP} &          \code{SEAICE\_ALLOW\_EVP} &
83            use EVP rather than LSR rheology solver \\            enable use of EVP rheology solver \\
84            \code{SEAICE\_ALLOW\_JFNK} &
85              enable use of JFNK rheology solver \\
86          \code{SEAICE\_EXTERNAL\_FLUXES} &          \code{SEAICE\_EXTERNAL\_FLUXES} &
87            use EXF-computed fluxes as starting point \\            use EXF-computed fluxes as starting point \\
88          \code{SEAICE\_MULTICATEGORY} &          \code{SEAICE\_ZETA\_SMOOTHREG} &
89            enable 8-category thermodynamics (by default undefined)\\            use differentialable regularization for viscosities \\
90          \code{SEAICE\_VARIABLE\_FREEZING\_POINT} &          \code{SEAICE\_VARIABLE\_FREEZING\_POINT} &
91            enable linear dependence of the freezing point on salinity            enable linear dependence of the freezing point on salinity
92            (by default undefined)\\            (by default undefined)\\
93          \code{ALLOW\_SEAICE\_FLOODING} &          \code{ALLOW\_SEAICE\_FLOODING} &
94            enable snow to ice conversion for submerged sea-ice \\            enable snow to ice conversion for submerged sea-ice \\
95          \code{SEAICE\_SALINITY} &          \code{SEAICE\_VARIABLE\_SALINITY} &
96            enable "salty" sea-ice (by default undefined) \\            enable sea-ice with variable salinity (by default undefined) \\
97          \code{SEAICE\_AGE} &          \code{SEAICE\_SITRACER} &
98            enable "age tracer" sea-ice (by default undefined) \\            enable sea-ice tracer package (by default undefined) \\
         \code{SEAICE\_CAP\_HEFF} &  
           enable capping of sea-ice thickness to MAX\_HEFF \\ \hline  
99          \code{SEAICE\_BICE\_STRESS} &          \code{SEAICE\_BICE\_STRESS} &
100            B-grid only for backward compatiblity: turn on ice-stress on            B-grid only for backward compatiblity: turn on ice-stress on
101            ocean\\            ocean\\
# Line 105  Table \ref{tab:pkg:seaice:cpp} summarize Line 105  Table \ref{tab:pkg:seaice:cpp} summarize
105        \hline        \hline
106      \end{tabular}      \end{tabular}
107    }    }
108    \caption{~}    \caption{Some of the most relevant CPP preprocessor flags in the
109        \code{seaice}-package.}
110  \end{table}  \end{table}
111    
112  %----------------------------------------------------------------------  %----------------------------------------------------------------------
# Line 139  Table~\ref{tab:pkg:seaice:runtimeparms} Line 140  Table~\ref{tab:pkg:seaice:runtimeparms}
140    over grid cell in meters; initializes variable \code{HSNOW};    over grid cell in meters; initializes variable \code{HSNOW};
141  \item[\code{HsaltFile}:] Initial salinity of sea ice averaged over grid  \item[\code{HsaltFile}:] Initial salinity of sea ice averaged over grid
142    cell in g/m$^2$; initializes variable \code{HSALT};    cell in g/m$^2$; initializes variable \code{HSALT};
 \item[\code{IceAgeFile}:] Initial ice age of sea ice averaged over grid  
   cell in seconds; initializes variable \code{ICEAGE};  
143  \end{description}  \end{description}
144    
145  %----------------------------------------------------------------------  %----------------------------------------------------------------------
# Line 187  improved: Line 186  improved:
186  \item the code has been rewritten for an Arakawa C-grid, both B- and  \item the code has been rewritten for an Arakawa C-grid, both B- and
187    C-grid variants are available; the C-grid code allows for no-slip    C-grid variants are available; the C-grid code allows for no-slip
188    and free-slip lateral boundary conditions;    and free-slip lateral boundary conditions;
189  \item two different solution methods for solving the nonlinear  \item three different solution methods for solving the nonlinear
190    momentum equations have been adopted: LSOR \citep{zhang97}, and EVP    momentum equations have been adopted: LSOR \citep{zhang97}, EVP
191    \citep{hun97};    \citep{hun97}, JFNK \citep{lemieux10,losch14:_jfnk};
192  \item ice-ocean stress can be formulated as in \citet{hibler87} or as in  \item ice-ocean stress can be formulated as in \citet{hibler87} or as in
193    \citet{cam08};    \citet{cam08};
194  \item ice variables are advected by sophisticated, conservative  \item ice variables are advected by sophisticated, conservative
# Line 296  air and ocean drag coefficients; $\rho_{ Line 295  air and ocean drag coefficients; $\rho_{
295  densities; and $R_{air/ocean}$ are rotation matrices that act on the  densities; and $R_{air/ocean}$ are rotation matrices that act on the
296  wind/current vectors.  wind/current vectors.
297    
298  \paragraph{Viscous-Plastic (VP) Rheology and LSOR solver \label{sec:pkg:seaice:VPdynamics}}~\\  \paragraph{Viscous-Plastic (VP) Rheology\label{sec:pkg:seaice:VPrheology}}~\\
299  %  %
300  For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can  For an isotropic system the stress tensor $\sigma_{ij}$ ($i,j=1,2$) can
301  be related to the ice strain rate and strength by a nonlinear  be related to the ice strain rate and strength by a nonlinear
# Line 349  recommended). For stress tensor computat Line 348  recommended). For stress tensor computat
348  = 2\,\Delta\zeta$ \citep{hibler95} is used so that the stress state  = 2\,\Delta\zeta$ \citep{hibler95} is used so that the stress state
349  always lies on the elliptic yield curve by definition.  always lies on the elliptic yield curve by definition.
350    
351  In the current implementation, the VP-model is integrated with the  Defining the CPP-flag \code{SEAICE\_ZETA\_SMOOTHREG} in
352  semi-implicit line successive over relaxation (LSOR)-solver of  \code{SEAICE\_OPTIONS.h} before compiling replaces the method for
353  \citet{zhang97}, which allows for long time steps that, in our case,  bounding $\zeta$ by a smooth (differentiable) expression:
354  are limited by the explicit treatment of the Coriolis term. The  \begin{equation}
355  explicit treatment of the Coriolis term does not represent a severe    \label{eq:zetaregsmooth}
356  limitation because it restricts the time step to approximately the    \begin{split}
357  same length as in the ocean model where the Coriolis term is also    \zeta &= \zeta_{\max}\tanh\left(\frac{P}{2\,\min(\Delta,\Delta_{\min})
358  treated explicitly.        \,\zeta_{\max}}\right)\\
359      &= \frac{P}{2\Delta^*}
360      \tanh\left(\frac{\Delta^*}{\min(\Delta,\Delta_{\min})}\right)
361      \end{split}
362    \end{equation}
363    where $\Delta_{\min}=10^{-20}\text{\,s}^{-1}$ is chosen to avoid divisions
364    by zero.
365    
366    \paragraph{LSR and  JFNK solver \label{sec:pkg:seaice:LSRJFNK}}~\\
367    %
368    % By default, the VP-model is integrated by a Pickwith the
369    % semi-implicit line successive over relaxation (LSOR)-solver of
370    % \citet{zhang97}, which allows for long time steps that, in our case,
371    % are limited by the explicit treatment of the Coriolis term. The
372    % explicit treatment of the Coriolis term does not represent a severe
373    % limitation because it restricts the time step to approximately the
374    % same length as in the ocean model where the Coriolis term is also
375    % treated explicitly.
376    
377    \newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}}
378    %
379    In the matrix notation, the discretized momentum equations can be
380    written as
381    \begin{equation}
382      \label{eq:matrixmom}
383      \mat{A}(\vek{x})\,\vek{x} = \vek{b}(\vek{x}).
384    \end{equation}
385    The solution vector $\vek{x}$ consists of the two velocity components
386    $u$ and $v$ that contain the velocity variables at all grid points and
387    at one time level. The standard (and default) method for solving
388    Eq.\,(\ref{eq:matrixmom}) in the sea ice component of the
389    \mbox{MITgcm}, as in many sea ice models, is an iterative Picard
390    solver: in the $k$-th iteration a linearized form
391    $\mat{A}(\vek{x}^{k-1})\,\vek{x}^{k} = \vek{b}(\vek{x}^{k-1})$ is
392    solved (in the case of the MITgcm it is a Line Successive (over)
393    Relaxation (LSR) algorithm \citep{zhang97}).  Picard solvers converge
394    slowly, but generally the iteration is terminated after only a few
395    non-linear steps \citep{zhang97, lemieux09} and the calculation
396    continues with the next time level. This method is the default method
397    in the MITgcm. The number of non-linear iteration steps or pseudo-time
398    steps can be controlled by the runtime parameter
399    \code{NPSEUDOTIMESTEPS} (default is 2).
400    
401    In order to overcome the poor convergence of the Picard-solver,
402    \citet{lemieux10} introduced a Jacobian-free Newton-Krylov solver for
403    the sea ice momentum equations. This solver is also implemented in the
404    MITgcm \citep{losch14:_jfnk}. The Newton method transforms minimizing
405    the residual $\vek{F}(\vek{x}) = \mat{A}(\vek{x})\,\vek{x} -
406    \vek{b}(\vek{x})$ to finding the roots of a multivariate Taylor
407    expansion of the residual \vek{F} around the previous ($k-1$) estimate
408    $\vek{x}^{k-1}$:
409    \begin{equation}
410      \label{eq:jfnktaylor}
411      \vek{F}(\vek{x}^{k-1}+\delta\vek{x}^{k}) =
412      \vek{F}(\vek{x}^{k-1}) + \vek{F}'(\vek{x}^{k-1})\,\delta\vek{x}^{k}
413    \end{equation}
414    with the Jacobian $\mat{J}\equiv\vek{F}'$. The root
415    $\vek{F}(\vek{x}^{k-1}+\delta\vek{x}^{k})=0$ is found by solving
416    \begin{equation}
417      \label{eq:jfnklin}
418      \mat{J}(\vek{x}^{k-1})\,\delta\vek{x}^{k} = -\vek{F}(\vek{x}^{k-1})
419    \end{equation}
420    for $\delta\vek{x}^{k}$. The next ($k$-th) estimate is given by
421    $\vek{x}^{k}=\vek{x}^{k-1}+a\,\delta\vek{x}^{k}$. In order to avoid
422    overshoots the factor $a$ is iteratively reduced in a line search
423    ($a=1, \frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \ldots$) until
424    $\|\vek{F}(\vek{x}^k)\| < \|\vek{F}(\vek{x}^{k-1})\|$, where
425    $\|\cdot\|=\int\cdot\,dx^2$ is the $L_2$-norm. In practice, the line
426    search is stopped at $a=\frac{1}{8}$. The line search starts after
427    $\code{SEAICE\_JFNK\_lsIter}$ non-linear Newton iterations (off by
428    default).
429    
430    
431    Forming the Jacobian $\mat{J}$ explicitly is often avoided as ``too
432    error prone and time consuming'' \citep{knoll04:_jfnk}. Instead,
433    Krylov methods only require the action of \mat{J} on an arbitrary
434    vector \vek{w} and hence allow a matrix free algorithm for solving
435    Eq.\,(\ref{eq:jfnklin}) \citep{knoll04:_jfnk}. The action of \mat{J}
436    can be approximated by a first-order Taylor series expansion:
437    \begin{equation}
438      \label{eq:jfnkjacvecfd}
439      \mat{J}(\vek{x}^{k-1})\,\vek{w} \approx
440      \frac{\vek{F}(\vek{x}^{k-1}+\epsilon\vek{w}) - \vek{F}(\vek{x}^{k-1})}
441      {\epsilon}
442    \end{equation}
443    or computed exactly with the help of automatic differentiation (AD)
444    tools. \code{SEAICE\_JFNKepsilon} sets the step size
445    $\epsilon$.
446    
447    We use the Flexible Generalized Minimum RESidual method
448    \citep[FGMRES,][]{saad93:_fgmres} with right-hand side preconditioning
449    to solve Eq.\,(\ref{eq:jfnklin}) iteratively starting from a first
450    guess of $\delta\vek{x}^{k}_{0} = 0$. For the preconditioning matrix
451    \mat{P} we choose a simplified form of the system matrix
452    $\mat{A}(\vek{x}^{k-1})$ \citep{lemieux10} where $\vek{x}^{k-1}$ is
453    the estimate of the previous Newton step $k-1$. The transformed
454    equation\,(\ref{eq:jfnklin}) becomes
455    \begin{equation}
456      \label{eq:jfnklinpc}
457      \mat{J}(\vek{x}^{k-1})\,\mat{P}^{-1}\delta\vek{z} =
458      -\vek{F}(\vek{x}^{k-1}),
459      \quad\text{with}\quad \delta\vek{z}=\mat{P}\delta\vek{x}^{k}.
460    \end{equation}
461    The Krylov method iteratively improves the approximate solution
462    to~(\ref{eq:jfnklinpc}) in subspace ($\vek{r}_0$,
463    $\mat{J}\mat{P}^{-1}\vek{r}_0$, $(\mat{J}\mat{P}^{-1})^2\vek{r}_0$,
464    \ldots, $(\mat{J}\mat{P}^{-1})^m\vek{r}_0$) with increasing $m$;
465    $\vek{r}_0 = -\vek{F}(\vek{x}^{k-1})
466    -\mat{J}(\vek{x}^{k-1})\,\delta\vek{x}^{k}_{0}$
467    %-\vek{F}(\vek{x}^{k-1})
468    %-\mat{J}(\vek{x}^{k-1})\,\mat{P}^{-1}\delta\vek{z}$
469    is the initial residual of
470    (\ref{eq:jfnklin}); $\vek{r}_0=-\vek{F}(\vek{x}^{k-1})$ with the first
471    guess $\delta\vek{x}^{k}_{0}=0$. We allow a Krylov-subspace of
472    dimension~$m=50$ and we do not use restarts. The preconditioning operation
473    involves applying $\mat{P}^{-1}$ to the basis vectors $\vek{v}_0,
474    \vek{v}_1, \vek{v}_2, \ldots, \vek{v}_m$ of the Krylov subspace. This
475    operation is approximated by solving the linear system
476    $\mat{P}\,\vek{w}=\vek{v}_i$. Because $\mat{P} \approx
477    \mat{A}(\vek{x}^{k-1})$, we can use the LSR-algorithm \citep{zhang97}
478    already implemented in the Picard solver. Each preconditioning
479    operation uses a fixed number of 10~LSR-iterations avoiding any
480    termination criterion. More details and results can be found in
481    \citet{lemieux10, losch14:_jfnk}.
482    
483    To use the JFNK-solver set \code{SEAICEuseJFNK = .TRUE.} in the
484    namelist file \code{data.seaice}; \code{SEAICE\_ALLOW\_JFNK} needs to
485    be defined in \code{SEAICE\_OPTIONS.h} and we recommend using a smooth
486    regularization of $\zeta$ by defining \code{SEAICE\_ZETA\_SMOOTHREG}
487    (see above) for better convergence.  The non-linear Newton iteration
488    is terminated when the $L_2$-norm of the residual is reduced by
489    $\gamma_{\mathrm{nl}}$ (runtime parameter \code{JFNKgamma\_nonlin =
490      1.e-4} will already lead to expensive simulations) with respect to
491    the initial norm: $\|\vek{F}(\vek{x}^k)\| <
492    \gamma_{\mathrm{nl}}\|\vek{F}(\vek{x}^0)\|$.  Within a non-linear
493    iteration, the linear FGMRES solver is terminated when the residual is
494    smaller than $\gamma_k\|\vek{F}(\vek{x}^{k-1})\|$ where $\gamma_k$ is
495    determined by
496    \begin{equation}
497      \label{eq:jfnkgammalin}
498      \gamma_k =
499      \begin{cases}
500        \gamma_0 &\text{for $\|\vek{F}(\vek{x}^{k-1})\| \geq r$},  \\
501        \max\left(\gamma_{\min},
502        \frac{\|\vek{F}(\vek{x}^{k-1})\|}{\|\vek{F}(\vek{x}^{k-2})\|}\right)  
503    %    \phi\left(\frac{\|\vek{F}(\vek{x}^{k-1})\|}{\|\vek{F}(\vek{x}^{k-2})\|}\right)^\alpha\right)  
504        &\text{for $\|\vek{F}(\vek{x}^{k-1})\| < r$,}
505      \end{cases}
506    \end{equation}
507    so that the linear tolerance parameter $\gamma_k$ decreases with the
508    non-linear Newton step as the non-linear solution is approached. This
509    inexact Newton method is generally more robust and computationally
510    more efficient than exact methods \citep[e.g.,][]{knoll04:_jfnk}.
511    % \footnote{The general idea behind
512    %   inexact Newton methods is this: The Krylov solver is ``only''
513    %   used to find an intermediate solution of the linear
514    %   equation\,(\ref{eq:jfnklin}) that is used to improve the approximation of
515    %   the actual equation\,(\ref{eq:matrixmom}). With the choice of a
516    %   relatively weak lower limit for FGMRES convergence
517    %   $\gamma_{\min}$ we make sure that the time spent in the FGMRES
518    %   solver is reduced at the cost of more Newton iterations. Newton
519    %   iterations are cheaper than Krylov iterations so that this choice
520    %   improves the overall efficiency.}
521    Typical parameter choices are
522    $\gamma_0=\code{JFNKgamma\_lin\_max}=0.99$,
523    $\gamma_{\min}=\code{JFNKgamma\_lin\_min}=0.1$, and $r =
524    \code{JFNKres\_tFac}\times\|\vek{F}(\vek{x}^{0})\|$ with
525    $\code{JFNKres\_tFac} = \frac{1}{2}$. We recommend a maximum number of
526    non-linear iterations $\code{SEAICEnewtonIterMax} = 100$ and a maximum
527    number of Krylov iterations $\code{SEAICEkrylovIterMax} = 50$, because
528    the Krylov subspace has a fixed dimension of 50.
529    
530  \paragraph{Elastic-Viscous-Plastic (EVP) Dynamics\label{sec:pkg:seaice:EVPdynamics}}~\\  \paragraph{Elastic-Viscous-Plastic (EVP) Dynamics\label{sec:pkg:seaice:EVPdynamics}}~\\
531  %  %

Legend:
Removed from v.1.19  
changed lines
  Added in v.1.20

  ViewVC Help
Powered by ViewVC 1.1.22