FAULTS a thin-plate, flat-Earth finite element for anelastic crustal deformation (and PLOTFAULTS) (a graphics program to plot both input and output datasets on Versatec) (and Faults2AI) (a graphics program to plot input and output in Adobe Illustrator PostScript) released August 1993 updated 5 March 1999 by Peter Bird Department of Earth and Space Sciences University of California Los Angeles, CA 90095-1567 (310) 825-1126 pbird@ess.ucla.edu Boilerplate ----------- (C) Copyright 1991, 1992, 1993, 1997, 1998, 1999 by Peter Bird and the Regents of the University of California. Colleagues are free to use, modify, and distribute this software at no charge, but subject to the following conditions: 1. If by using these programs you can obtain scientific results worthy of publication, you must acknowledge the source of the program in any seminars, proposals, and papers you prepare. 2. The programs and/or source code shall not be modified in a way which hides their author or origin. 3. These notices shall be distributed together with any copies of the source code or executable programs. While you may recover your actual costs of blank media and handling, no other fee or price shall be charged for copies you distribute. 4. By using this package, you assume an ethical obligation to help the author improve it, by reporting any bugs or defects together with a complete accounting of relevant circumstances. Purpose ------- To solve the momentum equation (stress-equillibrium equation) under the quasistatic (creeping-flow) approximation, for thin faulted plates of crust deforming anelastically on a flat Earth. To calculate long-term average velocities, long-term average fault slip rates, long-term average anelastic strain rates, and long-term average deviatoric stresses. (Here, "long-term" means averaged over several seismic cycles.) To provide either colored or black-and-white maps of the input data (finite element grid, elevation, heat flow, crustal thickness) and the computed results (velocity, slip rate, strain rate, stress). Limitations: ------------ In program FAULTS, the bottom of the model domain is the base of the crust (the Moho). If you wish to include mantle lithosphere in the domain, then you need program PLATES, available in directory .../neotec/plates. (It is especially important to include the strength of mantle lithosphere when modeling oceanic plates.) If the region you are modeling is so wide that a flat-Earth approximation is unacceptable, then you need programs ORBWEAVE and SHELLS, available in directory .../neotec/shells. (However, these programs use less sophisticated elements, and are therefore less accurate when modeling a local region like one state.) Programs FAULTS, PLATES, and SHELLS are all for "snapshot" models of instantaneous velocity at the present. They do not time-step or compute finite strains, so they cannot be used for historical simulations. To do this, obtain the finite element programs of the LARAMY group from directory .../paleotec. (Unfortunately, these programs do not permit fault elements within the plate- only a single subduction zone at the edge of the grid is allowed.) All of my finite-element programs are "thin-plate" programs which use the isostatic approximation in place of the vertical component of the momentum equation. Therefore, they do not model anisostatic effects like plate flexure in a trench, or around a seamount. All of my finite-element programs are anelastic. Therefore, they cannot be used to model earthquakes or post-earthquake relaxation, for which you will need a (nonlinear) viscoelastic program like Melosh's TECTON. References ---------- Bird, P. (1988) Formation of the Rocky Mountains, western United States: A continuum computer model, Science, v. 239, p. 1501-1507. Bird, P. (1989) New finite element techniques for modeling deformation histories of continents with stratified temperature-dependent rheologies, J. Geophys. Res., v. 94, p. 3967-3990. Bird, P. (1992) Deformation and uplift of North America in the Cenozoic era, in: K. R. Billingsley, H. U. Brown, III, and E. Derohanes (Ed.), Scientific Excellence in Supercomputing: the IBM 1990 Contest Prize Papers, v. 1, Baldwin Press, Athens, Georgia, p. 67-105. Bird, P., and X. Kong (1994) Computer simulations of California tectonics confirm very low strength of major faults, Geol. Soc. Am. Bull., v. 106, p. 159-174. Bird, P. (1996) Computer simulations of Alaskan neotectonics, Tectonics, v. 15, p. 225-236. Directory PATH listing for Volume FAULTS ---------------------------------------- \: (the "root" directory) READ_ME.8 (this file; different from READ-ME.1, ...) FillGrid.exe (executable of nodal-data utility, for WinNT, Win95) Faults2AI.exe (executable of graphics program for WinNT, Win95) +---SOURCE (a directory of source code in Fortran77) | FillGrid.f (utility program to add nodal data to a FE-grid file) | FAULTS.FOR (the finite element program) | PLOTFAULTS.FOR (the graphics program for Versatec plotter) | Faults2AI.f90 (graphics program with Adobe Illustrator output) | \---WORK (a directory of sample working files) NET1991.FEG (a grid of California, from DRAWGRID) ALIASES (renumbering of nodes needed for some graphics plots) BC1991.BCS (boundary conditions for this grid) LanModel.ai (generic Adobe Illustrator PostScript frame for Faults2AI) P98C813.IN (parameter values and plot-control switches) T98C813.OUT (text output from a run of FAULTS) V98C813.OUT (node-velocity output from a run of FAULTS) CALXY.DIG (a basemap file of California, from DIGITISE) Suggestions: ------------ These programs will not run under DOS because of its 640K memory limit. They can be run under Windows NT or Windows 95 or OS/2 on a fast 80486DX or Pentium PC (using 32-bit Fortran77 or Fortran 90 compilers), The executable of the graphics program (Faults2AI.exe) requires Win NT or Win 95 (or possibly OS/2?). My experience is that FAULTS takes about 5 minutes to converge with the sample grid given (1000 nodes) when run on an IBM or Cray vector computer. The graphics program needs about 20 seconds to warm up and then 5-15 seconds per plot after that. The times for FAULTS will rise as the square of the number of nodes, approximately. Therefore, you might prefer a mainframe or a fast workstation. (Memory requirements are not as serious; I think you can work with a grid of 1000 nodes in 5 Megabytes.) Because of this, both FAULTS and PLOTFAULTS are written as batch- job programs which take all input from files (actually, from Fortran I/O device numbers, which you assign to files) and leave all output in files. There is no interaction. Faults2AI can be run as a batch job (if you preconnect the input and output files to the Fortran device numbers), or it can be run as an interactive program which will prompt you for the names of needed files. The interaction is text-mode only for maximum portability. To ease the job of porting this code to different machines, I have tried to use standard Fortran77 throughout. One place where I have cheated is in PLOTFAULTS, where I used "internal writes" to convert floating-point numbers to text strings without the bother of writing them to an external file. If your system doesn't allow this, just write the numbers to a disk file using a FORMAT like (F8.3), and read them back with a format like (A8). Faults2AI.f90 is in Fortran 90, if you need to compile it for a system other than a PC. Another place is in the nodal-data utility program FillGrid, which uses some Fortran 90 code mixed in with inherited Fortran 77 code. If you read my code, you will quickly discover that I use up to four dummy parameters in subroutine calls: INPUT, OUTPUT, MODIFY, and WORK. These are never assigned any values or used; they are just place-holders to identify the groups of parameters which follow. INPUT-class parameters must have a value before the call, and are not changed. OUTPUT-class parameters are computed by the subprogram, overwriting any initial value. MODIFY-class parameters must have a tentative value before the call, but the subprogram may change this value. WORK-class parameters do not have useful values either before or after the call, but merely provide necessary workspace for the subprogram. This workspace is dimensioned in the main program, so that PARAMETER statements may be used the set all the dimensions consistently. If you are unsure what values to assign to the integer array dimensions in the PARAMETER statments, then just set them to 1, compile FAULTS (or PLOTFAULTS), and try to run it. The program will detect a problem, and notify you how large the parameter needs to be. Then you can change it, recompile, and try again. I wish there was a better way to do this! Since I have been working on an IBM mainframe, normal floating-point numbers have 4-byte (32-bit) precision, and DOUBLE PRECISION numbers have 8-byte (64-bit) precision. If you move to a system where the default is 64-bit precision, then by all means you should remove the DOUBLE PRECISION statements. At the heart of both programs is a solver of linear systems. Since this needs to be both accurate and fast, I didn't write my own. In FAULTS, I used routines from IBM's Engineering Sciences Subroutine Library (ESSL) which are hand-optimized for the IBM vector processors. In Faults2AI, I used an IMSL routine, LSLPB, from the MicroSoft version of the IMSL library that is supplied with DIGITAL Visual Fortran (Pro). You should search out the best solver for banded real linear systems on your particular platform, and substitute it. The routines to be supplied are replacements for DGBF and DGBS called from FAULTS (8-byte factor and solve a banded real linear system, respectively), and SPBF and SPBS called from PLOTFAULTS (4-byte factor and solve a banded real linear system, respectively). In searching out replacements, bear in mind that the linear systems in FAULTS do not have symmetric coefficient matrices. Now, you will quickly discover that each solver for banded linear systems assumes a different ideosyncratic storage system for the coefficient matrix! Therefore, I have made this easy to change. Instead of letting Fortran compute the storage locations of a particular element (particular row, particular column) of the coefficient matrix, I have done it myself with a statement function INDEXK. You only need to change this statement function when you change solvers. (Yes, it would be more modular to make INDEXK a subprogram function instead of a statment function, because then it would only have to be changed in one place. Unfortunately, that is not only slower because of the many calls performed, but doubly slower because it blocks vectorization of any loops from which it is called.) In PLOTFAULTS, the actual graphics are produced by calls to a library of Fortran77-compatible subprograms for the Versatec Spectrum 11" electrostatic plotter. This does reasonably good work in either grey- and-white or color (the former for publication, the latter for slides), and does it quickly and cheaply. If you have no access to such a device, you might want to supply dummy subroutines to divert these calls and make them do what you want. You might want to put your graphics on a monitor, for instance. Whatever you do, you will have to know what the parameters in my calls to the Versatec routines are supposed to represent. One source of information is the "Versaplot Random Enhanced Programmer's Reference Guide", from Verstec, 2710 Walsh Avenue, Santa Clara, CA 95051 (Form 5770-02). Or, you can call me, and I'll look the answers for you in my personal cheat-sheets (I don't actually own this manual, either.) Another (usually, better) option for graphics output is Faults2AI. This is a program that I wrote in Fortran 90 to convert finite-element input and output files to maps expressed as ASCII text files in the Adobe Illustrator dialect of the PostScript graphics language. I supply both source code (Faults2AI.f90) and an executable (Faults2AI.exe) which can be run on Windows NT or Windows 95 (or probably on OS/2, and possibly on Windows 3.1 if you have loaded the Win32 subsystem(?)). The output file(s) from this program are intended to be read by Adobe Illustrator 4 for Windows 3.1, from which they can be edited and printed. The editing option is especially nice: you can hide or emphasize groups of lines, change titles, add annotations and arrows, shrink maps and combine them, or redefine the custom colors to shades that you find more pleasing. In particular, you can redefine the custom color "foreground" from black to white and the custom color "background" from white to black, resulting in very nice slide copy (white letters + colors on black). (While running Plates2AI, you will notice that it must read a file called LanModel.ai. I provide this. It is a generic PostScript file containing Adobe Illustrators macros, custom colors, and custom shading patterns (but no actual plot). You can edit this file with Adobe Illustrator, but I would not advise it unless you really understand AI.) Whichever graphics program you choose, notice that most of the input values controlling the graphics are in the model parameter file (e.g., sample P98C813.IN provided), FOLLOWING the lines read by FAULTS. For each possible plot, there is one line. The leading T/F indicates whether that plot is to be produced (and this will affect which input files are required). The second item is the contour interval, if you want to specify it (leaving this blank means the program will choose). The third item is the variable value that should be associated with the central color (or gray shade) of the color bar; if you do not want to choose this, leave it blank and the program will decide. The 4th item is an integer (+1 or -1) which controls whether dark/blue colors or light/red colors are associated with lower variable values. Following the lines for each plot type are choices about pen weights, the location of the center of the map window, the map scale, and whether the plot is to be black/white/gray (for publication and laser printers) or colored (for slide-makers and/or inkjet printers). Notice that after the graphics block, there are 5 more lines read by Faults2AI which define the conic projection that was used to "flatten" the Earth into the (x,y) model plane. These lines are not read by either FAULTS or PLOTFAULTS, and their ONLY function is to allow Faults2AI to put latitude/longitude fiducial lines on the maps. Therefore, you do not need to be concerned about great accuracy. (If you used some other map projection than conic, you can probably still find an equivalent conic projection that will be close enough for this purpose. If you can't, then just use Adobe Illustrator to select and hide the meridians and parallels on the maps, since they are not accurate.) Before you begin to work with FAULTS, be sure to take advantage of all the documentation. Read Bird (1989) for a general foundation in the finite element techniques. Then, read the Appendix in Bird and Kong (1994) to learn how faults are modeled. You will also find some useful comments there on convergence, and how it is affected by the input parameter OKDELV (in sample file P98C813.IN). Then, read the many comment lines at the beginning of FAULTS. You will find explanations of the PARAMETER statements, the I/O device numbers, and a complete glossary of variables used in the main program. Finally, look at the sample text output file (T98C813.OUT) where comments are added to clarify the meaning of many input variables. The units in all of my files are SI. I suggest that you use the same. My programs are entirely free from units, until the very last steps of PLOTFAULTS and Faults2AI, where m/s velocities are converted to mm/year, and dimensions in m are converted to dimensions in km. If you don't mind recoding the legends in the output plots, you are free to use cgs or English or any other consistent set of units. Once you actually begin modeling a region of your own, you will need to replace subroutine MANTLE of FAULTS, which contains problem- specific estimates of the horizontal components of flow at the top of the mantle. Be sure to also replace the same subprogram in PLOTFAULTS and in Faults2AI (and to recompile the latter), or else your plots of basal drag will be wrong, although everything else in the solutions will be correct! If you are ready to begin preparing a finite-element grid for your own problem, then you should look at the other read_me files in this FTP site. However, here is a short "recipe" for creating a grid, in terms of the programs you should use, in order: Digitize -> DrawGrid -> Number -> FillGrid -> Faults2AI -> Adobe Illustrator. (At this point, you will only be able to plot the grid, elevation, heat-flow, and crustal-thickness nodal data.) When the grid is satisfactory, you can model deformation with FAULTS and then complete the other plot types (velocity, fault slip rate, strain rate, stress, etc.) with Faults2AI and Adobe Illustrator. GOOD LUCK! Peter Bird UCLA 29 April 1998