Priorities and Plans (pre July 2005)
Summary of priorities and plans for software development associated with the Short-Term Crustal Dynamics working group. This summary was developed before the July 11-14, 2005, workshop and is superceded by the workshop report.
Working Group Objectives
- Develop software for simulation of crustal dynamics on spatial scales ranging from meters to hundreds of kilometers and temporal scales ranging from milliseconds to thousands of years. The software will be open-source, well-documented, benchmarked, and maintained by members of the working group and the CIG community.
Brad Aagaard and Charles Williams have been working with Matt Knepley and Michael Aivazis to integrate their codes (EqSim and Lithomop) into the Pyre framework. Charles is maintaining functionality in Lithomop while making changes and is releasing updated versions. Brad decided a more coherent organization of EqSim was in order and chose to write the Pyre version from scratch while using the previous version, which is fully functional, for current projects. Brad and Charles identified a significant amount of commonality between the two codes and decided to coordinate their development efforts with a plan to merge their codes into a single suite of modules, which will be called PyLith. Both codes will gain additional functionality as a result of the merge.
- Matt Knepley added Python bindings to PETSc. This significantly reduces the amount of code required to write solvers in the Pyre framework using PETSc because nearly all of the code can be written in Python. Without the bindings, one would have to implement the solver in C, C++, or Fortran and write Python bindings.
- Lithomop development highlights
- Added ability to automatically compute gravitational prestresses. This avoids the common problem in modeling tectonic problems where unrealistically large deformations occur when gravity is “turned on”. Material models have been derived using this new procedure (linear elastic, linear Maxwell viscoelastic, and power-law Maxwell viscoelastic). The first two of these have been included in the code, and the third will be added when the code structure for material models has been firmly established.
- Top-level code has been reorganized in preparation of a parallel version (anticipated release date is early Fall 2005). Restructuring of code has led to much cleaner implementation of boundary conditions and initial conditions as well as material models. Numerous small bugs have been fixed.
- Use of PETSc for matrix/vector storage and operations. Using the PETSc linear solver reduced computation time by a factor of 2 over the previous implementation that used BLAS alone! PETSc options can be set from the command line; this provides access to PETSc debugging tools which reduces code development time.
- Updated code to use latest version of the Pyre framework (pythia-0.8). Pythia-0.8 includes many enhancements to the user interface for configuring simulations. Fixed compilation problems so code now runs under Darwin and Linux (and presumably other Pyre build procedure supported systems).
Although EqSim is not fully functional, the core finite-element routines are finished. Faults and boundary conditions still need to be implemented.
- EqSim development highlights
- Added ability to use different types of finite-elements. Currently 8-node hexehedral and 4-node tetrahedral elements are implemented. Higher order hexahedral and tetrahedral elements will be added soon.
- User can specify finite-element basis functions and/or integration points at runtime. The default values are the traditional basis functions and Gauss quadrature points.
- PETSc is used for matrix and vector storage and operations. This should reduce the memory usage as a result of the efficient block sparse matrix format in PETSc. We anticipate the computation time will also be significantly reduced by using the matrix and vector operations provided by PETSc compared with a “user” implementation.
- Code includes unit testing to verify modules function properly during development and continue to function whenever code is modified.
Short-Term Software Development Priorities (2005-2006)
- A portable build procedure that makes downloading and building CIG software easy. The build procedure must be able to support software composed of multiple packages and a wide variety of system architectures.
- Incorporate mesh generation capabilities into the framework that are suitable for geophysical problems. Such software must be able to handle complex geometry such as topography and fault surfaces and produce high quality meshes. Discretization size should be able to be set according to user-selectable parameters, such as shear-wave speed or distance from a fault.
- Complete development of a software package for managing parameters specified over a spatial domain (e.g., traction and displacement boundary conditions, material properties, fault slip and tractions). The package should also include the ability to transform between local cartesian coordinate systems and geographic coordinate systems to permit specifications of boundary conditions in either local or geographic coordinate systems. The package should also allow association of units with parameters. A preliminary package called spatialdata-0.1 written by Brad Aagaard, which is in the CIG cvs repository, could be used as a starting point.
- Creation of guidelines and style guides for CIG software development to facilitate (1) compatibility among software packages, (2) promote standardization of interfaces, and (3) ease code understanding by users. The guidelines should also outline expectations for code documentation.
SHORT-TERM CRUSTAL DYNAMICS
- Complete integration of Lithomop and EqSim into the Pyre framework. Begin merging Lithomop and EqSim into a cohesive suite of modules (PyLith) for modeling crustal dynamics.
Brad and Charles can write most of the EqSim/!Lithomop/PyLith code with guidance from CIG software developers.
However, there is an immediate need for a software package for preprocessing finite-element meshes. The package would include (1) importing of meshes from widely used mesh generation tools into common data structures, (2) partitioning a mesh and associated information across processors, (3) global refinement of elements (i.e., uniformly refine a mesh to reduce the node spacing by a factor of 2, 4, etc), and (4) creation of cohesive elements (for modeling dislocations across surfaces such as faults). Most of these capabilities exist, but they are neither in one coherent package nor in the Pyre framework (with the exception of Michael Aivazis’s out-of-date tetra package).
- With guidance from the CIG software development staff, the working group should define standard formats, data structures, and/or interfaces to allow exchange of information between modeling codes and modules. For example, it might be possible to define a simple data structure for exchange of finite-element mesh information between meshers and codes, and codes and visualization packages.
- Implement an Okada dislocation code in the framework. This widely used solution would benefit from a portable, efficient implementation. It could also serve as a hands-on example for people learning the Pyre framework. This could be done by someone in the working group with guidance from CIG software developers or by one of the CIG software developers with some guidance from the working group leaders.
- Create an archive of benchmark results for the various modeling codes. Features would include (1) complete specification of the benchmark problems, (2) archive of benchmark results from various codes, and (3) definition of metrics for determining “goodness of fit”.
- Extend !Lithomop/PyLith and/or other codes with similar capabilities to use spherical coordinates. Input and output would use spherical coordinates. This would facilitate user specification of boundary conditions and incorporation of gravity for domains where the curvature of the Earth must be included.
- Make numerical codes extensible/modular so that it is easy to replace constitutive equations with user defined ones. Implement the most common elastic, viscoelastic, viscoplastic, and porelastic constitutive equations.
Long-Term Software Development Priorities (2-3 years)
- Enhance software donated to CIG to permit simulation of multiple earthquake cycles with sufficient resolution to capture the buildup of strain in the crust, strain release in propagating ruptures that radiate seismic waves, and postseismic relaxation of the crust.
- Develop infrastructure necessary to couple crustal dynamics software to other CIG software for comprehensive simulation of Earth processes. Current candidates for coupling include CitCom, Snac, Snarc, and SpecFEM3D.
- Incorporate formalized data assimilation techniques into crustal dynamics software. This would promote integration of Earthscope data into crustal dynamics modeling.
- Develop or enhance codes to accommodate large deformation problems. Such codes would help bridge the gap between short-term and long-term crustal dynamics modeling.