help.txt

main: binary: interact
main: internal double prec, A matrix double prec
main: initializing on 2023-11-27 23:08:18
init: compiled on Nov 27 2023 22:36:41, running in serial

Interact: calculate fault stresses and displacements in a half space or in 2-D
          using a boundary element approach.

Program reads in patch dividing a fault geometry and either solves
for stresses and/or displacements in a one-step problem for various boundary
conditions (sign (i.e. direction) constrained or unconstrained), or for a simulated loading
cycle where each patch follows a frictional constitutive law (e.g. Navier-Coulomb)
and repetitive rupture is allowed during continuous "plate-tectonic" loading.

When stresses (stress tensor components) or pressure are referred to, positive values mean
extensional and compressional regimes, respectively (physics convention).

Output is written to files and in real-time to X11 if PGPLOT support was compiled in.
PGPlot support was compiled in. 

PetSc support was compiled in, providing limited access to parallel solves for one-step problems.
      For direct solve, use    "-pc_factor_mat_solver_type scalapack -mat_type scalapack" or
                           "-pc_factor_mat_solver_type elemental -mat_type elemental".
      For iterative solve "-ksp_type fgmres -pc_type none   -ksp_max_it 10000 -ksp_rtol 1.0e-8" or
                          "-ksp_type fgmres -pc_type jacobi -ksp_max_it 10000 -ksp_rtol 1.0e-8".
      Check the makefile for other solver options and MPI settings.
      When running a one-step computation, will also compute stress fields in parallel.

(1) The fault geometry is input via the file "geom.in",
    a list of fault patches.
    This file has the following ("patch") format for rectangular faults or point sources
    (free format ASCII list of parameters):

    x_0 y_0 z_0 strike_0 dip_0 hlength_0 hwidth_0 group_0 ...
    x_1 y_1 z_1 strike_1 dip_1 hlength_1 hwidth_1 group_1 ...
    ...
    x_N-1 y_N-1 z_N-1 strike_N-1 dip_N-1 hlength_N-1 hwidth_N-1 group_N-1 ...

    Here, N is the total number of patches and

    - x,y,z are the coordinates of the center of the rectangular fault patch or
      in the case of a point source, the point source location.

    - strike and dip are fault plane orientation angles in degree
      strike is defined in degrees clockwise from North
      dip    is defined in degrees downward from the horizontal, 90 degrees means vertical fault
      note that right now, the geometrical rake has to be 90 or 0 degrees, i.e. no tilted rectangles
      slip can, however, have a rake as given by different along-strike and along-dip values, see below

    - hlength (L) and hwidth (W) are the patch's -->HALF<-- length and half width in 
      strike and dip direction, respectively. The patch area is thus 4LW.

    - group is the patch's fault group, a number that is used to define faults that consist of
      several patches (discretizing fault planes by rectangular elements) (0 <= group <= N-1))
    - ... at the end mean possible additional input, see below.

    Since the "ALLOW_NON_3DQUAD_GEOM" flag was used during runtime, the program can deal with
    some additional patch geometries besides rectangular. Those are selected as follows:

    - point source in half-space
      If ONLY fault half-length is set to a negative value, a point source will be assumed instead
      of a rectangular fault patch. In this case, width should be set to the `fault' area
      and half-length is equivalent to -aspect_ratio, where aspect_ratio is some equivalent L/W.

      WARNING: Not properly tested yet.

    - triangular in half-space:
      If BOTH fault half-width and length are negative, then the patch is a triangular element.
      In this case, x, y, z, have no meaning but will be reassigned from the centroid.
      strike and dip are global projection directions, and length and width will be sqrt(area) subsequently.
      The next nine numbers in the input geometry line will be the coordinates of the three nodes of the
      triangle, however. FE-style counterclockwise numbering, input is then in the format:

      999 999 999 999 999 -1 -1 group_0 x_x^1 x_y^1 x_z^1 x_x^2 x_y^2 x_z^2 x_x^3 x_y^3 x_z^3

      where exponents indicate the local number of the node, and 999 are place holder values, not used.

      WARNING: This is experimental and not fully tested.

    - segments in two dimensions in the x-y plane
      If fault half-width is zero, then dip should be 90, z=0, and program will run in 2-D mode.
      In this case, all z coordinates should still be specified (e.g. in the grid output requests)
      but z should be set to zero. Computation is plane-strain unless changed by -ps.
      For plane stress, all output values of u[Z] (displacements into stress free direction) are
      given as strains e_zz.
      Calculations are performed in a full plane by default, can be changed with the -hp option.

    If all your patches are rectangular in 3-D, you could recompile
    without the ALLOW_NON_3DQUAD_GEOM set to save memory and enhance speed.


    Examples:

    0 0 -1 30 90 4 2 0

    for a rectangular, vertical fault patch centered at depth -1, striking 30 degrees from North
    with length 8 and width 4, patch group number is 0.

    0 0 -1 30 90 -1 32 0

    is the same kind of fault, only as a point source. Note that 32 is the fault area from above.

    999 999 999 999 999 -1 -1 0   2.0 3.4641 0  -2.0 -3.4641 0  0 0 -4

    would be for a triangular element with roughly the same geometry.

    Please note that there are various programs to convert to and from the above patch format,
    e.g. patch2xyz to go to GMT style xyz coordinates for the edges of each patch,
    patch2vtk to go from patch format to (paraview) VTK format, or
    points2patch to convert a set of four points in FE ordering to the above patch format.


(2) Boundary conditions and the operational mode are input via the "bc.in" file
    which has the following format (free format ASCII list of parameters):

    operational_mode:
      1: one step calculation
      2: loading simulation, and
      3: loading simulation and x window plotting.


 IF ONE STEP CALCULATION IS CHOSEN, the following input fields should be

     print_bulk_fields

      If print_bulk_fields is set to 0, no stress or displacement fields will be output.
      If print_bulk_fields is set to 1, the program will read in the boundaries from the next line

       xmin xmax n ymin ymax m zmin zmax o

       where n, m, and o are the number of samples between the given limits.

         If o is set to -1, the output will be in a plane with the
          average strike and dip direction of fault group 0. x/y min/max are then the
          limits in strike and dip direction.
         If o is set to -2, the output will be in a plane with the
          average strike and normal direction of fault group 0. x/y min/max are then the
          limits in strike and normal direction.
         In these cases, the coordinates in the output files will correspond to the local,
         rotated systems, e.g. for -2, x will go along strike, and y along the normal direction
         with respect to the average fault patch geometry, with origin in the center of the fault
         group. HOWEVER, the displacements and stresses will still be given in the old system!
         If the -pc switch is set, the global coordinates will be printed to "plane.xyz"

      If print_bulk_fields is set to 2, the program will read the output locations from file "oloc.dat"
      This file has the x, y, and z coordinates in ASCII as an unformatted list.


     patch_number boundary_code boundary_value

      patch_number runs from 0 to N-1, where N is determined from the number of patches
      (fully read lines, unformatted) in the geometry file.

      If patch_number < 0, then the line should instead read

       -increment start_fault stop_fault boundary_code boundary_value

      and boundary code boundary_code will be assigned with value
      boundary_value from start_fault to stop_fault with increments increment.

      If start_fault < 0 and stop_fault < 0, will select all patches.


     If the patches for which boundary conditions are specified are triangular,
     then the strike and dip specified in the geometry input will define global directions
     the prescribed slipor stress values refer. The slip boundary conditions will then be converted to the triangle local
     reference frame internally, and converted back to the global frame on output for fault properties.

     The patch_number ... boundary_value line can be repeated as often as necessary (say, twice for each fault
     if the strike and dip movement modes are to be activated the same time, but not for triangles yet, for 
     which specifying either strike or dip stress will enforce the other shear component to be zero.)

   Possible boundary conditions (as indicated by their integer codes) are:

   (A) 0, 1, 2: slip on fault specified (if called several times, will add up)
    0 means strike, 1 dip, and 2 normal direction

    For strike: positive values of slip mean left-lateral fault, negative right-lateral;
    	this corresponds to a resulting drop resp. increase in the strike component of stress.
    	(We use stress values with signs, that is a negative stress drop corresponds to a
    	reduction of positive shear stress, a positive stress drop to reduction of shear
    	stress with negative sign. Both times, the absolute value of the shear stress is
    	reduced).
    For dip:    positive values of slip mean thrust fault, negative normal fault;
    	this corresponds to a resulting drop resp. increase in the dip component of stress.
    For normal: positive values of slip mean explosive source, negative implosive;
    	this corresponds to a resulting drop resp. increase in the normal component of stress.

    For triangles, these are all referring to the global frame specified.

   (B) 10, 11, 12 or 20, 21, 22 or 30, 31, 32: stresses are specified

    where the following codes indicate the degree of freedom for the faults to
     achieve that stress starting from 0
    10	strike slip, either way
    11	strike slip, left lateral   (stress drop bc usually < 0 for activation)
    12	strike slip, right lateral  (stress drop bc usually > 0 for activation)
    20	dip direction, either way
    21	dip direction, thrust       (stress drop bc usually < 0 for activation)
    22	dip direction, normal       (stress drop bc usually > 0 for activation)
    30	normal direction, either way
    31	normal direction, explosive (stress drop bc usually < 0 for activation)
    32	normal direction, implosive (stress drop bc usually > 0 for activation)

    If any of the strike or dip modes (10, 11, 12 or 20, 21, or 22) are given as 
    (110, 111, 112 or 120, 121, or 122) instead (add 100) this means
    that the target stress should be corrected by the change in normal
    stress times the dynamic coefficient of friction (makes sense for Coulomb type rupture).

   (C) 50, 51, 52: read in background loading (see -l, NOT -s), reduce resolved stresses on faults

    50	such that faults are shear-stress free after slip in strike and dip directions
    51	or such that the shear stress is given by the static friction coefficient
      	times the normal stress, after slip in strike and dip direction
    52	or such that faults are shear stress free after slip in strike direction only
    For the above two cases, the BC value (the last number) is meaningless.


    NOTE: The final stresses, e.g. in flt.dat, after slip will correspond to the stress drops
    needed to reduce the background load as specified above. The stresses will NOT include
    the resolved background stresses since those are only used to initialize the problem
    as boundary conditions. I.e. add the resolved background shear stress to the shear stress
    specified in flt.dat stress to get shear stress free conditions for the 50 case, for example.
    For debugging, the resolved shear, dip, and normal stresses for each patch will be written to "rstress.dat".

    Mixed boundary conditions: If you want to mix specified slip boundary conditions with stress
    boundary conditions, slip has to be completely specified first such that subsequent stress conditions
    can take the stress changes induced by the slip into account.

    WARNING: Stress boundary conditions will NOT work for point sources, and are
    unreliable for triangular sources right now.

    Example:
    1 1
    -5 5 51 -5 5 51 -1 -1 1
    -1 -1 -1 10 1


 IF LOADING SIMULATION IS CHOSEN (with or without X output)

   the first three numbers are:

    time_step print_interval stop_time

   where time runs from 0 to stop_time in base time steps of time_step
   and the program will print to fault files every print_interval intervals.

   The time_step will normally be adjusted (but see -ct) based on the
   program's prediction of the next rupture time or change in shear stress sign.
   The max step is given by print_interval to ensure frequent
   enough output to the flt.*.dat files.
   The minimum time step was set to 1.42857e-16.
   Time has to monotonously increase from zero.

   If X windows output is chosen, the next line (two numbers) should read

    x_plot_interval x_scroll_interval

   where x_plot_interval and x_scroll_intervals give the timing for the map view
   and time vs. x plots, respectively.

   All following lines assign active modes to faults and are of format:

    patch_number activation_mode_code

   where fault numbers should be given from 0 to N-1, N being the number of patches.
   If patch_number < 0 then give

    -increment start_fault stop_fault activation_mode_code

   instead. If start_fault<0 and stop_fault<0 will select all patches, as above.

   Activation_mode_codes are:
      0		inactive(default)
     10		strike slip,   either way
     110	strike slip,   Coulomb correction, either way
     11		strike slip,   left lateral
     111	strike slip,   Coulomb correction, left lateral
     12		strike slip,   right lateral
     112	strike slip,   Coulomb correction, right lateral
     20		dip direction, either way
     120	dip direction, Coulomb correction, either way
     21		dip direction, thrust
     121	dip direction, Coulomb correction, thrust
     22		dip direction, normal
     122	dip direction, Coulomb correction, normal
     30		normal direction, either way (not implemented yet)
     31		normal direction, explosive  (not implemented yet)
     32		normal direction, implosive  (not implemented yet)
     40		maximum stress direction in strike/dip plane
     140	maximum stress direction in strike/dip plane, Coulomb correction

   If Coulomb correction is set, the modes that involve slip in strike and dip direction will
   result in a correction for normal stress changes, assuming that the target
   stress drop is modified by the product of normal stress change during slip and
   the dynamic friction coefficient.


 NOTICE:
 Program maintains normal (opening) modes. If this is not needed, recompile
 with the NO_OPENING_MODES switch set to save memory.


PARAMETERS:

 Elastic and frictional law properties are hard-coded and should be set in the "properties.h" file.
 Whenever you want to change those, you will have to recompile the program. This was done to improve
 the execution speed, assuming that material properties will not change between different model runs.

 In the following, we will refer to two different matrices, A and I. The interaction, or I, matrix is
 assembled once for the loading simulations and holds all possible stress influence coefficients for all
 combinations of patch slips for a given geometry. I is typically held in memory and can be quite big, depending
 on the problem. Hence, there is the option to store a reduced, sparse version of I. If patches are activated
 (slipping) during a loading simulation (according to their frictional law), a system of equations
 A x = b is set up to solve for the slip vector x that corresponds to the target stress drop. This A matrix
 is also the only matrix that is assembled for a one-step calculation as we specifically know the patch/boundary
 conditions we have to consider in this case. There are several possible modes for solving completely
 unconstrained A x = b systems; constrained systems (where some entries in x have to be >=0 for unidirectionally
 constrained slip) are always solved with a direct implementation of Lawson and Hanson's solver.

 The following parameters can be set to non-default values on the command line:

 -c  value      sets the cohesion term of the friction law to value, by default: 5
                If the cohesion=0 at all times, recompile with NO_COHESION set
                to improve speed. (Stress drop is set to 1, which is 0.0001 times mu).

 -pr value      set the background pressure to value "value".
                The background pressure amplitude is 10 by default, both for one-step and loading
                simulations. The pressure is constant with depth (see properties.h).

 -ms value      set minimum stress drop to value "value". (Characteristic stress drops are
                set to 1.) The minimum stress drop is zero by default, ie.
                allowing stress drops of all sizes (loading simulation).

 -ei value      set matrix cutoff value for sparse storage of I or A in absolute values
                             (5.000000e-05 by default.)
                (only used if I matrix is stored as sparse, see -nd option, or for sparse Ax=b solver.)

 -wc value      will use value as the maximum fraction of singular value for an SVD solution (default: 1.00e-15)



OPTIONS:

 The following operational modes (switches) can be switched on or off by
 using the toggle options listed below (type, eg, interact -s)
 Please read carefully which switches are ON by default. If those
 active switches are selected, they will be switched OFF instead.

 -s  read in initial stress values for each fault from file "fsi.in", format:

     s_strike^0 s_dip^0 s_normal^0
     s_strike^1 s_dip^1 s_normal^1
     ...

     (repeated for each patch). Here, the s_i are the initial values of the stress 
     tensor resolved on the fault slip directions. All values are ADDED to the
     homogeneous background stress as set in the program, ie. use all zeros for normal
     initial stress state. (loading simulation)
 --> WARNING: This option will only work for the loading simulation mode, for one-step calculations
     specify fault stresses as individual boundary conditions, or use -l to get a resolved background
     stress.

     ON by default, if switch is set will be OFF.

 -f  read in fault properties from file "fp.in", so far only parameters for the friction law
     are supported. Format:

     mu_static^0 mu_dynamic^0
     mu_static^1 mu_dynamic^1
     ...

     (repeated for each fault). Static and dynamic coefficients of friction for
     Coulomb.
     ON by default, if switch is set will be OFF.

 -i  suppress all interactions between fault patches except the self-effect of slip on any
     patches within the group that the patch belongs to and the effect of slip
     on the number 0 master group (this should be the main fault).
     (Meant as a experimental feature for LOADING SIMULATIONS ONLY!)
     OFF by default, if switch is set will be ON. See -ni for suppressed interactions in one step.

 -ni suppress all interactions between faults except the interactions of fault patches

 -ct use constant time steps for the loading simulation
     OFF by default, if switch is set will be ON.

 -b  suppress output of all fault stress and slip data. Normally,
     individual output files for each
     group of faults are generated when the number of groups is smaller than 50.
     If -b is set, this maximum number will be set to zero.

 -p  prints information about the cumulative slip of every patch in a group to the
     slipline.*.dat files if the fault files are used (see -b). (loading simulation)
     OFF by default, if switch is set will be ON.

 -w  whole fault activation mode (loading simulation)
     If set, one patch will activate all patches in the group for sliding.
     Also, if one patch is switched off because of low normal stress, all
     patches in the group will be switched off as well (*unless* -v is set).
     OFF by default, if switch is set will be ON.

 -k  keep slipping mode (loading simulation)
     If set, the first critical patches can trigger slip on other patches.
     If so, the initial patches will keep slipping, and the solution is obtained
     at each iteration until no more triggering takes place. If not set, slip
     on triggering patches is exhaustive, and iteration checks for triggering step
     by step (faster, since old solution doesn't have to be removed from stress
     field). The exhaustive method is affected by tiny differences in critical stress!
     If keep slipping, critical stress is -1.000000e-15, else -1.000000e-15.
     ON by default, if switch is set will be OFF.

 -v  whole fault deactivation mode (loading simulation)
     If set (and -w is *not* set), the deactivation (low compressive normal stress) of a patch
     will lead to a deactivation of the whole fault group.
     If -w is on (whole fault activation mode), this is the default. If you use -v in this
     case, the whole fault de-activations will be suppressed.
     OFF by default, if switch is set will be ON.

 -l  reads a/b factors for the linear background loading stress relationship.
     If set, the file "smlin.in" is read, it should hold 6 pairs of a_i b_i factors
     where the *non-hydrostatic* background stress at time t is given by s_i=a_i+b_i*t
     s_i = {s_xx, s_xy, s_xz, s_yy, s_yz, s_zz}, so that the format is

     a_xx b_xx
     a_xy b_xy
     ...
     a_zz b_zz

     To these factors, a depth independent hydrostatic pressure of 10
     will be added. You can change the amplitude of the pressure (e.g. to zero) with -pr.

     If no factors are read in (if "smlin.in" is not found), the code will use simple shear loading 
     with b_xy=1 plus the depth independent hydrostatic pressure of 10.

     (Modify pressure amp. with -pr). If you want to specify resolved background
      stresses in a one-step calculation, simply use the a_i factors to specify the stress 
      at time t=0 but mind the hydrostatic part.

     ON by default, if switch is set will be OFF.

 -nd switches to sparse matrix I storage, no matter what the I matrix size is.
     By default, the program goes from dense to sparse only if the
     matrix is bigger than 8000 MB.
     The sparse matrix storage implies that matrix elements whose absolute
     value is less than 5.000e-05
     are discarded. This behavior can be changed by specifying
     a different cutoff value using the -ei parameter option.
     (Make sure to use -si if you want to keep the matrix.)
     OFF by default, if switch is set will be ON.
     note that sparse here refers to only the storage mode of the I matrix,
     not the solution method of A.x=b. (see -ss).

 -ci checks the interaction matrix for positive Coulomb stress feedback loops, once calculated
     This may be useful for loading simulations.
     OFF by default, if switch is set will be ON.

 -si saves the interaction (I) matrix during a loading simulation,
     if I is written to a file during runtime and not held in memory.
     Else, "/tmp/i" and "/tmp/i.hdr" will be deleted.
     OFF by default, if switch is set will be ON.

 -oi reads old interaction (I) matrix from files "/tmp/i.dat" and "/tmp/i.hdr".
     Saves time when the geometry
     does not change between subsequent model runs (this is not automatically checked!)
     (Make sure to use -si if you want to keep the matrix.)
     OFF by default, if switch is set will be ON.

 -sa saves the A x = b A matrix to "a.dat" and "a.dat.hdr" during a one-step calculation.
     Binary format, with dimensions (n x m) and precision (p, in bytes) in the header file as n p m
     OFF by default, if switch is set will be ON.

 -oa reads A matrix from files "a.dat" and "a.dat.hdr" during a one-step calculation.
     Saves time when A
     does not change between subsequent model runs (this is not automatically checked!)
     (Make sure to use -sa if you want to keep the matrix.)
     NOTE that A is only constant for a constant geometry and the boundary conditions are
     either pure slip or pure stress components.
     If you are specifying a Coulomb stress change, elements of A will be modified
     by the boundary conditions, and not only b as one might think!
     OFF by default, if switch is set will be ON.

 -r  attempts to restart a loading simulation model run based on previous results
     WARNING: all result files will be overwritten NOT appended, so save before.
     For this to work, "geom.in" and "bc.in" (and possibly "fsi.in" and "fp.in"
     should not be substantially changed.

     Slip events from the previous run in the format as in the output file "events.bin" will be read in
     from the restart event file "restart.events.bin"

 -sn will print NaN values in displacement and stress outputs of the one-step calculation.
     Normally, those locations (end patch/segments) yielding one or more inf values in the
     stress matrix or the displacement vector will be omitted during output.
     OFF by default, if switch is set will be ON.

 -pc will print the global coordinates if fault-local planes are chosen for bulk stress
     and displacement output (-1 or -2 for o in zmin zmax o), writes xyz to "plane.xyz"
     ON by default, if switch is set will be OFF.

 -sv uses SVD solver for the unconstrained A x = b system.
     This approach should reliably find the least-squares minimum norm solution
     if the A matrix is singular.
     Note that LU can be orders of magnitude faster than SVD but use caution.
     If non-negative solutions are sought (as for specified slip directions),
     the program will automatically default back to the NNLS solver.

 -sl uses LU solver for the unconstrained A x = b system (default).
     LU can be orders of magnitude faster than SVD, but use caution.
     If non-negative solutions are sought (as for specified slip directions),
     the program will automatically default back to the NNLS solver.

 -d  debug output, if you want extra checks, compile with -DDEBUG
     OFF by default, if switch is set will be ON.

 -fpetsc    force Petsc solvers even for serial runs (else LAPACK for LU)

 -ps Change the default 2-D elastic approximation for segments from plane-strain to plane-stress.
     OFF by default, if switch is set will be ON.

 -hp run a 2-D plane strain calculation in a half-plane (y<=0) instead of plane
     OFF by default, if switch is set will be ON.

 -h  prints out this help message and exits to the operating system

(C) Thorsten Becker, thwbecker@post.harvard.edu, 1999 - 2023)
    interact - boundary element code for elastic half-spaces
    3-D quad dislocationS based on Okada (BSSA, 1992).
    3-D triangular dislocations based on Nikkhoo and Walter (GJI, 2015).
    2-D segment slip solution from Crouch and Starfield (1973).
    May include routines based on copyrighted software of others.
    Distributed under the GNU public license, see "COPYING".