Index: /issm/trunk/bugs/LarourJGR11b.tex =================================================================== --- /issm/trunk/bugs/LarourJGR11b.tex (revision 8037) +++ /issm/trunk/bugs/LarourJGR11b.tex (revision 8037) @@ -0,0 +1,325 @@ +%Page 3 is blank + +%Preamble{{{1 +\documentclass[jgr]{agutex} + +%pictures +\usepackage[dvips]{graphicx} +\setkeys{Gin}{draft=false} + +%Additional packages +\usepackage{amsmath} %for the multline environment +%}}} +%Authors and affiliations {{{1 +\authorrunninghead{Larour et al.} \titlerunninghead{Pine Island Glacier Sensitivity Analysis\copyright2011. All rights reserved} + +\begin{document} + +\title{Sensitivity Analysis of Pine Island Glacier ice flow using ISSM and DAKOTA.} + +\authors{E. Larour\altaffilmark{1}, J. Schiermeier \altaffilmark{1}, E. Rignot \altaffilmark{2,1}, H. Seroussi\altaffilmark{1,3}, M. Morlighem\altaffilmark{1,3}} + +\altaffiltext{1}{Jet Propulsion Laboratory - California Institute of technology, 4800 Oak Grove +Drive, Pasadena, CA 91109-8099, USA.} +\altaffiltext{2}{University of California, Irvine, Department of Earth System Science, Croul Hall, +Irvine, CA 92697-3100, USA.} +\altaffiltext{3}{Laboratoire MSSMat - CNRS U.M.R 8579, \'Ecole Centrale Paris, Grande Voie des +Vignes, 92295 Ch\^atenay-Malabry Cedex, FRANCE.} +%}}} +\begin{abstract}%{{{1 +Assessing output errors of ice flow models is a major challenge that needs to be addressed if we are to increase +our confidence level in projections of mass balance in Antarctica and Greenland. Some of the major constraints +in ice flow models include geometry (thickness and surface elevation of the ice sheet) and boundary conditions (geothermal +flux, basal drag, surface temperature). These constraints are either measured, in which case they carry errors +due to instruments, or inferred using inverse methods (such as for basal drag which is inverted using InSAR +surface velocities) in which case they carry additional errors generated by the inversion process itself. +In both cases, these input errors will result in uncertainties that propagate across a model, and that influence output +results. In order to estimate the resulting error margins on model outputs, we developed a framework based on the +DAKOTA software, which we interfaced to the Ice Sheet System Model. We present results on the Pine Island Glacier, +Antarctica, for which we evaluate error margins of computed mass fluxes for all major tributaries of the ice stream, +given currently known error margins on +thickness, basal drag and rigidity inputs. Our results suggest errors in these inputs transfer linearly into output +errors. This provides an easy way for calibrating measurement requirements for field campaigns collecting data such +as bedrock or surface topography. This type of model also provides a pathway for helping quantify uncertainties +in projections of future sea level rise. +\end{abstract} %}}} + +\begin{article} +\section{Introduction} %{{{1 +%}}} +\section{Model} %{{{1 + +\subsection{Ice flow model} + +The model is based on the two dimensional (2D) MacAyeal implementation \citep{MacAyeal1989,MacAyeal1992b} of the Shallow +Shelf Approximation (SSA) for the ice shelf, and the Shelfy Stream Approximation for the ice sheet, using a Continuous Galerkin +Finite Element Method to solve the stress-equilibrium equations. The Shelfy Stream Approximation is valid for fast flow, where +shear stress are confined mainly near the bed. This model is applicable to the bulk of Pine Island Glacier, where most of +the flow is above 60 m/yr. For the remainder of the basin, we assume the Shelfy Stream Approximation to hold true. + +Inverse method for model setup. \\ +Anisotropic mesh. \\ +Mass flux computation across a flux gate. \\ + +\subsection{Sampling and local reliability model} +The DAKOTA (Design Analysis Kit for Optimization and Terascale Applications) \citep{Eldred2008} toolkit provides +an interface between analysis codes and iterative systems analysis methods. Such iterative methods include +uncertainty quantification, sensitivity analysis, design of experiments, parameter studies, and optimization. For +the current study, the uncertainty quantification using sampling and reliability methods was of most interest. + +In order to perform a sampling analysis for uncertainty quantification, the uncertain input variables are defined +with a statistical distribution. These distributions can be normal, uniform, etc. Repeated analyses are run with +input variables generated from within the distributions, and statistics are calculated on the output responses. Typically +means, standard deviations, and cumulative distribution functions (CDFs) are calculated, even though the output +responses need not be (and in general are not) normal distributions. + +Generation of the input variables can be performed in a number of ways. Probably the most common is traditional +Monte Carlo (MC), +where the variables are generated randomly according to the probably distributions. This has the disadvantage +that values in the tails may be neglected, since those probabilities are very low; however, values in the tails +are often critical when performing uncertainty quantification. A second way of generating the input variables +is Latin Hypercube Sampling (LHS), which is a binned approach \citep{Swiler2004}. The n-dimensional variable +space is divided into bins according +to the given distributions, such that one and only one sample occurs in each bin. This has the advantage of forcing +samples into the tails and is also a more efficient method of sampling for a given level of statistical accuracy. +There is overhead in the computation of the LHS samples, but this is generally negligible compared to the cost of +the analyses. + +For a large number of input variables, the cost of sampling analysis to decrease the confidence intervals of the +output responses to desired levels may be prohibitive, since the number of samples required increases geometrically +with the number of input variables. In this case, a local reliability analysis may be performed +to determine +the input variables that have the most significant effects on the output variables. In essence, this method +calculates a finite difference partial derivative for each output response with respect to each input variable at +their baseline values, which requires only $n+1$ solutions. The +variables that have the largest effects can then be studied further, with sampling or parameter methods, and +those with little or no effect might not be of interest. + +Local reliability methods work best when the problem is linear in the neighborhood of the baseline solution. The +finite-difference step size is typically defined by the user, so if the function is not linear in the area, the +value of the secant will change. If the variable is of primary interest, one-dimensional parameter studies +(or different step sizes) can be used to ascertain its behavior. The local reliability method also considers +only one variable at a time, but n-dimensional parameter studies can be used to ascertain if any input variables +of interest have a coupling effect. + +For the local reliability method, given a response $r$ that is a function of $J$ multiple input variables $X_i$: +\begin{equation} + r = r \left( X_1 , X_2, ... , X_J \right) + \label{eq:responsefunction} +\end{equation} +the sensitivities $\theta_i$ are defined as: +\begin{equation} + \theta_i = {{\partial r} \over {\partial X_i}} + \label{eq:sensitivities} +\end{equation} + +The mean of the output responses is assumed to be the baseline value. If each of the input variables is +independent, the variance $\sigma_r^2$ of the output response can be computed from the well-known error +propagation equation: +\begin{equation} + \sigma_r^2 = \sum_{i=1}^J \theta_i^2 \cdot \sigma_i^2 + \label{eq:variance} +\end{equation} +where the $\sigma_i^2$ are the variances of each input variable. Importance factors $IF_i$ for each input +variable may be calculated by dividing each term of the equation by $\sigma_r^2$: +\begin{equation} + IF_i = {{\theta_i^2 \cdot \sigma_i^2} \over \sigma_r^2} + \label{eq:importancefactor} +\end{equation} +These importance factors provide unitless quantities that add up to unity and therefore can be used to +rank the contribution of the input variables. + +\subsection{Mesh partitioning} + +Sampling, as well as local reliabilty methods are based on updates of input variables, according to their statistical +distribution. In order to use such methods on a two-dimensional variable field, such as thickness or basal drag, which covers +the entire mesh domain, each field must be distributed over a partition of equal area surfaces. Such a partition is shown +in Fig. \ref{fig:Partition}, for the Pine Island Glacier, where partition edges are shown in red, and mesh edges in blue. + +For each partition surface, a specific statistical distribution is specified according to the field being sampled. For example, +if thickness is being considered, error margins on thickness measurements from Ground Penetrating Radar can be used to specify +the 3$\sigma$ standard deviation for that particular partition. Also, the average value of thickness over the partition nodes +can be used to specify the mean value of the input variable. Each node of the mesh that belongs to the partition area will +behave according to this statistical distribution. Sampling will therefore be carried out, not over the entire field, but one +field partition at a time. In effect, this solves the problem of sampling an entire two-dimensional field into by distributing +the field over an entire partition, and considering each partition as a new sampled variable. + +In order to partition a mesh into equal area surfaces, for which at least one node must be present, partitioning algorithms were used. +Three packages were considered in our model: +\begin{itemize} +\item MeTiS: a Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices \citep{Karypis1998}. +\item CHACO: a Software for Partitioning Graphs \citep{Hendrickson1995}. +\item SCOTCH: a Software Package for Static Mapping by Dual Recursive Bipartitioning of Process and Architecture Graphs \citep{Pellegrini1996}. +\end{itemize} + +Mesh partitioning techniques: \\ +- talk about differences between all partitioners \\ +- talk about different types of bisections \\ +- talk about anisotrpic meshing \\ + +For a distributed input variable, such as the thickness mentioned above, the domain must be divided into a +number of discrete regions to vary. The finite element mesh provides a convenient discretization of the +domain; however, varying the input variable over each finite element node or element would be prohibitive for very +large problems. In addition, there is the problem of physical size. For meshes with differently sized +elements, some variables would have an inordinate contribution to the response given the physical +areas over which they extend. + +The partitioners investigated for this project all have the goal of reducing parallel computing time +for matrix solutions in finite-element analysis or other computations areas. They +allow a wide variety of parallel machine architecture, but most perform the partitioning based strictly +on a nodal map of the finite element connectivity. Since the thicknesses and other variables are defined +at the nodes, this worked out well, but they do not consider the physical geography of the +nodes. In order to consider the size of the elements, each element area was divided by its number of nodes, +and that area was added to the weighting of each node for the partitioning. + +The effect of this weighting can be shown with respect to the anisotropic mesh from Fig. \ref{fig:Partition}. +In Fig. \ref{fig:PartHist1}, the mesh has been partitioned without the equal-area weighting. As +a result, the partitions, shown in red, vary widely in area, and the accompanying histograms show +that while the partitions each have nearly the same number of nodes, the areas vary by well over an +order of magnitude. In contrast, in Fig. \ref{fig:PartHist2}, the mesh has been partitioned with +the equal-are weighting. Consequently the partitions, shown in red, look visually similar, and +the accompanying histograms show that the partitions have a different number of nodes by over two +orders of magnitude while the areas are approximately the same. + +For the following partitioners, a mesh of the Pine Island Glacier was created with 955177 nodes +and 1904404 triangular elements in order to test the performance on larger models. The mesh was +isotropic, and with such small elements the area of each element was approximately the same. The +result is that non-weighted and weighted partitions did not differ by much, though they did +affect how the partitioner was run. + +MeTiS is a software package for partitioning unstructured graphs and meshes \citep{Karypis1998}. +It was already available as an external package within ISSM, so it was investigated first. +The METIS\_PartMeshNodal module is used to partition a mesh into equal-size parts, which refers +to the number of nodes. It can produce both an element map and a nodal map for the partitions, +but it does not consider nodal weighting, and some of the partitions were discontinuous. An +example of a discontinous partition is shown in Fig. \ref{fig:Metis}. + +The problem with discontinuous partitions is not from a numerical point of view, but rather a +physical one. Since the various partitioned areas are having sensitivities calculated, it +is not apparent with discontinuous partitions which part or parts is causing the effects. This +becomes more important if an aircraft or spacecraft is being assigned to gather additional data +from the areas with the highest influence. + +Chaco is an older package designed to partition graphs \citep{Hendrickson1995}. It has three +primary global partitioning methods: Kernighan-Lin, inertial, and spectral, all of which can be used +recursively. K-L had difficulty coarsening and tended to produce discontinuous partitions unless no +vertex weights were used; inertial produced fast, good, and straight partition boundaries; and +spectral computed eigenvalues, which led to very long and very discontinuous partitions. + +Chaco was the only partitioner investigated that uses geometrical information, specifically the coordinates +of the nodes. For its inertial method, it computes the principal axis using the nodal locations and weights, +and then divides the structure into equal masses by planes orthogonal to the principal axis. One, three, or +seven planes, denoted by bisection, quadrisection, or octasection, respectively, may be used. Examples of the +three types of division are shown in Fig. \ref{fig:Chaco}. This leads +to a banded appearance for the global partitioning, but when coupled with the Kernighan-Lin local optimization, +can significantly improve. The process is recursive depending on the number of partitions requested. +Fig. \ref{fig:Chaco} shows the recursion that occurs with a larger number of partitions. + +SCOTCH is a more recent package for static mapping, partitioning, and sparse matrix block ordering +\citep{Pellegrini2008}. The gmap module, which maps a source graph onto a target graph, was used for +the partitioning. As shown in Fig. \ref{fig:Scotch}, without vertex weights, SCOTCH produced good partitions +with no discontinuous regions; however, with vertex weights, there were some very small discontinuous areas. + +Both Chaco and SCOTCH allow edge weighting as well as vertex weighting. Edge weighting has the appeal +of being able to force partition boundaries along certain element edges, for example geographical +features. However, this functionality has not yet been utilized. + +Since all three of the partitioners have source code available, tight interfaces were written to pass data +back and forth within memory, rather than reading and writing to disk files. This greatly increased +the performance. Note that for the types of uncertainty quantification analysis described previously, +sampling and local reliability, only one partitioning needs to occur for the entire analysis. The finite +element mesh is not changing, and therefore neither is the nodal weighting. + +\subsection{DAKOTA wrapping} +Wrapping DAKOTA around ice flow model +%}}} +\section{Validation} %{{{1 +Partition independent results. +Hand verification of importance factors. +%}}} +\section{Data} %{{{1 +%}}} +\section{Results} %{{{1 + +\subsection{Low and high band filter} +\subsection{Simple transfer function for error margins} +\subsection{Linear or nonlinear scaling of error margins} +\subsection{Importance factors different according to inputs} +\subsection{Coupling between tributaries} +Entrainment by some tributaries of other tributaries, upstream and downstream influence. +Zero importance factors for basal drag on ice shelves. + + +%}}} +\section{Discussion} %{{{1 +%}}} +\section{Conclusions}%{{{1 +%}}} +\begin{acknowledgments} %{{{1 + This work was performed at the Jet Propulsion Laboratory, California Institute of Technology, + at the Department of Earth System Science, University of California Irvine, + and at Laboratoire MSSMat, \'Ecole Centrale Paris, + under a contract with the National Aeronautics and Space Administration, Cryospheric Sciences + Program. +\end{acknowledgments} +%}}} +%References {{{1 +\begin{thebibliography}{9} +\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi + +\bibitem[{{\it Eldred et~al.\/}(2008)}]{Eldred2008} +Eldred, M.~S., et~al., {DAKOTA}, a {M}ultilevel {P}arallel {O}bject-{O}riented + {F}ramework for {D}esign {O}ptimization, {P}arameter {E}stimation, + {U}ncertainty {Q}uantification, and {S}ensitivity {A}nalysis, {V}ersion 4.2 + {U}ser's {M}anual, {T}echnical {R}eport {SAND}2006-6337, {\it Tech. rep.\/}, + Sandia National Laboratories, PO Box 5800, Albuquerque, NM 87185, 2008. + +\bibitem[{{\it Hendrickson and Leland\/}(1995)}]{Hendrickson1995} +Hendrickson, B., and R.~Leland, The {C}haco user's guide, version 2.0, + {T}echnical {R}eport {SAND}-95-2344, {\it Tech. rep.\/}, Sandia National + Laboratories, Albuquerque, NM 87185-1110, 1995. + +\bibitem[{{\it Karypis and Kumar\/}(1998)}]{Karypis1998} +Karypis, G., and V.~Kumar, {\it A {S}oftware {P}ackage for {P}artitioning + {U}nstructured {G}raphs, {P}artitioning {M}eshes, and {C}omputing + {F}ill-{R}educing {O}rderings of {S}parse {M}atrices\/}, University of + {M}innesota, 1998. + +\bibitem[{{\it MacAyeal\/}(1989)}]{MacAyeal1989} +MacAyeal, D., Large-scale ice flow over a viscous basal sediment - {T}heory and + application to ice stream-{B}, {A}ntarctica, {\it J. Geophys. Res.\/}, {\it + 94\/}, 4071--4087, 1989. + +\bibitem[{{\it MacAyeal\/}(1992)}]{MacAyeal1992b} +MacAyeal, D., The basal stress distribution of {I}ce {S}tream {E}, + {A}ntarctica, {I}nferred by {C}ontrol {M}ethods, {\it J. Geophys. Res.\/}, + {\it 97\/}, 595--603, 1992. + +\bibitem[{{\it Pellegrini\/}(2008)}]{Pellegrini2008} +Pellegrini, F., {SCOTCH} and {LIBSCOTCH} 5.1 {U}ser's {G}uide (version 5.1.1), + {\it Tech. rep.\/}, ScAlApplix project, INRIA Bordeaux Sud-Ouest, ENSEIRB \& + LaBRI, UMR CNRS 5800, Universit\'e Bordeaux I, 351 cours de la Lib\'eration, + 33405 Talence, France, 2008. + +\bibitem[{{\it Pellegrini and Roman\/}(1996)}]{Pellegrini1996} +Pellegrini, F., and J.~Roman, {SCOTCH}: {A} {S}oftware {P}ackage for {S}tatic + {M}apping by {D}ual {R}ecursive {B}ipartitioning of {P}rocess and + {A}rchitecture {G}raphs, in {\it Proceedings of the {I}nternational + {C}onference and {E}xhibition on {H}igh-{P}erformance {C}omputing and + {N}etworking\/}, HPCN Europe 1996, pp. 493--498, Springer-Verlag, London, UK, + 1996. + +\bibitem[{{\it Rignot\/}(2008)}]{Rignot2008} +Rignot, E., {PALSAR} studies of ice sheet motion in {A}ntarctica, in {\it + {ALOS} {PI} Symposium, Nov. 3-7 2008\/}, 2008. + +\bibitem[{{\it Swiler and Wyss\/}(2004)}]{Swiler2004} +Swiler, L.~P., and G.~D. Wyss, A {U}ser's {G}uide to {S}andia's {L}atin + {H}ypercube {S}ampling {S}oftware: {LHS} {UNIX} {L}ibrary/{S}tandalone + {V}ersion, {T}echnical {R}eport {SAND}2004-2439, {\it Tech. rep.\/}, Sandia + National Laboratories, PO Box 5800, Albuquerque, NM 87185-0847, 2004. + +\end{thebibliography} +%}}} +\end{article} +\end{document}