GeoFEST Documentation

• New!! Documentation for special release 4.5(G). This version of the program was developed in order to demonstrate and study the use of solution-driven adaptive mesh refinement (using Pyramid) in large viscoelastic finite element problems. 4.5(G) was used in March 2005 to study performance of a 16 million element mesh distributed over a cluster of 490 Intel Pentium processors.
  
  
• The GeoFEST Design Document describes the physical model underlying GeoFEST, including the quasi-static mathematical equations for visco-elastic materials, and a finite element approximate solution.
  
  
• The GeoFEST Source Structure and Subroutines describes the code structure and the functional routines used by GeoFEST.
  
  
• The User's Guide appears below.

There is also a detailed example explaining the format of an input file for GeoFEST.

User's Guide

Introduction
In order to simulate viscoelastic stress and flow in a realistic model of the earth's crust and upper mantle, the modeling technique must be able to accommodate a vertically layered structure with imbedded faults which may cut ar arbitrary angles. Stress and displacement features will vary most rapidly near the faults and particularly near fault-terminations. These features argue for fully three-dimensional finite element modeling in the time domain. Two-dimensional modeling, semianalytical techniques, finite difference and semispectral methods either cannot model significant features or geometry of interest, or require gross oversampling in regions of little interest, leading to impossible computational requirements.

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.

Features
The primary quantity computed by GeoFEST is the displacement at each point in a domain. The stress tensor is also computed as a necessary byproduct. The computational domain represents a region of the earth's crust and possibly underlying mantle. It is typically a square or rectangular domain in map view, with a flat upper free surface and constant depth, but the domain may deviate from this. The only requirement is that it be a bounded 3D domain with appropriate surface boundary conditions to render the problem well defined. These boundary conditions may be specified as surface tractions and/or displacements, which are usually specified on all surfaces and possibly on interior surfaces such as faults. Free surfaces have zero surface traction by definition. Faults are interior surfaces, and may have associated dislocation increments at set times, when a stress threshhold is exceeded, or according to some other fault friction/triggering model. The solid domain may contain layers or other distributions of material with associated rhelogical properites. Currently supported materials are isotropic, newtonian elastic, newtonian viscoelastic, and nonnewtonian power-law viscosity.

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.

Theory of operation
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.

Input/Output
The user (perhaps via higher-level tools) is required to specify the model domain, constraints, temporal effects, and other solution controls in a single input file. Note that the model domain is a discretized version of the continuous domain described above. Therefore the volume is defined by its tetrahedral/hexahedral elements (or quad/triangular in 2-D), faults and other surfaces are sets of tetrahedral facets (triangles), and many constraints are specified on nodes. Higher-level constructs such as lines or surfaces may not have explicit descriptive records (although there is a construct for element groups sharing material properties; read on). This input file may be roughly characterized by the following brief list of items, which is in file order:

• An output filename
  
• Arbitrary comments
  
• Number of processors and node decomposition (unsupported)
  
• Number of free displacement dimensions
  
• Number of time periods (0 for pure elastic; else number of steady epochs)
  
• Flags: "Save shape functions", "Solver option"
  
• Node free/fixed axis flags
  
• Node coordinates
  
• Node imposed conditions: displacement, velocity
  
• Properties for clusters of elements: type, integration rule, material
  
• Element records (nodes, property ID) (may be a fault element)
  
• Split nodes
  
• Surface tractions
  
• List of nodes and elements for results printing
  
• Number of time groups, plus reform interval, save time, quake time
  
• Time group start times, implicitness parameter, time step
  
• Output times (list, or flag for every earthquake)
  
• Restart file name (optional) (i.e., start where another run left off, using this file)
  
• Save-state file name (optional) (i.e., save this run so another run can start where we left off).
  

Outputs come in three flavors:

1. Standard Out: comments, summary info, progress messages.
   
2. Named output file (the name is in the first line of the input): the requested displacement and stress information at the requested times.
   
3. 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.