Example Problem in 2D

To illustrate how to specify parameters in a text file, here is an example problem test2d.in (included in the problems directory). This example illustrates a simple 2D rupture problem based on the SCEC Rupture Code Verification Group TPV3 (this is a horizontal slice of the 3D simulation at hypocentral depth). The initial stress and friction parameters are homogeneous, with the exception of a nucleation patch at the center of the fault and strong frictional barriers at the external boundaries of the fault. The simulation saves several fields, both on-fault and off-fault.

[fdfault.problem]
test2d
data/
1000
0
0
0.3
50
4

[fdfault.domain]
2
2
801 802 1
1 2 1
801
401 401
1
1
slipweak
4
elastic

[fdfault.fields]
0. 0. 0. 0. 0. 0.
none
none
0

[fdfault.block000]
2.67 32.04 32.04
0. 0.
40. 20.
absorbing
absorbing
absorbing
none
none
none
none
none

[fdfault.block010]
2.67 32.04 32.04
0. 20.
40. 20.
absorbing
absorbing
none
absorbing
none
none
none
none

[fdfault.operator]
0.

[fdfault.interface0]
y
0 0 0
0 1 0

[fdfault.friction]
2
constant 0. 0. 0. 0. 0. -120. 70. 0.
boxcar 0. 20. 1.5 0. 0. 0. 11.6 0.
none

[fdfault.slipweak]
3
constant 0. 0. 0. 0. 0. 0.4 0.677 0.525 0. 0. 0.
boxcar 0. 2.5 2.5 0. 0. 0. 10000. 0. 0. 0. 0.
boxcar 0. 37.5 2.5 0. 0. 0. 10000. 0. 0. 0. 0.
none

[fdfault.outputlist]
vx
vx
0 1000 100
0 800 1
0 801 1
0 0 1
V
V
0 1000 10
0 800 1
401 401 1
0 0 1
sxy
sxy
0 1000 1
600 600 1
501 501 1
0 0 1


[fdfault.frontlist]
0

This model is fairly simple, so use of a text input file rather than a python script is a reasonable choice. The following shows the model geometry:

The input file first sets the problem variables. The simulation and output directories are named, then the time stepping information is set up. The simulation will take 1000 time steps, and use a CFL ratio of 0.3. Output will be written to the console every 50 time steps, and the integration order is 4th order.

Next, the simulation domain is set up. The problem will be a 2D mode 2 fracture. The total domain will have 801 grid points in the \({x}\)-direction (along strike) and 802 grid points in the \({y}\)-direction (across the fault). Note that the extra grid point in the \({y}\)-direction accounts for the shared grid point along the fault; this will maintain the same grid resolution in both directions. We will use 2 blocks in the \({y}\)-direction, so since this is a 2D problem the number of blocks in each cartesian direction is 1 2 1. Following this, the exact number of grid points in each block along each direction is entered. Since there is one block in the \({x}\)-direction, we must give one number (801) for the \({x}\) grid length. Since there are two blocks in the \({y}\)-direction, we give two numbers (401 401). The \({z}\)-direction just has a 1 to signify that there is only one grid point in that direction.

Next in the domain, we specify the interfaces. There is a single slip-weakening frictional interface in this simulation, so we give 1 for the number of interfaces and slipweak for the interface type. Finally, we set the finite difference order to be 4, and set the material properties to elastic (this is needed to correctly set up the blocks, fields, and simulation output).

Information on the fields comes in the following section. This contains information on data stored at every grid point: the initial stress fields, heterogeneous material properties (located here because of how the code stores the heterogeneous field values), and whether or not the code will save the full plastic strain tensor for plasticity problems. Because this particular simulation chooses not to use any of these capabilities, all of these are set to zero, none, or false. For more information on how these are used, see the sections on code input.

The simulation blocks are set up in the next section. The simulation has two blocks: block000 and block010. The three digit indicies are the x/y/z coordinates of each block, so 000 is the first block in the \({y}\) direction, and 010 is the second block in the \({y}\) direction. This sets the material properties, geometry, and boundary conditions on each of the blocks. The blocks both have the same material properties: \({\rho = 2.67}\) MPa s^2 / km / m (note the funny units, which are used to allow slips to be measured in meters while lengths are measured in kilometers, and the elastic moduli to be measured in GPa while the stresses are in MPa for calculational convenience), \({\lambda = 32.04}\) GPa, and \({G = 23.04}\) GPa. The blocks are each 40 km x 20 km, and they align at \({y = 20}\) km on the fault plane. All external boundaries (the input order is left, right, front, back) are absorbing, while the shared boundary at the fault has boundary conditions none (back for the 000 block, front for the 010 block), since the friction law will set these boundary conditions. This problem does not use any non-rectangular grids, so the section where filenames would be given for external curves/surfaces all give values of none.

This simulation does not use artificial dissipation, though it includes the operator section for completeness.

The three lists that follow describe the fault friction. The problem has a single interface (interface0), with a normal direction in computational space in the \({y}\)-direction, and it joins together blocks 000 and 010. All of this is specified in the interface section. Next come the information on interface stresses and friction parameters, which are in the friction and slipweak sections. Note that no number is required for these section – the code simply finds the first occurrence of friction and slipweak following the interface0 section.

Under the friction section, the traction on the fault is modified with two “perturbations.” The first sets the constant initial background stress to be \({\sigma_n = -120}\) MPa and \({\tau_h = 70}\) MPa. Since this is a 2D problem, the vertical interface traction is ignored. Next is a boxcar function to nucleate the rupture by overstressing a fault patch. The boxcar is centered at \({x=20}\) km and has a half width of 1.5 km, and the patch includes an additional 11.6 MPa of horizontal traction (this ensures that the initial stress, which is the sum of all perturbations, exceeds the static friction level). No file is used to further modify the initial tractions, so none is given for the load file entry.

The friction parameters are changed in the slipweak section using three perturbations. The first entry sets the constant background values to \({d_c = 0.4}\) m, \({\mu_s = 0.677}\), and \({\mu_d = 0.525}\). The simulation has no cohesion or forced rupture. The other two perturbations increase the initial static friction to be high enough to be unbreakable barriers for the outer 5 km at each end of the domain; these prevent the rupture from hitting the edge of the simulation domain. No parameter files is used to further modify the friction parameters.

The final sections set up simulation output. The simulation saves three output units to disk. These output units illustrate the different ways that output can be saved to disk.

The first ouput unit saves the full 2D grid values of \({v_x}\) (the horizontal particle velocity) every 100 time steps. The first entry gives a name for the output unit (used to keep track of the different output units, so each one should have a unique name), and the second entry is the field to be saved (vx here, to denote the particle velocity, and note that this is case sensitive). Next are three numbers describing time output: start end stride. The simulation will start saving data at the 0 time step, stop saving at the 1000 time step (i.e. the end of the simulation), and save every 100 time steps (chosen to be large as the files can be quite large if we save every grid point). The next three sets of numbers are the start, end, and stride for the x, y, and z directions. Thus, we start at grid point 0, end at grid point 800, and take every point for x; start at 0, end at 801, and take every point for y; and start and end at 0 for z (since this is a 2D simulation).

The second output unit saves the slip velocity on the fault. The name and field are both V, and since the slip velocity is only defined on the fault, you will get an error if you try to choose values that are not on the fault interface. Here I pick values on the fault using y indices of 401 for both the start and end values. The x indices go from 0 to 800, and the time values go from 0 to 1000, saving every 10 points.

Finally, we save the shear stress at a single grid point for all time steps. The field sxy indicates which stress tensor value we want, and the time values go from 0 to 1000, saving every time step. Note that the spatial indices have the same starting and ending values, ensuring that only a single grid point is saved.

The output list must end with a blank line before any additional sections may be specified on the input file. The only one that comes afterwards is the frontlist, which determines the time that all grid points on the fault rupture. This is turned off for this simulation.

This shows all of the different sections that need to be set up to specify a rupture problem. This is a simple 2D simulation with homoeneous properties (save the nucleation point and the barriers to rupture), so it can easily be written as a single text file. More complicated problems are better handled using the Python module, which is described in the next example.