Performance optimization¶

In this notebook, we'll show you how to make icepack run faster if you're finding that simulation speed is becoming a barrier for you. The solvers in icepack have been designed to give an answer as reliably as possible, even with somewhat outlandish input data. For example, we default to solving linear systems using a direct factorization of the matrix because this approach has better guarantees of giving a useful result. If you're willing to spend the time tweaking solver options, you can get a faster simulation using iterative linear solvers. In this demo we'll show you how.

As a test case, we'll use the MISMIP+ experiment as describe in Asay-Davis et al. (2016). MISMIP+ is a standard test problem for evaluating glaciological flow models and we'll revisit it again in the guides that follow.

The MISMIP+ experimental protocol defines a profile for the ice bed that is a 6th-order polynomial in $x$ and an exponential in $y$. The bed shape was designed to create a bench or bedrock high in the $x$ direction that the ice could ground on and, given enough of a meltwater kick, to retreat off of. The shape in the $y$ direction is designed to funnel ice off of the side walls and create geometric constrictions.

The experiment further sets the default values for the ice fluidity and friction, but leaves it to modelers to set the exact friction law. We'll use our form of the Schoof sliding law from the tutorials on ice streams rather than the Weertman sliding law. This sliding law crucially includes a ramping factor that smoothly dials the friction down towards the grounding line. Without this ramping factor, the system would experience a sudden shock across the grounding line, which often manifests in numerical artifacts.

The experiment specifies an accumulation rate of 30 cm / year in the initial spin-up phase. In later guides we'll look at the next phase of the experiment, which uses a pre-defined function of the depth below the water line to set sub-ice shelf melt rate.

For now, we'll focus on the 2D ice stream model, which is more than enough to illustrate the different kinds of solver options you can choose. Comparing the results with those of the hybrid model is a worthwhile exercise though -- it tells us the degree to which resolving vertical shear is important for the problem we're solving.

First stage¶

We don't have a great way to know ahead of time what the steady state of the simulation will be, so instead we'll start with a very blunt initial guess -- a constant ice thickness of 100m everywhere. Then we'll spin this up to steady state, like we've done before for other synthetic problems.

Since we'll be running the same simulation several times, we'll wrap it up in a Python function. This function will take in the solver object, the initial state of the system, the length of the simulation, and the timestep.

The initial state of a constant ice thickness is very weird and unphysical. A consequence of this is that the diagnostic equations become horribly ill-conditioned. So to start things off we'll propagate the model forward by a few hundred years in order to get rid of some of the wilder transients.

For this initial phase of the experiment, we'll use the default nonlinear solver in icepack. This solver is slow but it's more resilient against weird input data than the alternatives. We've shown how to pass otions to the flow solver before in order to specify which parts of the boundary have ice inflow and which parts are the side walls. Here we're passing some extra keyword arguments: 'diagnostic_solver_type', which can be either 'icepack' or 'petsc'; and 'diagnostic_solver_parameters', which specifies a list of options to pass to PETSc describing how to solver linear or nonlinear systems.

Now the interesting part begins.

Experiment¶

Now that we have a more-or-less reasonable state to start from, we can experiment with solver options. Our experiment will consist of running the model forward by 2700 years. This still isn't quite enough to get close to steady state. Near the shear margins, the ice thickness can get very thin -- on the order of 1-10 m -- after a few thousand years. At the low spatial resolution we're using for this run, the solver can even get negative ice thicknesses, although this unphysical feature can be eliminated by clamping the thickness from below at every timestep. We'll show how to increase the resolution in the demos that follow.

Baseline: icepack solver¶

First, let's continue using the default solver from icepack again. We'll then try to beat this time by choosing better solvers from PETSc.

And now let's visualize the results:

PETSc with LU¶

Now let's try using PETSc's Scalable Nonlinear Equation Solver, or SNES. To switch to PETSc's solver instead of icepack's, pass the key-value pair 'diagnostic_solver_type': 'petsc' when initializing the flow solver. You can then pass a number of parameters dictating how to solver linear and nonlinear systems. Here, we've chosen to use a direct factorization method to solver linear systems. This is the exact same method that the icepack solver used. We've also specified a special type of line search -- 'cp' for critical point. This line search procedure assumes that the residual equations are actually the gradient of some nonlinear functional, and in our case they are -- this is exactly the action functional for the ice stream equations. For more information, you can consult the PETSc documentation.

To summarize, we're using the same algorithms with the PETSc solver as we are in the icepack solver. This experiment is more a test of how much we save by doing all of our nonlinear equation solving with PETSc.

Using PETSc's nonlinear solvers shaved off about 1/8 of the total simulation time on the machine I'm using. One likely explanation for some of this performance gain is that round-trips from Python to C and back are relatively expensive. Since the icepack solver is calling out to scipy to do a line search, this round trip happens more frequently than with PETSc, which makes a single call into the C code.

PETSc with multigrid¶

Now let's see if we can do even better. The direct solver that we used before requires the same amount of computational work no matter what. If you already have a good guess for the solution, say by taking the value from the previous timestep, there's no way to take advantage of that information with direct solvers. Iterative solvers, on the other hand, can converge faster with a good initial guess for the solution. This is true of both the diagnostic equation for the ice velocity and the prognostic equation for the ice thickness. The value of the velocity search direction from the previous Newton step or the ice thickness from the previous timestep can be a good starting guess and we should try to use this fact.

The diagnostic equations of ice flow are the derivative of a convex action functional, so the Hessian is a symmetric and positive-definite linear operator. That means that we can use the conjugate gradient algorithm, which has very favorable convergence properties for symmetric and positive-definite systems. Choosing a good preconditioner is just as important as choosing the iterative algorithm. Here we'll be a little blunt and use PETSc's multigrid preconditioner. The multigrid method can be a very efficient option for symmetric, positive elliptic systems like the linearization of the diagnostic equations. But there may be other choices that work as well or better and, if you are concerned about performance, you should try others. PETSc offers loads of different preconditioners; to see an exhaustive list, you can consult their documentation.

The prognostic equations don't have the nice symmetry properties of the diagnostic equations. We'll instead use the generalized minimum residual (GMRES) algorithm together with an incomplete LU preconditioner.

This shaves off a further 1/3 of the total simulation time. There's a lot of room for dialing in the multigrid solver as well; for example, you can decide to use V- or W-cycles, and you can set the number of levels. You might also be able to speed up the prognostic solver using a multigrid method, but this is likely to require more hand-tuning than for the diagnostic problem.

Conclusion¶

There are a huge number of ways to accelerate nonlinear solver performance beyond what we've shown here. PETSc implements just about anything you could dream of and then some. To see what else is out there, you can consult the full list of linear solvers and preconditioners. You can access these different solvers and preconditioners by passing the right solver parameters to the diagnostic or prognostic solver.

Making choices about solution methods is often a tradeoff between robustness and speed. For this simulation, we had to use the more robust but slower solver in icepack for the first 900 years of the spin-up because the equivalent line search solvers in PETSc would fail. (They would actually work fine on some machines but fail on others, because computing is hard.) The trust region method from PETSc does work for the first 900 simulation years as well, but can be slower on more reasonable inputs. If you are interested in tuning performance, you should expect to use multiple solution methods in the same simulation, especially if the system can meander through multiple regimes.

We've focused here on how to get the most performance out of a single simulation at a fixed resolution. Things get much more interesting when we need to get an accurate result and are free to increase the resolution as need be. For these types of problems, a lot of the work can be done at coarser resolution and knowing exactly how much is both an art and a science.