SISIPHUS: Scalable Ice-Sheet Solvers and Infrastructure for Petascale, High-Resolution, Unstructured Simulations

SISIPHUS is a model and software infrastructure for climate simulations related to the dynamics of ice sheets -- glaciers that straddle
land and sea as ground ice slides into the ocean. Understanding the dynamics of ice sheets is crucial for accurate predictions of long-term
sea levels changes, as the melt of the ice near its grounding line is expected to be the main contributor to the rise of global sea levels.

The dynamics of ice sheets unfolds on the timescale of decades or even centuries, but at the same time, can undergo sudden changes
when the glaciers protruding from land into the ocean disintegrate in a matter of weeks. These sudden changes are conjectured
to be related to instabilities in the ice sheets associated with the geometry of their grounding lines. Thus, accurate predictions of the future
state of ice sheets and their contribution to the rise of global sea levels necessitates long-term simulations
of the ice sheet dynamics, calculations of its (quasi)stationary states and their bifurcations.

Ice sheet models describe a component of the global climate system and are coupled to other models describing the atmosphere and the ocean.
These other models evolve on much smaller timescales making the global system very stiff because of the presence of rapidly-moving gravity waves
in the ocean and other fast modes. Thus, in order to be able to span decade- and century-long scales of ice sheet evolution, we propose an implicit
formulation of ice sheet dynamics, which allows us to side-step the prohibitively short stiff time scales and still to yield accurate long-term predictions.
Rapid changes in the state of ice sheets can be identified via a stability analysis of their quasi-stationary states and the attendant bifurcations.

Presentations

​Mesh-based tools for land ice simulations (presentation to 2010 Land Ice Working Group, Breckenridge, CO, June 2010)

Nonlinear Stokesian model of ice sheets

Ice sheets are modeled as an incompressible Stokesian fluid with a nonlinear power law rheology (rendering it non-Newtonian) under a free-surface
evolution, constrained by the nonlinear slip boundary condition at the bedrock boundary, as well as by the atmosphere and ocean models.
[Details of the model to follow.]

Scalable software infrastructure

The simulation infrastructure is based upon scalable linear and nonlinear solvers encapsulated by the PETSc library
(Portable Extensible Toolkit for Scientific Computation) and the ITAPS collection of services ​https://trac.mcs.anl.gov/projects/ITAPS,
delivering scalable meshing, mesh (re)distribution and mesh data coupling services.
The interrelationship between the different components of the SISIPHUS simulation infrastructure are illustrated in Figure 1 below

Figure 1. Conceptual picture of the SISPHUS simulation infrastructure. PETSc components are in green, ITAPS in yellow, climate models
(including the SISIPHUS ice model) in blue. Existing components are dark, components to be added or extended with new capabilities are in lighter colors.

• Model defines the core nonlinear equations to be used in different simulation scenarios.
• The actual equations to be solved are defined in terms of the model according to the simulation needs:
  • Spin up: equilibration of the core equations with the other models (ocean, atmosphere) fixed in their current states.
  • Transient: time stepping through the ice sheet dynamics above the stiff time scales.
  • Bifurcation: parameter continuation and/or bifurcation identification and tracking
• Inverse: data assimilation and inverse problem solution to identify the model parameters.

• Model and the overlaid simulation scenario define the nonlinear system to be solved (generally, one per time or parameter step).
  • Nonlinear residuals and the (action of) the (approximate) Jacobian of the system are formed by the model.
  • SNES (Newton's method) solves the nonlinear system.
  • KSP carries out the inner linear (iterative) solve using the supplied Jacobian action.
  • PC preconditions the system to ensure rapid KSP convergence.

• [Need to expand on the mesh movement etc]

Physics-based preconditoners

The need for an implicit formulation of the dynamics of the Stokesian ice sheet model places a significant burden on the preconditioning of the linearized system.
Indeed, to be effective, the implicit solve, be it within the (pseudo) time stepping of parameter stepping context, has to converge within a few Newton iterations and each inner linear iteration has to converge quickly. Since the linearization of the nonlinear Stokes model is diagonally-submissive due to the presence of the incompressibility constraint, the system is a saddle-point, which is undesirable from the spectral point of view for the convergence of Krylov subspace methods (e.g., GMRES). Moreover, traditional and readily available preconditioners for elliptic systems, most notably, multigrid, perform very poorly on such systems.

The answer lies in a class of factorization-based preconditioners, that separate the velocity and pressure degrees of freedom.
The diagonally-dominant velocity-velocity coupling block is then effectively preconditioned with multigrid, while the (zero) pressure-pressure
subblock after factorization is replaced by the Schur complement matrix S. However, S is generally both dense and has a dense inverse.

The physics of the problem dictates that, in the infinite-dimensional setting, the inverse of the Schur complement is a singular integral operator
with a known kernel related to the Green's function of the Stokes operator. This circumstance can be very profitably exploited for preconditioning the Schur complement block, if a low-order discretization of the infinite-dimensional inverse is used.
Moreover, only the action of the discretized integral operator is necessary and the matrix is never formed. Even further, applications of singular integral operators, such as the one based on the fast multipole method (FMM), perform very well on modern hardware architectures because of their sharply lower memory bandwidth requirements. In particular, FMM can approach the peak performance of general purpose GPUs (GPGPUs).

Thus a combination of physics-based splitting of degrees of freedom (velocity vs pressure), subsequent block factorization of the Jacobian (formation of the Schur complement) and the fast application of singular integral operators can result in a very effective preconditioner for the Stokesian ice sheet model, especially on the emerging architectures. Furthermore, this motif (splitting-factorization-fast application) is pervasive throughout scientific computing and can be used in many different settings. In particular, the coupling of the ice sheet model to the atmosphere and ocean models results in a similar saddle-point structure that too can be treated via a physics-based splitting. We are developing a reusable software component infrastructure within PETSc that can greatly simplify the construction of such physics-based preconditioners. We summarize this section with the following list.

• Diagonally submissive (saddle-point) systems result from incompressibility constraints (velocity vs. pressure).
• Velocity-velocity block can be effectively treated with multigrid, but not the full coupled system.
• Factorization and separation of degrees of freedom results in a dense Schur complement block for the pressure-pressure coupling.
• Physically, Schur complement S and its inverse are singular integral operators that can be applied without forming the matrix.
• S preconditioner application can be done effectively with the fast multipole method (FMM)
• FMM has excellent (low) memory bandwidth requirements and can approach peak performance on GPU-like architectures.
• The described factorization-based scheme is ubiquitous
• Reusable software components for this class of preconditioners is being developed withing PETSc and can apply to a wide range of problems.

Mesh and Geometry Tools

Under Construction

---

To learn more about Trac, see the WelcomeToTrac page.