In this example, we will simulate the evolution of a square fault in an elastic full space. The fault will be loaded plate motion and the relationship between slip rate and traction on the fault will be governed by rate and state friction.
First, let's import our tools!
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import tectosaur.mesh.mesh_gen
import tectosaur as tct
import tectosaur.qd as qd
qd.configure(
gpu_idx = 0, # Which GPU to use if there are multiple. Best to leave as 0.
fast_plot = True, # Let's make fast, inexpensive figures. Set to false for higher resolution plots with latex fonts.
)
plt.style.use('default')
Next, we'll create a square triangulated fault mesh:
In [ ]:
n_fault = 25
L = 1.0
W = 1.0
m = tectosaur.mesh.mesh_gen.make_rect(n_fault, n_fault, [
[-L, 0, W], [-L, 0, -W],
[L, 0, -W], [L, 0, W]
])
n_tris = m[1].shape[0]
Now, let's define all the relevant material and numerical parameters needed:
In [ ]:
qd_cfg = dict(
# Material properties
sm = 3e10, # Shear modulus (Pa)
pr = 0.25, # Poisson ratio
density = 2670, # rock density (kg/m^3)
# Frictional properties
Dc = 0.000002, # state evolution length scale (m)
f0 = 0.6, # baseline coefficient of friction
V0 = 1e-6, # when V = V0, f = f0, V is (m/s)
a = np.ones(n_tris * 3) * 0.010, # The rate-state a parameter
b = np.ones(n_tris * 3) * 0.015, # The rate-state b parameter
# Boundary conditions
plate_rate = 1e-9, # (m/s), equivalent to ~31.5 mm/yr
# This is only necessary because this is a full space model and there's no concept of depth or gravity
additional_normal_stress = 50e6,
# numerical preferences
timestep_tol = 1e-4, # error tolerance for the RK45 time stepper
tectosaur_cfg = dict(
quad_coincident_order = 6, # Order of coincident quadrature
quad_edgeadj_order = 6, # Order of edge adjacent element quadrature
quad_vertadj_order = 6, # Order of vertex adjacent element quadrature
quad_near_order = 5, # Order of nearfield, but not touching, element quadrature
quad_near_threshold = 2.5, # The distance threshold for defining near vs farfield
quad_far_order = 2, # Order of farfield quadrature
quad_mass_order = 4, # Order of gaussian quadrature used for the mass operator
float_type = np.float32, # 32 bit or 64 bit floating point?
use_fmm = False, # Should we use the FMM for the farfield?
fmm_mac = 2.5, # The FMM P2P vs M2P distance criteria
pts_per_cell = 100, # The maximum number of points per cell in the FMM tree?
fmm_order = 2, # The order of the FMM solid harmonic approximation
log_level = 'INFO' # The logging level to display, use DEBUG for lots of info, INFO for less, and WARNING for almost none.
)
)
The model object is the main implementation of the QD equations. We create a FullspaceModel. There is also a TopoModel for use when we want to have a free surface, possibly with non-planar topography.
In [ ]:
model = qd.FullspaceModel(m, qd_cfg)
There are two relevant length scales for QD simulations:
In [ ]:
qd.print_length_scales(model)
We will use an initial state variable value of 0.7 everywhere on the fault and an initial slip of 0. As a result, the initial stress state is zero everywhere.
In [ ]:
init_state = np.full(model.m.n_tris('fault') * 3, 0.7)
init_slip = np.zeros(model.m.n_tris('fault') * 9)
init_conditions = np.concatenate((init_slip, init_state))
Next, we create the Integrator object. This object encapsulates the logic that runs the Runge-Kutta time stepper and stores the results to a folder. By default, every 100 time steps, a file with the latest results will be output to a folder prefixed with data. If data0 already exists, data1 will be used and so on. To create an Integrator, we provide a model object that describes the frictional and elastic behavior of the model we're dealing with. We also provide the initial conditions, a big vector consisting of the all the slip degrees of freedom and then all the frictional state degrees of freedom. Because there are three basis functions per triangle and three components of slip and one component of state per basis function, we have 3 * 4 * n_tris total degrees of freedom.
In [ ]:
print(init_conditions.shape, 3 * 4 * n_tris)
The most important part of the Integrator object are the two lists h_t and h_y. Each element in h_y is an array with the same shape as init_conditions containing the value of each degree of freedom at that time step.
In [ ]:
integrator = qd.Integrator(model, (0, init_conditions))
Finally, we will start the simulation. We'll run for 1000 time steps and plot every 10th time step.
At each time step, many plots will be produced. The columns are the x, y and z components of each field. The rows are:
In [ ]:
from IPython.display import clear_output
def display_fnc(integrator):
print(len(integrator.h_t))
print(integrator.h_t[-1])
if len(integrator.h_t) % 10 == 0:
clear_output(wait = True)
t = integrator.h_t[-1]
y = integrator.h_y[-1]
integrator.model.display(t, y)
integrator.integrate(
n_steps = 1000,
display_interval = 1,
display_fnc = display_fnc
)
Finally, a quick description of how the data is stored. Everything from a given simulation is stored in the same data# folder (for example, data0). First there is a file, initial_data.npy that stores a tuple consisting of (mesh, qd_cfg, init_conditions). This is essentially all the information needed to restart the simulation from the beginning.
In [ ]:
import numpy as np
mesh, qd_cfg, init_conditions = np.load('data0/initial_data.npy')
print(mesh)
print(cfg)
print(init_conditions)
Then, the time history and slip/state history is stored each 100 time steps in a file named after the last time step (e.g. data0/300.npy).
In [ ]:
ts, ys = np.load('data0/300.npy')
print(len(ts))
print(len(ys), ys[0].shape)
print(ts[:5])
print(ys[:5])
That should be sufficient to get started analyze the results, but check out the fullspace_qd_plotter.ipynb notebook for some built-in tools to make nice figures and videos of the results.