24 June 1998 TO: Fellow Earth Scientists and Hackers FROM: Peter Bird Department of Earth and Space Sciences University of California Los Angeles, CA 90045 (310) 825-1126 pbird@ess.ucla.edu RE: Introduction to software in the "VERSATEC" directory THIS FILE CONTAINS A GENERAL INTRODUCTION TO THE SOFTWARE I AM OFFERING, WITH INSTRUCTIONS ON HOW TO OBTAIN RELATED PROGRAMS. THIS IS FOLLOWED BY DETAILED NOTES ON THE INDIVIDUAL FILES OF THE "VERSATEC" DIRECTORY. I ASSUME THAT YOU ALREADY HAVE "LARAMY". I currently offer four thin-plate finite-element programs: LARAMY finite-strain no faults flat-Earth crust/mantle FAULTS neotectonics with faults flat-Earth crust only PLATES neotectonics with faults flat-Earth crust/mantle SHELLS neotectonics with faults round-Earth crust/mantle The finite element program LARAMY and the associated graphics programs of the VERSATEC directory may be copied and used without charge or restriction, except that (1) any user who obtains publishable results is expected to acknowledge the source of the program, and (2) no charge may be made for copies of the software which are passed on to third parties. It was developed with support from the National Science Foundation and the University of California, which places it in the public domain. Also, I would be very pleased to see my 6 years of programming effort bring some benefit to others besides myself. Finally, I have no desire to be personally responsible for running simulations of every interesting problem for the rest of my career! The LARAMY program was created to permit simulations of continental deformation in which a crust of either one or two rheological layers (either one of which may enclose the brittle- ductile transition) overlies a mantle lithosphere and/or fluid asthenosphere and/or subducting slabs. While it has been customized for simulations of the Tertiary history of North America, the problem- specific data are mostly enclosed in DATA statements in two BLOCK DATA subprograms which can be replaced with relative ease. Its only fundamental limitations are that it is (1) a continuum program, with no fault discontinuities, but only an equivalent frictional plasticity; (2) rigorously isostatic, and neglects all bending stresses; and (3) written for a flat Earth. The VERSATEC directory contains programs needed to post-process the results of LARAMY runs to produce either black & white or colored contour maps of the variables on a Versatec electrostatic plotter. Before attempting to work with Laramy, it is essential to gain an overview of the program's structure, abilities, and limitations. Readers are STRONGLY ADVISED to skim the following two articles before proceeding further: 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. The first of these references contains figures which give some idea of the capabilities of the VERSATEC post-processing programs. (Actually, these images where created with the GRAFIC package and photographed off of a terminal; the VERSATEC package has the same capabilities, except that the backgrounds will be white.) The software needed to run the VERSATEC programs includes a FORTRAN-77 compiler (vectorizing, if possible), a linear-equation solver (which I obtained from IBM's ESSL), and the proprietary "Versatec Random" software supplied with the Versatec plotter. However, with a little time and ingenuity, the program can be adapted to work with any low-level graphics package, such as DISSPLA. About the compiler; it shouldn't matter which one you use, but I have had mysterious bugs under some versions of IBM VS FORTRAN prior to 2.4.0 (for the 3090/VF); version 2.4.0 seem to be OK. You should be aware of useful programs in other related subdiretories. Here is an overview. (Note: Numbers 0 - 5 tie to the Read_me.n files.) 0. Laramy2AI (added 1997) is a graphics program that takes output from LARAMY and produces single or multiple maps of various scalar, vector, or tensor variables-- as .AI files readable by Adobe Illustrator for Windows. This allows the maps to be previewed and edited(!) before printing, and also allows the printing to take place on a wide variety of hardware (including hardware run by service bureaus)! The .AI files produced can be read by Adobe Illustrator 4 for Windows 3.1, and by Adobe Illustrator 7 for Windows 95 or Windows NT. They probably can be read by Adobe Illustrator on Macs, too. The processing program Laramy2AI.exe runs under Windows 95 or Windows NT on Intel (or compatible) processors. Source code is also provided in Fortran 90 so that you can move this utility to Mac, Unix, or other systems. 2.GRAFIC is a two-part graphics package that post-processes finite element results to produce colored contour maps of all variables, with overlays to show vector and tensor directions and/or state lines, etc. The number-crunching part is run as a batch job, and produces graphics metabolism. The second part is run interactively from a graphics terminal to view the images at various scales, with or without the overlays, etc. These images can be photographed from the screen to make slides. This package requires IBM's Graphical Data Display Manager. 3. VERSATEC (this directory) is comparable to GRAFIC, but produces either B/W or color output on a Versatec electrostatic plotter (for publications). It runs as a single batch job and is not interactive at all. 4. FIXER is a fully interactive package for correcting small parts of a finite element grid which become excessively strained, so that a run can continue. It can also be used to fine-tune the files representing initial conditions. (However, it is not a bells-and-whistles grid-generator like some expensive engineering programs offer; just an editor. You must have already created and input a grid with the same topology as the final product you want.) A color graphics monitor attached to your mainframe is required, but not a mouse. You also need IBM's Graphical Data Display Manager. 5. MAPPER is a fully interactive package for making plate reconstructions of oceanic plates which are no longer on the Earth's surface (based on marine magnetic anomaly data and finite rotations) and saving the results in the form of DATA statements that can then be compiled into LARAMY in order to represent complex boundary conditions on the bottom of continents. This program also interactively displays the boundary conditions as they would be applied to a representative finite-element grid, so that unexpected or erroneous BC's can be corrected. A color graphics monitor attached to your mainframe is required, but not a mouse. You also need IBM's Graphical Data Display Manager. The hitch with accessory programs #2-#5 is that, while they have been written in pure FORTRAN77, they make thousands of calls to proprietary graphics packages which may not be available to you. My packages GRAFIC, FIXER, and MAPPER make calls to subroutines of the IBM Graphical Data Display Manager software (sometimes known as "that G*D-D*M* software"); my VERSATEC program calls the proprietary package that is supplied with that brand of plotter. My advice to you is that if you only want to plot maps of output (and you don't have a Versatec or IBM's Graphical Data Display Manager), then you should use Laramy2AI (in a different subdirectory) instead. If you need to run FIXER or MAPPER (and you don't have GDDM), it will probably be easiest and cheapest to hire a programmer to provide dummy interface subroutines that accept these calls and translate them into further calls to your local graphics software (rather than programming these tools from scratch). Or, if you have finite-element post-processors that you are satisfied with, another simple alternative is to write an "output translator" which reads the packed-numbers (device-9) output from LARAMY and puts it into the format required by your existing package. Please understand that I will expect you to communicate with me if you find significant bugs or opportunities for program improvement, and that if you use these programs to obtain publishable results I will expect the source of the program to be acknowledged in your manuscript. I also hope that you will remember me when you have reprints to distribute. =========================================================== NOTES ON THE FILES OF THE "VERSATEC" SUBDIRECTORY VERSCOMP = Fortran77 source code for the program DRAW, plus IBM Job Control Language statements needed to compile it into a load module. PLOT = IBM Job Control Language used to run the precompiled program DRAW (see above) as a batch job and send the plot to the Versatec plotter. BASEMAP = a file of digitized state lines in North America, useful as a geographic reference to overlay on displays; also a model for similar files that could be created for other regions. ----------------------------------------------------- VERSCOMP This file has two logical parts: the "outer wrapper" of IBM Job Control Language statements, and the program DRAW. JCL Statements. =============== The JCL statements all begin with '//'. If you wish to use these programs on a non-IBM system, you will begin by stripping these out and replacing them by commands to your own operating system. Their function is to compile the program, link it with the necessary graphics and matrix-algebra packages, and store it as a load module on disk. -IRAAGPB and EFF9GPB are two of my account numbers. Change them to yours everywhere(!) they occur. -I begin with a dummy job step, running a dummy program //CLEAR EXEC PGM=IEFBR14 //DD1 DD DSN=xxxxxxx.VERSMOD,DISP=(OLD,DELETE) The point of this odd procedure is actually in the following line, which instructs the system to scratch any existing load modules (from a previous compilation). Due to a bug in IBM operating systems, multiple uncatalogued datasets may be created unless you begin with this step. However, on the VERY FIRST run of VERSCOMP, these two lines: //CLEAR EXEC PGM=IEFBR14 //DD1 DD DSN=xxxxxxx.VERSMOD,DISP=(OLD,DELETE) should be omitted, because VERSMOD does not yet exist. -Check with your local consultants for the local name of the Fortran compiler. At UCLA this program FORTVS resides in a dataset called APP1.FORTVS.COMPILER, which has the logical name STEPLIB. However, the name of the dataset may be different at your institution. (I recommend version 2.4.0 or later if you use the IBM VS Fortran compiler; earlier versions have bugs which cause mysterious errors.) At the end of the compilation, the compiled program DRAW is saved in a temporary disk dataset (//SYSLIN DD DSN=&&LOADSET...). It is not yet ready to run because it lacks Fortran, algebraic, and graphical support routines. -At the end of the file, there is some further JCL. This defines the LKED (link-edit) step, where the compiled program DRAW is linked with necessary support routines to create load module VERSMOD(GO). Necessary Fortran77 routines are supplied under the name APP1.FORTVS.LIBRARY at UCLA; see your local consultants for the equivalent dataset name at your institution. -A second source of code needed during link-editing is IBM's Engineering & Scientific Subroutine Library, or some equivalent source of a linear-system solver. As described in the first page of DRAW, routines DGBF and DGBS are used to factor and solve real, asymmetric, banded linear systems. If you have ESSL at your institution (either the scalar or the vector version is OK), then identify the dataset containing its load modules in the place where I have "APP1.ESSLV". If you don't have ESSL, your work will be a little harder. Substitute appropriate calls to some local library. If your computer has 4-byte precision as the default, then use double-precision routines as I have. If your computer (a Cray?) has 8-byte numbers as the default, double precision is not needed and will be very costly. Unfortunately, different solvers assume different storage schemes for the banded coefficient matrix! To preserve portability, I have always referred to the coefficient matrix as a "vector" array, with a single subscript. The index to this vector (i.e., position in memory) is computed from the two logical (row, column) indeces by a statement function INDEXK which appears in many subroutines. Thus, if you change the banded-equation solver called by SOLVER, then you should also write a new statement function INDEXK and copy it into each subprogram that currently has such a statement (the name is always the same). [Incidentally, the incorrect compilation of DRAW by some earlier version of VS Fortran seemed to be connected with improper vectorization of loops which invoked this statement function for indirect addressing. ] -Link-editing also involves some routines from something called "APP1.VTEC21.FORT77". This is our local name for Fortran77-compatible load modules of the Versatec Random plotting package. If you have this at your institution, substitute the local name of this partitioned dataset of load modules. If you don't have it, or you can't convince your institution to get it, then the best alternative is to get only the Versatec Random manual, and then to write substitute Fortran subroutines that simulate Versatec Random functions by using some other graphics package you do have. (If you are adding Fortran routines to translate to another graphics package, this DD statement is the place to identify the location of that other dataset of load modules to the link-editor.) -The link-editor requires its own tiny dataset of one control statement, which appears as: //SYSIN DD * ENTRY MAIN /* I don't know why the link-editor cannot figure this out for itself, but it can't! Program DRAW. ============= Program DRAW resembles a stripped-down version of LARAMY, with the following main differences. While it reads the same parameter- input file as LARAMY, it reads further down into the file to find the plot-control parameters which specify which variables are to be plotted and how. It reads packed-numbers files on device 8 like LARAMY, but it does not produce any new ones. It has no subroutines to calculate the velocities of the nodes or to advance the variables through time, since LARAMY has already done this and stored the arrays in the dataset which will be connected to device 8. And, it has extra graphical subroutines called from the report-writing subroutine REPORT. The most important are PAINT (which does colored contour maps) and ETCH (which draws the finite element grids). The logical heart of the program is subroutine CONTEL which contours the variables within each finite element and then colors in the spaces between the contours. (Because the logic in CONTEL is very complex and was previously buggy, I have left in some debugging WRITE statements which are surrounded by comment lines ("C$$$$$$$$$$$$$$$$$$$$"). These WRITE statements are currently inactive (because DUMPIT=.FALSE.) but you can make use of them if you encounter any problems in the future.) If you wish to change the color palette used for contouring, this can be done by adjusting the DATA statements in subroutine PAINT. However, I don't recommend this. ----------------------------------------------------- PLOT This is a short file of IBM Job Control Language statements used to run the pre-compiled program DRAW and create hard-copy plots. If you are on an IBM MVS system, you can run this with only minor changes; if your system is another type, you may wish to ignore this file completely and start from scratch. However, you should still read the comments below about how the operation of DRAW is controlled by the input file (//FT05F001 DD DSN=xxxxxxx.IN9828). -IRAAGPB and EFF9GPB are two of my account numbers. Change them to yours everywhere(!) they occur. -The plot-control dataset is read from Fortran device 5 (FT05F001). It is actually file IN9828, which is available in the LARAMY directory. This single unified parameter file combines physical and graphics-control parameters. The graphics commands begin at about line 79, and are all fixed-format (so don't enter them in any other columns than the ones that they appear in presently). The first plot parameter is an integer which specifies which timestep is to be plotted (or, you can choose to plot all of them by entering 0). Following that are a list of 24 variables which you may choose to plot. For each variable: Use "T" (True) to plot it or "F" (False) to skip it. The next number following on the line is the desired contour interval (always positive). Do not choose a tiny interval or the program will run out of workspace and "give up on" many of the finite elements, leaving them blank (it will also get very expensive). If you leave this blank, the program will assign a contour interval based on the number of contours you typically prefer. (This is set in an input line further down in the file; try 10.) The third entry on the line, also optional, is a "normal" value of the variable which, if specified, will be "pinned" at a mid-range color (yellow or yellow-green) or grey-level in all plots. (This is to prevent one undesirable result of auto-scaling, in which the predominant color or shade of a plot changes from one timestep to the next, even though the modal value of the variable has not changed.) Finally, there is an integer switch, which is given a value of +1 or - 1. This determines which end of the color scale is associated with high values of the variable. A setting of +1 means that red and white shades are high, a setting of -1 means that blue and black shades are high. Specifying plot scale is a little confusing. If you have used cgs (centimeter units) in your input, then this is just the denominator of the ordinary dimensionless map scale (e.g., 10,000,000). But, if your input was in SI (meters) or feet or some other units, the scale becomes the number of these units that you want to fall into each centimeter of the plot. You can also specify several pen weights, which are expressed in number of pixels. If you do very large (wall-map) plots, you might wish to increase these. One of the best features of this program is the ability to switch from color to B/W with a single logical switch (toward the end of the parameter input file). If you don't want color, you get a set of customized grey-levels which are textured so that they should be distinguishable even after reproduction in journals!!! I have one piece of advice: It is very cumbersome to try to increase the size of the dots in these shading patterns, because you have to redefine the patterns as hex strings and recompile. Don't bother. Instead, reduce the size (increase the scale) of the whole plot until the existing dot sizes are relatively large enough for reproduction in your chosen journal. You will note that there are logical options to plot a basemap of state lines, and also to retrodeform this basemap for times before the present. If you do choose to plot state lines, then you will have to connect a dataset such as BASEMAP (provided) to Fortran device 11 with a DD statement, as I have done in this JCL file PLOT. An interesting option in this program is that the statelines may be retro-deformed to their previous positions if you desire. The way that this is done is the final (present-day) positions of the nodes of the finite-element grid is read from Fortran device 12. Then, the difference in the positions at any timestep currently being plotted is the paleo-displacement. This is interpolated to move all features of the basemap back to their paleo-positions. The result is very elegant; as one flips through a series of plots from past to present, the basemap slowly unflexes from a distorted to a familiar geometry. This option is useful because most geologists think in a reference frame that follows the outcrop, not an abstract one fixed to some distant ridid part of the plate. However, if a computation is not finished, the final node locations may not yet be available! Therefore, there is a logical switch "T" or "F" in the input file just after the switch which selects for inclusion of a basemap, and you can turn off the retrodeformation feature. If you are not using a basemap, then Fortran devices 11 and 12 do not need to be defined in PLOT. If you are not retro-deforming, Fortran device 11 needs to be connected, but not device 12. My preference is to leave the lines in PLOT as a reminder, but to "DUMMY" out the datasets I'm not using. Look at the current contents of PLOT to see an example of this. -Line printer output is produced on Fortran device 6 (//FT06F001 DD ...). This is in the same format as in LARAMY, and is sufficient for getting a quick look at the behaviour of the variables if the plotter is down. It also echoes the values of the plot-control parameters, so you should check that they are being read correctly if you have any troubles. -Our particular plotting package is designed to produce a lot of boring output on Fortran device 7 for each plot describing the plotter setup. Since this program makes a separate plot for each variable and each timestep, this can get voluminous. I have gotten rid of this by assigning FT07F001 to a temporary system dataset which is not kept. -Fortran device 8 (FT08F001) is the logical name of the packed- numbers file produced by LARAMY on its Fortran device 9. If you do not have one of these yet, you can also use an initial-conditions file such as CORD11N from the LARAMY package, as these are in the same format. (However, many of the variables will have constant values, producing very dull plots. Try plotting elevation, crustal thickness, or heat flow.) -Fortran device 11 (FT11F001) is the logical name of the dataset containing the basemap of state lines (if requested in your plots). It is assigned to the dataset BASEMAP (provided) which is describe below. -If you request state lines on your plots and if you want them retro-deformed, then you must supply the final positions of the nodes of the finite-element grid to Fortran device 12 (FT12F001). You can make such a file by copy the first section (node locations) of the last report produced by LARAMY on its device 9. I have not included such a file, however. ----------------------------------------------------- BASEMAP This is a very simple file containing digitised statelines in the 48 conterminous United States, plus the northern 2/3 of Mexico, plus the southern 2/3 of Canada. If you choose (by setting parameters in the Fortran device-5 input file IN9828) to have state lines included in plots, then this dataset is required when you use PLOT (provided) to run the compiled version of program DRAW (see VERSCOMP above). The state lines are expressed as a sequence of straight-line segments connecting points. The file consists of a string of point positions, and for each a logical indicator (T/F). If the indicator is T, then the "pen is down" when moving to the point, and a segment is drawn. If the indicator is F, then the "pen is up" and we are beginning a new line by moving to its initial point. You can easily write a little program to take digitizer output and put it into this format, for any other region where you need a basemap. The points are expressed as (latitude, longitude) coordinate pairs. The units are degrees. Latitude is positive in the Northern hemisphere. Longitude is positive East of Greenwich, England (the prime meridian). It is not important whether a point is indentified with a positive or a negative longitude (i.e., 245 = -115). ----------------------------------------------------- GOOD LUCK !!!!!!!!