Background and rationale

Adaptive mesh refinement (AMR) has long been identified as a key technology that would aid in the accurate and efficient numerical solution of a number of geodynamics applications. It has been discussed in the geodynamics community for several years and has been a continuous topic on the task list of CIG since its inception. Yet relatively little has happened in this direction so far. Only recently have there been attempts to use AMR in geodynamics: CIG sponsored a workshop on AMR technique in Boulder in October 2007, and a collaboration between George Biros, Omar Ghattas, Mike Gurnis, and Shijie Zhong's groups is currently developing a parallel adaptive mantle convection solver.

One of the reasons for the slow adoption of AMR techniques in geodynamics is the relatively steep initial hurdle: codes have to provide the data structures and algorithms to deal with adaptive meshes, finite elements have to be able to deal with hanging nodes, etc. To do so efficiently and in sufficient generality adds several 10,000 lines of code to finite element programs, too much for the average student to do within the time frame of a dissertation. On the other hand, there are libraries that provide the infrastructure code on which applications supporting AMR can rapidly be built; most of the Boulder workshop was in fact devoted to discussions of one of these libraries, deal.II (see www.dealii.org), a very large open source finite element library supporting AMR, a wide variety of finite elements, and tools to quickly build linear solvers and the linear algebra that underlies it. The author of this proposal is also a maintainer and the principal author of deal.II.

Herein, we propose to continue developing deal.II into a library with wide applicability to geodynamics. In particular, we propose to extend the already existing tutorial of deal.II — currently containing 27 extensively documented programs that show how deal.II can be used for a wide variety of problems — by a suite of programs that demonstrate using this technology for geodynamics applications. In consultation with the geodynamics community, the following areas have been identified where AMR would be of interest:

• Mantle convection: A simple mantle convection application was already demonstrated at the Boulder workshop. Since thermally driven convection shows significant localization, AMR is destined to improve the available accuracy significantly.
• Magma dynamics: For the Magma Dynamics Demonstration Suite currently being built, Spiegelman and co-workers have already identified a number of prototypical equations that are known to benefit from AMR. This includes most importantly the McKenzie equations in one of its forms to simulate the migration of a mixture of solid and molten rock.
• Incompressible Stokes flow: Stokes flow is used to model a number of applications in geodynamics, from being a part of the Boussinesq equations used in mantle convection and thermo-chemical convection to a formulation of incompressible elasticity used in describing antishear problems. A flexible and extendible Stokes solver may therefore be the basis for applications in a number of fields.

deal.II already contains several tutorial programs that are also related to geodynamics applications, for example a quasistatic elasticity solver, and that may also be extended.

The goal of the work proposed here is to write programs for each of the topics above. Continuing in the style of the existing tutorial programs — an extensive introduction explaining the background and formulation of an application as well as the concepts of the numerical scheme used in its solution; detailed comments throughout the code explaining implementation details; and a section showing numerical results — we intend to provide the resulting programs as well-documented applications solving model problems. In particular, they are aimed at the following goals:

• Starting points: The existing tutorial of deal.II has proven to be an excellent starting point for graduate students and researchers to jump-start developing their own applications. By providing programs that are already close to the targeted application, first results can often be obtained very quickly, both maintaining the initial enthusiasm during development as well as allowing to spend research time on implementing application specific behavior rather than using months of work on basic infrastructure code supporting AMR.

  Supporting this point is the fact that although there are currently at least 135 publications presenting results obtained with deal.II, we are aware of only a handful of applications that have been built with deal.II from scratch; all others have started as modifications of one of the tutorial programs.

• Training: The tutorial programs we propose to write will provide students and researchers with a reference implementation of current numerical technology such as AMR, higher order elements, sophisticated linear and nonlinear solvers, stabilization techniques, etc. Providing these as starting points for further development by others will also serve the goal of training a new generation of geodynamicists in modern numerical algorithms.

• Extending equations and formulations: In deal.II, it is fairly simple to extend a set of equations by another equation, for example an additional advected quantity that enters the existing equations as a right hand side or in one of the coefficients. Since applications typically use blocked matrices rather than the one-big-matrix-for-everything approach, it is also not complicated to find suitable linear solvers for augmented equations. Consequently, deal.II is a good tool for trying out more complex formulations of problems, or more complete models and their effects on the accuracy of solutions.

• Rapid prototyping and benchmarking: deal.II provides many interchangeable components that allow rapid prototyping of finite element kinds and orders, stabilization techniques, or linear solvers. For example, typically only a few lines of code have to be changed to replace low-order by high-order elements. Through this, it becomes relatively simple to try out higher order elements, a different block elimination solver, or a different stabilization technique. In turn, this may help in benchmarking applications both regarding computing times to solve as well as concerning the accuracy of numerical solutions.

  The applications proposed herein will already have been benchmarked for correctness. Existing tutorial programs typically employ simpler rather than more complicated solver schemes for exposition but frequently suggest more complicated schemes including hints on how they might be implemented in an appendix.

• Try algorithms: The rapid prototyping abilities of deal.II may also help in determining best algorithms on the scale of programs to which deal.II is applicable, and then to implement this particular algorithm (without the ability to change it easily) in a dedicated program that can run on larger scale machines. For example, a small mantle convection code built on deal.II may be used to determine whether second order elements are useful for this purpose. If so, then one may use this result to implement second, rather than first, order elements in the mantle convection code that Ghattas and Zhong are building and that may run on 1000s of processors, a range currently unattainable by deal.II.

Proposal

We propose to create tutorial programs in the following areas (tentative timeline in parentheses, with a justification given below):

• An incompressible Stokes flow solver, with test cases relevant to applications in geodynamics such as a corner flow to simulate mid-ocean ridges (4/1/2008 - 5/15/2008);
• A Boussinesq solver for thermally driven convection, with application to a simple mantle convection model, based on the code already presented at the Boulder workshop (5/16/2008 - 8/31/2008);
• A program for a mantle convection problem that incorporates more realistic nonlinear material behavior and suitable nonlinear solvers (9/1/2008 - 11/30/2008);
• A solver for thermochemically driven convection, with a setup that mimics the benchmark used in the comparison paper by van Keken (12/1/2008 - 2/28/2009);
• A solver for magma dynamics based on a suitable formulation of the McKenzie equations (3/1/2009 - 6/30/2009);
• A program for quasistatic elastic or viscoelastic behavior in domains that induce singularities and would therefore benefit from adaptivity (7/1/2009 - 10/31/2009).

All programs will be developed in consultation with scientists working in the relevant fields, in particular with Magali Billen, Peter van Keken, Scott King, Laurent Montesi, and Marc Spiegelman, all of whom have expressed interest in this. Development of each application will include a 3-5 day visit at the institution of one of the mentioned colleagues to clarify equations, boundary conditions, benchmark descriptions, and physically reasonable setups. Other researchers may be consulted as necessary. Proponent is also willing to host graduate students for a certain period of time to train them in the use of deal.II and the newly developed programs.

Each of the programs will include state-of-the-art adaptive finite element schemes, modern solver techniques, and extensive documentation that includes an introduction explaining the background and formulation of an application as well as the concepts of the numerical scheme used in its solution; detailed comments throughout the code explaining implementation details; and a section showing numerical results. The resulting programs and accompanying documentation will be released as part of the open source software deal.II, both as part of the publicly accessible Subversion repository as well as future releases. CIG may, at its discretion, choose to also host these programs on its website. Experience has shown that developing an advanced tutorial application including documentation, verification of correctness, benchmarking, and implementation of realistic test cases takes on average on the order of four weeks of work. Taking into account that not all of the programs are entirely independent of each other, this proposal asks for 5.5 months of salary over the course of 1 1/2 years, from 4/1/2008 to 10/31/2009, as well as support for travel to the home institutions of abovementioned colleagues. The tentative timeline given above reflects a roughly 1/3 effort.