In [1]:
# Make plots interactive in the notebook
%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
import os
import sys
# Add QDYN source directory to PATH
# Go up in the directory tree
upup = [os.pardir]*2
qdyn_dir = os.path.join(*upup)
# Get QDYN src directory
src_dir = os.path.abspath(
os.path.join(
os.path.join(os.path.abspath(""), qdyn_dir), "src")
)
# Append src directory to Python path
sys.path.append(src_dir)
# Get QDYN plotting library directory
plot_dir = os.path.abspath(
os.path.join(
os.path.join(os.path.abspath(""), qdyn_dir), "utils", "post_processing")
)
# Append plotting library directory to Python path
sys.path.append(plot_dir)
# Import QDYN wrapper and plotting library
from pyqdyn import qdyn
import plot_functions as qdyn_plot
To prepare a simulation, the global simulation and mesh parameters will have to be specified. This is done in three steps:
In this simulation, the only heterogeneity stems from a lateral variation in the direct effect parameter $a$, which is chosen such that the asperity has $(a-b) < 0$, and such that the matrix has $(a - b) > 0$.
In [2]:
# Instantiate the QDYN class object
p = qdyn()
# Predefine parameters
t_yr = 3600 * 24 * 365.0 # seconds per year
Lasp = 7 # Length of asperity / nucleation length
L = 5 # Length of fault / nucleation length
ab_ratio = 0.8 # a/b of asperity
cab_ratio = 1 - ab_ratio
resolution = 7 # Mesh resolution / process zone width
# Get the settings dict
set_dict = p.set_dict
""" Step 1: Define simulation/mesh parameters """
# Global simulation parameters
set_dict["MESHDIM"] = 1 # Simulation dimensionality (1D fault in 2D medium)
set_dict["FINITE"] = 0 # Periodic fault
set_dict["TMAX"] = 15*t_yr # Maximum simulation time [s]
set_dict["NTOUT"] = 100 # Save output every N steps
set_dict["NXOUT"] = 2 # Snapshot resolution (every N elements)
set_dict["V_PL"] = 1e-9 # Plate velocity
set_dict["MU"] = 3e10 # Shear modulus
set_dict["W"] = 50e3 # Loading distance [m]
set_dict["SIGMA"] = 1e8 # Effective normal stress [Pa]
set_dict["ACC"] = 1e-7 # Solver accuracy
set_dict["SOLVER"] = 2 # Solver type (Runge-Kutta)
# Setting some (default) RSF parameter values
set_dict["SET_DICT_RSF"]["A"] = 0.9e-2 # Direct effect (will be overwritten later)
set_dict["SET_DICT_RSF"]["B"] = 1e-2 # Evolution effect
set_dict["SET_DICT_RSF"]["DC"] = 4e-4 # Characteristic slip distance
set_dict["SET_DICT_RSF"]["V_SS"] = set_dict["V_PL"] # Reference velocity [m/s]
set_dict["SET_DICT_RSF"]["TH_0"] = set_dict["SET_DICT_RSF"]["DC"] / set_dict["V_PL"] # Initial state [s]
# Compute relevant length scales:
# Process zone width [m]
Lb = set_dict["MU"] * set_dict["SET_DICT_RSF"]["DC"] / (set_dict["SET_DICT_RSF"]["B"] * set_dict["SIGMA"])
# Nucleation length [m]
Lc = Lb / cab_ratio
# Length of asperity [m]
Lasp *= Lc
# Fault length [m]
L *= Lasp
# Find next power of two for number of mesh elements
N = int(np.power(2, np.ceil(np.log2(resolution * L / Lb))))
# Spatial coordinate for mesh
x = np.linspace(-L/2, L/2, N, dtype=float)
# Set mesh size and fault length
set_dict["N"] = N
set_dict["L"] = L
# Set time series output node to the middle of the fault
set_dict["IC"] = N // 2
""" Step 2: Set (default) parameter values and generate mesh """
p.settings(set_dict)
p.render_mesh()
""" Step 3: override default mesh values """
# Distribute direct effect a over mesh according to some arbitrary function
p.mesh_dict["A"] = set_dict["SET_DICT_RSF"]["B"] * (1 + cab_ratio*(1 - 2*np.exp(-(2*x/Lasp)**6)))
# Write input to qdyn.in
p.write_input()
Out[2]:
To see the effect of setting a heterogeneous value of a over the mesh, we can plot $(a-b)$ versus position on the fault:
In [3]:
plt.clf()
plt.plot(x, p.mesh_dict["A"] - p.mesh_dict["B"])
plt.axhline(0, ls=":", c="k")
plt.xlabel("position [m]")
plt.ylabel("(a-b) [-]")
plt.tight_layout()
plt.show()
As desired, the asperity is defined by $(a-b) < 0$, embedded in a stable matrix with $(a-b) > 0$.
The p.write() command writes a qdyn.in file to the current working directory, which is read by QDYN at the start of the simulation. To do this, call p.run(). Note that in this notebook, the screen output (stdout) is captured by the console, so you won't see any output here.
In [4]:
p.run()
Out[4]:
During the simulation, output is flushed to disk every NTOUT time steps. This output can be reloaded without re-running the simulation, so you only have to call p.run() again if you made any changes to the input parameters. To read/process the output, call:
In [5]:
p.read_output()
Out[5]:
For this tutorial, we will use an auxiliary library of functions (plot_functions.py) that handle the plotting logistics. To get a general impression of how our fault behaved, we plot the time series of the shear stress $\tau$, state $\theta$, and the maximum slip rate $v_{max}$ recorded on the fault.
In [6]:
# Time series of stress, state, and maximum slip rate on the fault
qdyn_plot.timeseries(p.ot[0])
The simulations typically take a few cycles to "warm-up" and to converge to a stable limit cycle. After the warm-up the behaviour of the fault is independent of our choice for initial values. In cases where the fault slip behaviour is chaotic, no stable limit cycle may ever be attained. In our case, a stable limit cycle is attained after about 4 years. To better see what is going on during each cycle, we plot the evolution of the slip rate on the fault from 4 years onwards:
In [7]:
# Spatio-temporal evolution of slip rates
qdyn_plot.slip_profile(p.ox, warm_up=4*t_yr)
In this plot, the warmer colours indicate higher slip rates. The asperity is positioned in the centre of the fault (around $x = 0$). As the fault is progressively loaded, nucleation starts at the centre and the "crack" grows laterally until it enters the matrix, after which it quickly decelerates and the rupture ceases.
This type of plots is perhaps not immediately intuitive, so it helps to look at an animation to get a better sense of the slip evolution on this fault.
In [8]:
# This will take a minute or two...
qdyn_plot.animation_slip(p.ox, warm_up=4*t_yr)
Out[8]:
In [ ]: