Click here to download the complete GeoFEST installation package, with validation test problems and makefile.
Finite element modeling in three dimensions allows faithful modeling of complex faulting geometry, inhomogeneous materials, realistic viscous flow, and a wide variety of fault slip models and boundary conditions. While there are particular problems which are more efficiently expressed using alternative approaches, finite elements represent the most generally useful method for inhomogeneous elastostatic and viscoelastic problems. Because finite elements conform to (nearly) any surface geometry and support wide variations in mesh density, solutions may be made arbitrarily accurate with high computational efficiency. This flexibility comes with a price tag for the user or the toolbuilder, that of generating and adapting the mesh of elements upon which the solution is computed. When such generation tools are primitive, researchers may spend substantially more time creating a mesh than solving their problem and interpreting the results. Therefore we describe automated tools for creating and adapting the mesh, including using an initial coarse mesh to generate a solution, whose computable error characteristics inform a further mesh generation cycle and produce an efficient and accurate solution.
We describe the supported features of the GeoFEST code as of this writing in the next section. The precise description of the codes' operation and mathematical specification follows in the Theory of Operation section. The entries in input and output files are listed next, enabling a user to set up an input file with a text editor, and run Visco on that basis. Then the existing tools for processing the input and output are described, allowing a more automated generation of an input mesh, and visualization of the output. Finally some examples are given.
Elastostatic solutions are supported, such as computing the displacements and stresses immediately caused by a specified slip distribution on a fault or finding the interior displacement and stress distribution due to a surface traction or displacement. These solutions are not time-dependent.
Viscoelastic solutions are also supported, in which the material flows and relaxes in response to imposed stress, such as an earthquake event. One may compute the viscoelastic reponse to a single event, or to multiple events in a sequence. The sequence may be user- specified, or may entail fault slips that are dynamically determined to occur according to a fault triggering model. Location-specific body forces are supported.
Boundary conditions and solutions apply to a finite-element discretized approximation to this domain. The domain is defined internally as a mesh of space-filling tetrahedral or hexahedral elements, with three components of displacement at each mesh node constituting the solution. Stress is computed for each element, and is element-wise constant for the current linear tetrahedral element type. Surface nodes carry special boundary conditions such as tractions or specified displacements. Nodes on faults are special split-nodes, that define screw or tensile dislocation on the fault without perturbing the mesh geometry. Temporal evolution is by discrete time steps using an implicit solution technique, allowing large time steps without numerical instability.
The code may be used with all nodes, elements, faults, boundary conditions and time history control created or modified by word processor, constituting the only needed input file. For large meshes hand-construction becomes impractical, and we have sought to support tools for automated mesh generation. Initial attempts have focused on having the user specify fault rectangles, materials, and mesh density, and then semiautomated tools produce the ascii input file.
Other supported features include: Specification of temporal epochs, each with differing steady boundary conditions Boundary velocity condition (steady change in the displacement, imposed) Controls for generating output on a subset of nodes/elements Control of implicit integration parameter Ability to shortcut temporal advance by fixing the sparse system for several time steps Control of checkpointing, saving state and allowing restart Checks a control file at each iteration, allowing clean interruption by the user.
It is worth noting that a central theme to finite elements is the set of shape functions for individual elements. For each linear tetrahedron these are quite simple: for each node, define the function that is unity at that node and linearly declines to zero at the entire opposite face. We may also speak of global shape functions: for each element that shares a particular node, define the function that is unity at that node and declines linearly to zero at the opposite faces of each of the sharing tetrahedra. Shape functions so defined are called interpolatory: if we determine a list of displacements at each node, we may interpolate those values to any point interior to any element. Simply multiply each node displacement by that node's shape function, and sum over those product terms that are nonzero (these will be terms from the nodes defining the enclosing tetrahedron of the interpolation point). Such interpolated displacements will be continuous across all element boundaries (tetrahedral facets).
Within the volume of an isolated finite element, appropriate derivatives of the element shape functions weight the spatial derivatives of the stress tensor (defined from the strain and generalized Hooke's law, relying on locally defined element rheology) to enforce local equilibrium. The shape functions constitute a finite degrees-of-freedom approximation to the continuous system, and a volume integral enforces the equilibrium in a weighted average sense. Elements sharing a node contribute such weighted average terms to an equation for a single displacement. Away from boundary conditions and absent body forces, this set of displacement equation terms is set to zero. Body forces and surface tractions add forcing terms to the right-hand side.
The ensemble system of equations so defined is sparse, and if boundary conditions are sufficient the system is closed and solvable by standard sparse matrix techniques (currently by direct Crout factorization, but also possible to do by iterative techniques). Closed boundary conditions are usually easy to obtain; the user does have to ensure that enough of the boundary is constrained that there are no unconstrained body translations or rotations permitted to the domain.
The solution is the elastostatic solution for the posed problem in terms of displacements at each node. Local linear combinations of these displacements with shape function derivatives yield the stress tensor in each element.
The elastostatic solution is required for any viscous relaxation computation. Once the static step is complete, the time evolution of quasi-static viscous relaxation may begin by computing the viscoplastic strain rate, which is directly determined by the stress and the viscosity parameters. Conceptually this rate adds a force term to the right-hand-side of a sparse system similar to the elastostatic equilibrium. In practice, to obtain the advantages of implicit temporal development, terms are rearranged to modify the sparse equation coefficients as well. Each time step involves a solution to a sparse equation system, of similar cost as the elastostatic solution.
Faults are specified as either fault elements or as split nodes. With fault elements the user specifies fault properties such as failure criterion in an input record similar to that for a finite element. With split nodes the user specifies for each node on the fault its direction and amount of slip. The fault will slip at the initial time step, and also add the same increment of slip at each fault failure event, specified by the user. The split node formalism may be considered to represent the screw dislocation at a node as a separate entity from the displacement at that node, but is implemented here as an equivalent increment in the stress affecting the nodes immediately adjacent.
Outputs come in three flavors:
Standard Out: comments, summary info, progress messages.
Named output file (the name is in the first line of the input): the requested displacement and
stress information at the requested times.
Restart file: contains simulation state. This is in ascii, but is intended for machine reading.
The input file affords considerable power to the user, but is very terse and cryptic. Without a full description of its options and features mistakes and misunderstandings inevitably arise. As a guide to the proper construction and formatting of a GeoFEST input file, we provide below a link to a hypertext-annotated template. The example is for a two dimensional problem, in order to facilitate easy visualization of the domain, however the generalization to 3-D is straightforward. Note the convention used in the input file that many vertical lists use the special flag "0 0" to indicate termination.
Example Format of Input File For GeoFEST Finite Element Program
Physics Technical Design Document
GeoFEST Source Structure and Subroutines
Please direct inquiries, requests, suggestions and bug reports to: lyzenga@aig.jpl.nasa.gov.