```python
from tectosaur2.nb_config import setup

setup()
```

# Quasidynamic earthquake cycles on two faults

## Summary 

* Much like the last section, we simulate quasidynamic rupture evolution but this time on two faults, one of which is curved.
* Most of the code here follows almost exactly from [the previous section on strike-slip/antiplane earthquake cycles](part6_qd).
* Realistic tectonic loading is a challenge when there is more than one fault. 
* The boundary integral approach used here is barely more complex with two faults than with one fault.

## Loading a multi-fault system

In the single fault system from the last section, we were able to gloss over loading by using a slip deficit or backslip formulation. Because the physical system is invariant with respect to a rigid body displacement, any tectonic motion that has already been fully accomodated by fault slip can be ignored. Put differently, in a single fault system with a purely elastic rheology, the stresses from tectonic motion can be equivalently described by backwards slip on the fault.

\begin{equation}
u_{\textrm{total}} = u_{\textrm{loading}} + u_{\textrm{fault}}
\end{equation}

There are a number of different options for loading a fault system:
- **backslip**: as discussed, we simply apply a backwards slip rate to account for distant plate motion. Once the stresses from this backwards slip accumulate, the fault slips forward.
- **farfield displacement**: apply a constant velocity boundary condition at a distance from the fault.
- **stress**: Either on a lower boundary of the domain or on distant side boundaries, apply a constant traction/stress loading rate. 
- **block model**: this is a more advanced version of backslip loading where each fault-bounded subregion of the domain has its own block rotation/motion rate. These rates are assumed to be the long term slip rate. A corresponding kinematic slip deficit is applied to each fault. 

Here, I do something a bit like the block model solution by just applying a backslip of 31 mm of loading per year (1 nm of loading per second) to both faults equally. This would model a real-world situation with a plate boundary that has 62 mm per year of strike-slip motion. Exactly 50% of that motion is accomodated by each of the two faults.


```python
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from tectosaur2 import (
    gauss_rule,
    integrate_term,
    panelize_symbolic_surface,
    concat_meshes,
)
from tectosaur2.laplace2d import double_layer, hypersingular

```


```python
shear_modulus = 3.2e10

quad_rule = gauss_rule(5)
t = sp.var("t")

fault_bottom = 40000
fault1 = panelize_symbolic_surface(
    t, -5000 + t * 0, fault_bottom * (t + 1) * -0.5,
    quad_rule,
    n_panels=200
)

fault2_x = 5000 + (t * 10000 + 10000) + (t ** 2 * 5000 - 5000)
fault2 = panelize_symbolic_surface(
    t, fault2_x, fault_bottom * (t + 1) * -0.5,
    quad_rule,
    n_panels=200 
)
faults = concat_meshes((fault1, fault2))

surf_half_L = 100000
free = panelize_symbolic_surface(
    t, -t * surf_half_L, 0 * t,
    quad_rule,
    n_panels=200
)
```


```python
print(
    f"The free surface mesh has {free.n_panels} panels with a total of {free.n_pts} points."
)
print(
    f"The fault mesh has {faults.n_panels} panels with a total of {faults.n_pts} points."
)
```

    The free surface mesh has 200 panels with a total of 1000 points.
    The fault mesh has 400 panels with a total of 2000 points.



```python
plt.plot(free.pts[:,0]/1000, free.pts[:,1]/1000, 'k-')
plt.plot(fault1.pts[:,0]/1000, fault1.pts[:,1]/1000, 'r-')
plt.plot(fault2.pts[:,0]/1000, fault2.pts[:,1]/1000, 'r-')
plt.xlabel(r'$x ~ \mathrm{(km)}$')
plt.ylabel(r'$y ~ \mathrm{(km)}$')
plt.axis('scaled')
plt.xlim([-100, 100])
plt.ylim([-80, 20])
plt.show()
```


```python
singularities = np.array(
    [
        [-surf_half_L, 0],
        [surf_half_L, 0],
        [5000, 0],
        [-5000, 0]
    ]
)

(free_disp_to_free_disp, fault_slip_to_free_disp), report = integrate_term(
    double_layer, free.pts, free, faults, tol=1e-11, singularities=singularities, return_report=True
)
fault_slip_to_free_disp = fault_slip_to_free_disp[:, 0, :, 0]
free_disp_solve_mat = (
    np.eye(free_disp_to_free_disp.shape[0]) + free_disp_to_free_disp[:, 0, :, 0]
)
```


```python
(free_disp_to_fault_stress, fault_slip_to_fault_stress), report = integrate_term(
    hypersingular,
    faults.pts,
    free,
    faults,
    tol=1e-10,
    singularities=singularities,
    return_report=True,
)
fault_slip_to_fault_stress *= shear_modulus
free_disp_to_fault_stress *= shear_modulus

fault_slip_to_fault_traction = np.sum(
    fault_slip_to_fault_stress[:,:,:,0] * faults.normals[:, :, None], axis=1
)
free_disp_to_fault_traction = np.sum(
    free_disp_to_fault_stress[:,:,:,0] * faults.normals[:, :, None], axis=1
)
```

    /Users/tbent/Dropbox/active/eq/tectosaur2/tectosaur2/integrate.py:208: UserWarning: Some expanded integrals reached maximum expansion order. These integrals may be inaccurate.
      warnings.warn(



```python
A = -fault_slip_to_fault_traction
B = -free_disp_to_fault_traction
C = fault_slip_to_free_disp
Dinv = np.linalg.inv(free_disp_solve_mat)

total_fault_slip_to_fault_traction = A - B.dot(Dinv.dot(C))
```


```python
analytical_sd = -np.arctan(-fault_bottom / free.pts[:,0]) / np.pi
numerical_sd = Dinv.dot(-C).sum(axis=1)
rp = (faults.pts[:,1] + 1) ** 2
analytical_ft = -(shear_modulus / (2 * np.pi)) * ((1.0 / (faults.pts[:,1] + fault_bottom)) - (1.0 / (faults.pts[:,1] - fault_bottom)))
numerical_ft = total_fault_slip_to_fault_traction.sum(axis=1)
plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.title('Stress')
plt.plot(faults.pts[:,1], np.log10(np.abs(numerical_ft)), 'b-')
plt.xlim([-fault_bottom, 0])
plt.xlabel('$\mathrm{depth ~~ (m)}$')
plt.ylabel('$t$')
plt.subplot(1,2,2)
plt.title('Displacement')
plt.plot(free.pts[:,0], numerical_sd, 'b-')
plt.xlim([-100000, 100000])
plt.xlabel('$x ~~ \mathrm{(m)}$')
plt.ylabel('$u$')
plt.show()
```

## Rate and state friction


```python
siay = 31556952 # seconds in a year
density = 2670  # rock density (kg/m^3)
cs = np.sqrt(shear_modulus / density)  # Shear wave speed (m/s)
eta = shear_modulus / (2 * cs)  # The radiation damping coefficient (kg / (m^2 * s))
Vp = 1e-9  # Rate of plate motion
sigma_n = 50e6  # Normal stress (Pa)

# parameters describing "a", the coefficient of the direct velocity strengthening effect
a0 = 0.01
amax = 0.025
fy = faults.pts[:, 1]
H = 15000
h = 3000
fp_a = np.where(
    fy > -H, a0, np.where(fy > -(H + h), a0 + (amax - a0) * (fy + H) / -h, amax)
)

# state-based velocity weakening effect
fp_b = 0.015

# state evolution length scale (m)
fp_Dc = 0.008

# baseline coefficient of friction
fp_f0 = 0.6

# if V = V0 and at steady state, then f = f0, units are m/s
fp_V0 = 1e-6
```


```python
plt.figure(figsize=(3, 5))
plt.plot(fp_a, fy/1000, label='a')
plt.plot(np.full(fy.shape[0], fp_b), fy/1000, label='b')
plt.xlim([0, 0.03])
plt.ylabel('depth')
plt.legend()
plt.show()
```


```python
# The regularized form of the aging law.
def aging_law(fp_a, V, state):
    return (fp_b * fp_V0 / fp_Dc) * (np.exp((fp_f0 - state) / fp_b) - (V / fp_V0))

def qd_equation(fp_a, shear_stress, V, state):
    # The regularized rate and state friction equation
    F = sigma_n * fp_a * np.arcsinh(V / (2 * fp_V0) * np.exp(state / fp_a))

    # The full shear stress balance:
    return shear_stress - eta * V - F

def qd_equation_dV(fp_a, V, state):
    # First calculate the derivative of the friction law with respect to velocity
    # This is helpful for equation solving using Newton's method
    expsa = np.exp(state / fp_a)
    Q = (V * expsa) / (2 * fp_V0)
    dFdV = fp_a * expsa * sigma_n / (2 * fp_V0 * np.sqrt(1 + Q * Q))

    # The derivative of the full shear stress balance.
    return -eta - dFdV


def solve_friction(fp_a, shear, V_old, state, tol=1e-10):
    V = V_old
    max_iter = 150
    for i in range(max_iter):
        # Newton's method step!
        f = qd_equation(fp_a, shear, V, state)
        dfdv = qd_equation_dV(fp_a, V, state)
        step = f / dfdv

        # We know that slip velocity should not be negative so any step that
        # would lead to negative slip velocities is "bad". In those cases, we
        # cut the step size in half iteratively until the step is no longer
        # bad. This is a backtracking line search.
        while True:
            bad = step > V
            if np.any(bad):
                step[bad] *= 0.5
            else:
                break

        # Take the step and check if everything has converged.
        Vn = V - step
        if np.max(np.abs(step) / Vn) < tol:
            break
        V = Vn
        if i == max_iter - 1:
            raise Exception("Failed to converge.")

    # Return the solution and the number of iterations required.
    return Vn, i
```


```python
for fm in [fault1, fault2]:
    mesh_L = np.max(np.abs(np.diff(fm.pts[:, 1])))
    Lb = shear_modulus * fp_Dc / (sigma_n * fp_b)
    hstar = (np.pi * shear_modulus * fp_Dc) / (sigma_n * (fp_b - fp_a))
    print(mesh_L, Lb, np.min(hstar[hstar > 0]))
```

    53.84693101057201 341.3333333333333 3216.990877275949
    53.84693101057201 341.3333333333333 3216.990877275949


## Quasidynamic earthquake cycle derivatives


```python
from scipy.optimize import fsolve
import copy

init_state_scalar = fsolve(lambda S: aging_law(fp_a, Vp, S), 0.7)[0]
init_state = np.full(faults.n_pts, init_state_scalar)
tau_amax = -qd_equation(amax, 0, Vp, init_state_scalar)

init_traction = np.full(faults.n_pts, tau_amax)
init_slip_deficit = np.zeros(faults.n_pts)
init_conditions = np.concatenate((init_slip_deficit, init_state))
```


```python
class SystemState:

    V_old = np.full(faults.n_pts, Vp)
    state = None

    def calc(self, t, y, verbose=False):
        # Separate the slip_deficit and state sub components of the
        # time integration state.
        slip_deficit = y[: init_slip_deficit.shape[0]]
        state = y[init_slip_deficit.shape[0] :]

        # If the state values are bad, then the adaptive integrator probably
        # took a bad step.
        if np.any((state < 0) | (state > 2.0)):
            return False

        # The big three lines solving for quasistatic shear stress, slip rate
        # and state evolution
        tau_qs = init_traction - total_fault_slip_to_fault_traction.dot(slip_deficit)
        V = solve_friction(fp_a, tau_qs, self.V_old, state)[0]
        dstatedt = aging_law(fp_a, V, state)
        if not np.all(np.isfinite(V)):
            return False
        self.V_old = V

        slip_deficit_rate = Vp - V
        out = slip_deficit, state, tau_qs, V, slip_deficit_rate, dstatedt
        self.state = out
        return self.state


def plot_system_state(SS):
    """This is just a helper function that creates some rough plots of the
    current state to help with debugging"""
    slip_deficit, state, tau_qs, slip_deficit_rate, dstatedt = SS

    plt.figure(figsize=(15, 9))
    plt.subplot(2, 3, 1)
    plt.title("slip")
    plt.plot(faults.pts[:, 1], slip_deficit)

    plt.subplot(2, 3, 2)
    plt.title("state")
    plt.plot(faults.pts[:, 1], state)

    plt.subplot(2, 3, 3)
    plt.title("tau qs")
    plt.plot(faults.pts[:, 1], tau_qs)

    plt.subplot(2, 3, 4)
    plt.title("slip rate")
    plt.plot(faults.pts[:, 1], slip_deficit_rate)

    plt.subplot(2, 3, 6)
    plt.title("dstatedt")
    plt.plot(faults.pts[:, 1], dstatedt)
    plt.tight_layout()

    plt.show()


def calc_derivatives(state, t, y):
    """
    This helper function calculates the system state and then extracts the
    relevant derivatives that the integrator needs. It also intentionally
    returns infinite derivatives when the `y` vector provided by the integrator
    is invalid.
    """
    if not np.all(np.isfinite(y)):
        return np.inf * y
    state = state.calc(t, y)
    if not state:
        return np.inf * y
    derivatives = np.concatenate((state[-2], state[-1]))
    return derivatives
```

## Integrating through time


```python
%%time
from scipy.integrate import RK45

# We use a 5th order adaptive Runge Kutta method and pass the derivative function to it
# the relative tolerance will be 1e-11 to make sure that even
state = SystemState()
derivs = lambda t, y: calc_derivatives(state, t, y)
integrator = RK45
atol = Vp * 1e-6
rtol = 1e-11
rk = integrator(derivs, 0, init_conditions, 1e50, atol=atol, rtol=rtol)

# Set the initial time step to one day.
rk.h_abs = 60 * 60 * 24

# Integrate for 1000 years.
max_T = 1000 * siay

n_steps = 1000000
t_history = [0]
y_history = [init_conditions.copy()]
for i in range(n_steps):

    # Take a time step and store the result
    if rk.step() != None:
        raise Exception("TIME STEPPING FAILED")
        break
    t_history.append(rk.t)
    y_history.append(rk.y.copy())

    # Print the time every 5000 steps
    if i % 5000 == 0:
        print(f"step={i}, time={rk.t / siay} yrs, step={(rk.t - t_history[-2]) / siay}")
    if rk.t > max_T:
        break

y_history = np.array(y_history)
t_history = np.array(t_history)
```

    step=0, time=0.002737907006988508 yrs, step=0.002737907006988508
    step=5000, time=161.2388012289919 yrs, step=6.068323795478569e-11
    step=10000, time=161.238801483585 yrs, step=3.200375945922213e-11
    step=15000, time=177.06735503165402 yrs, step=0.05782401363757997
    step=20000, time=187.60702659014385 yrs, step=8.29256996752649e-11
    step=25000, time=187.6070269264095 yrs, step=2.786049654906221e-10
    step=30000, time=213.4153146695212 yrs, step=1.3348499105890855e-10
    step=35000, time=229.49165911338014 yrs, step=1.8458825569115088e-10
    step=40000, time=298.5155537603046 yrs, step=1.1030568652139828e-10
    step=45000, time=298.51555399934733 yrs, step=3.251751197178755e-11
    step=50000, time=298.5168884817506 yrs, step=6.51359998844683e-06
    step=55000, time=323.01280183052285 yrs, step=1.2281707123916784e-10
    step=60000, time=323.01280207887106 yrs, step=3.3484575524851865e-11
    step=65000, time=323.2817940109216 yrs, step=0.001729718400292269
    step=70000, time=351.52626665410605 yrs, step=1.0785478482910092e-08
    step=75000, time=368.9516542294202 yrs, step=1.2976784052681758e-10
    step=80000, time=368.951654876993 yrs, step=6.515590688770813e-11
    step=85000, time=428.9915445603336 yrs, step=0.0014206634043130257
    step=90000, time=446.2939275656298 yrs, step=6.588120455250637e-11
    step=95000, time=446.29392801977554 yrs, step=9.456068304806992e-10
    step=100000, time=469.8994523682489 yrs, step=1.0256917809688377e-10
    step=105000, time=469.89945268723017 yrs, step=5.1798341561007304e-11
    step=110000, time=469.9812354270467 yrs, step=0.00034767164051762437
    step=115000, time=500.90614938849876 yrs, step=1.439715864624497e-10
    step=120000, time=573.9064649517662 yrs, step=1.424001081887202e-10
    step=125000, time=573.9064651906277 yrs, step=2.5022769435539118e-11
    step=130000, time=597.5674917048766 yrs, step=0.0012973735857483035
    step=135000, time=597.6230358626846 yrs, step=5.391379308333549e-11
    step=140000, time=597.6230366893858 yrs, step=4.587749495737104e-09
    step=145000, time=625.2112327038584 yrs, step=1.1906969997104363e-10
    step=150000, time=642.3455004578122 yrs, step=1.2898210138995284e-10
    step=155000, time=704.8417989867936 yrs, step=0.003113930973116129
    step=160000, time=716.3612963754118 yrs, step=2.9978970144993726e-11
    step=165000, time=716.3612966987445 yrs, step=3.9843018386249724e-10
    step=170000, time=743.0069768079197 yrs, step=1.6222491102653862e-10
    step=175000, time=743.0069772376432 yrs, step=4.6660816435353134e-11
    step=180000, time=743.0070823319467 yrs, step=6.836717438689674e-07
    step=185000, time=772.4229461837367 yrs, step=1.1205848921132735e-10
    step=190000, time=790.5550252030326 yrs, step=9.598105764163313e-11
    step=195000, time=790.5550257931098 yrs, step=9.694812119469745e-11
    step=200000, time=790.5770082962302 yrs, step=0.00011880289499414454
    step=205000, time=843.4093589520743 yrs, step=1.3986156636192637e-10
    step=210000, time=843.4093595399053 yrs, step=5.7540281407326665e-11
    step=215000, time=865.8601485323139 yrs, step=4.5113514750450234e-10
    step=220000, time=903.1263077569464 yrs, step=1.2438854951289735e-10
    step=225000, time=903.1263082887026 yrs, step=8.413452911659529e-11
    step=230000, time=927.4523165504868 yrs, step=0.008045238903279458
    step=235000, time=947.0642174585033 yrs, step=1.7528026899290684e-11
    step=240000, time=947.0685830652857 yrs, step=1.768826485736451e-05
    step=245000, time=972.5636479272499 yrs, step=4.078590535048743e-10
    step=250000, time=988.09794803525 yrs, step=8.195863612220059e-11
    step=255000, time=988.0979484325295 yrs, step=2.7077779485800782e-11
    CPU times: user 2h 15min 3s, sys: 36min 52s, total: 2h 51min 55s
    Wall time: 1h 48min 42s


## Plotting the results


```python
derivs_history = np.diff(y_history, axis=0) / np.diff(t_history)[:, None]
max_vel = np.max(np.abs(derivs_history), axis=1)
plt.plot(t_history[1:] / siay, np.log10(max_vel))
plt.xlabel('$t ~~ \mathrm{(yrs)}$')
plt.ylabel('$\log_{10}(V)$')
plt.show()
```


```python
for j in range(2):
    fm = fault1 if j == 0 else fault2
    plt.figure(figsize=(10, 4))
    last_plt_t = -1000
    last_plt_slip = init_slip_deficit
    event_times = []
    for i in range(len(y_history) - 1):
        y = y_history[i]
        t = t_history[i]
        slip_deficit_full = y[: init_slip_deficit.shape[0]]
        fault1_sd = slip_deficit_full[:fault1.n_pts]
        fault2_sd = slip_deficit_full[fault1.n_pts:]
        slip_deficit = fault1_sd if j == 0 else fault2_sd
        should_plot = False

        # Plot a red line every three second if the slip rate is over 0.1 mm/s.
        if (
            max_vel[i] >= 0.0001 and t - last_plt_t > 3
        ):
            if len(event_times) == 0 or t - event_times[-1] > siay:
                event_times.append(t)
            should_plot = True
            color = "r"

        # Plot a blue line every fifteen years during the interseismic period
        if t - last_plt_t > 15 * siay:
            should_plot = True
            color = "b"

        if should_plot:
            # Convert from slip deficit to slip:
            slip = -slip_deficit + Vp * t
            plt.plot(slip, fm.pts[:,1] / 1000.0, color + "-", linewidth=0.5)
            last_plt_t = t
            last_plt_slip = slip
    plt.xlim([0, np.max(last_plt_slip)])
    plt.ylim([-40, 0])
    plt.ylabel(r"$\textrm{z (km)}$")
    plt.xlabel(r"$\textrm{slip (m)}$")
    plt.tight_layout()
    plt.savefig("halfspace.png", dpi=300)
    plt.show()
```
