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 |
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\\ |
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 |
%---------------------------------------------------------------------- |
%---------------------------------------------------------------------- |
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 |
%---------------------------------------------------------------------- |
%---------------------------------------------------------------------- |
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 |
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 |
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 |
% |
% |