Slip on a planar fault in a halfspace

This is the first and simplest example of Tectosaur. Here, we'll solve for the halfspace surface displacement caused by a Gaussian slip field on a planar fault. we'll first solve the problem using Okada dislocations. Then, we will solve it using Tectosaur and compare the two solutions.

To start out, let's import the necessary modules. We use standard scientific Python packages: numpy, scipy, matplotlib. The okada_wrapper module is a simple Python wrapper around the original Okada Fortran code. Finally, import Tectosaur!

Now, let's make a mesh! Here, the tct.make_rect function produces a triangulated rectangle mesh, with n points on both the x and y dimensions and with the corners specified. And, we plot the mesh in the x-y plane just to make sure everything worked right! The plotting code might help understanding the layout of the mesh structure. The surf object is a tuple. The first element is the array of mesh points with shape (n_points, n_dims = 3). The second element is the array of triangles with shape (n_triangles, n_corners = 3).

And we'll do the same thing in the x-z plane for the fault.

And let's specify the physical setup of the problem. First off, we use a shear modulus of 1.0 and Poisson ratio of 0.25. While 1.0 is an absurdly small shear modulus, there is nothing numerically different between using a shear modulus of 1e0 and 1e10 -- it simply multiplies any stresses or tractions.

We specify and plot a Gaussian slip pulse centered at (0.0, -1.0).

Let's solve for the surface displacement using Okada dislocations. I discretize the fault into a 30x30 rectangular grid, and then for each observation point in the surface mesh, I loop over the fault elements and calculate the elastic effect.

Now, that we've set up and solved our problem using Okada dislocations, we get to solve using Tectosaur! We'll build up the integral equation and then solve using GMRES, an iterative linear solver.

I'll do a broad strokes overview of the Symmetric Galerking Boundary Element Method (SGBEM) approach to solving this problem. For a more detailed introduction to the method, look at the book Symmetric Galerkin Boundary Element Method (Sutradhar, Paulino and Gray 2008). To start, let's start from Somigliana's identity for a domain with a crack:

$$u(x) + \int_{S} T^{*}(x,y) u(y) dS + \int_{F} T^{*}(x,y) s(y) dF = \int_{S} U^*(x,y) t(y) dS$$

where $S$ is the free surface, $F$ is the fault, $u$ is the displacement, $s$ is the fault slip, $t$ is the surface traction and $U^*$ and $T^*$ are the respective Green's functions.

We can simplify this equation given that we know $t = 0$ on $S$ and that $s$ on $F$ is given:

$$u(x) + \int_{S} T^{*}(x,y) u(y) dS + \int_{F} T^{*}(x,y) s(y) dF = 0$$

The form of the second two integrals is identical, so instead we can use a single integral:

$$u(x) + \int_{\hat{S}} T^{*}(x,y) \hat{u}(y) d\hat{S} = 0$$

where $\hat{S} = S \cup F$ and

$$\hat{u}(x) = \begin{cases}u(x) & x \in S\\s(x) & x \in F\end{cases}$$

At the moment, we have an integral equation that gives us the displacement anywhere in the volume as a function of fault slip and surface displacement. We'd like to transform that to an integral equation that allows solving for the surface displacement from the fault slip. The first step is the enforce the integral equation on the surface multiplied by a test function:

$$\int_{S}\phi(x)\big[\frac{u(x)}{2} + \int_{\hat{S}} T^{*}(x,y) \hat{u}(y) d\hat{S}\big]dS = 0$$

(For reasons to do with boundary limits, we gain a factor of $\frac{1}{2}$.)

Then, discretizing all the fields with linear basis functions ($\phi_i(x)$) over a triangulated mesh and choosing the test functions to be equal to those basis functions results in the standard SGBEM. We get

$$\begin{align}M_{ij}u_{j} + T_{ik}\hat{u}_k = 0\\ M_{ij} = \int_{S}\phi_i(x)\phi_j(x)dS \\ T_{ik} = \int_S\phi_i(x)\int_{\hat{S}}T^*(x,y)\phi_k(y) d\hat{S} \\ \hat{u}(x) = \sum_{j}\phi_j(x)\hat{u}_j\end{align}$$

To put this into practice, we'll create these two operators using Tectosaur. First we'll create the $T$ operator.

Creating the $M$ "mass" operator is much simpler. We simply provide a Gauss quadrature order and the mesh.

Let's check out the matrix. It's very sparse! That's because $\phi_i(x)\phi_j(x)$ is only nonzero when the two basis functions are defined on the same triangle.

Wait, why didn't we look at the matrix for T? That's because it's partially a matrix-free operator representation. The nearfield is stored in a sparse matrix form. But, the farfield portions of the matrix are never stored. Because we can use a very low order quadrature to generate the farfield elements, the memory bandwidth requirements of storing those matrix elements are greater than the computational cost of recomputing each element every time it's needed. This is particularly true because we make heavy of GPUs where doing enormous amounts of simple arithmetic is very computationally cheap. A nice side effect of this design choice is that Tectosaur doesn't use very much RAM!

Ok, let's put together the mass and T operators into a single operator. Because the T operator is matrix free, we have to use SumOp and MultOp so that the summation and multiplication is done whenever a matrix vector product is needed, rather than right now on the matrix elements.

Next, we need to set up the boundary condition on the fault. We already have the gauss_slip_fnc for fault slip! The question then is how to calculate the coefficients of the basis functions, $\hat{u}_j\phi_j(x)$. For linear basis functions, we can just evaluate gauss_slip_fnc at the corners of each triangle. The values at the corners are the degrees of freedom (DOFs).

We'll calculate the slip for each DOF in the mesh. And we create a list of boundary conditions constraints to impose on our linear system using tct.all_bc_constraints.

Because we are using a linear basis for the displacement on the free surface, there are several DOFs for most points in the mesh -- each triangle that touches a point has a DOF at that point. We'd like to impose continuity of displacement on the free surface, which requires equality between all the DOFs that share a point. But, we also need to ensure that there's a discontinuity across the fault. So, we also pass information about the fault mesh!

So, what do we do with all these constraints? We map from the original unconstrained DOFs to a new set of constrained DOFs. The tct.build_constraint_matrix function uniquely defines a new (smaller) set of DOFs and constructs a matrix to transform from one set of DOFs to the another. Suppose we have 2000 original unconstrained DOFs, and 700 independent constraints. There will be 1300 constrained DOFs. So, the constraint matrix cm will be shaped (2000, 1300). If any of those constraints have a non-zero constant offset (e.g. $x + y = 3$), then the c_rhs vector will contain those values. To transform from constrained DOFs to unconstrained DOFs: cm.dot(x) + c_rhs. So, if we start out with the linear system:

$$Ax = b$$

then in terms of the constrained DOFs, we have

$$A(Cy + r) = b$$

where $C$ is the constraint matrix, cm and $r$ is the vector of offsets, c_rhs. Next, we gain regain symmetry and square matrices by multiplying by $C^T$:

$$C^TA(Cy + r) = C^Tb$$

And rearranging, the final constrained linear system will be:

$$C^TACy = C^Tb - C^TAr$$

These next few lines construct this constrained linear system. First, $C^Tb = 0$ for this problem, so we ignore that term. What's left?

Building $C$:

And building $C^TAr$:

And things are getting exciting! Before we solve the linear system, we provide scipy.sparse.linalg with the info needed for our custom matrix vector product implementation.

And solve!

And the last step is to calculate the unconstrained solution from the constrained solution.

You did it! You solved for the surface displacement with Tectosaur!!!

Let's plot up the solution. I want to make some matplotlib.tricontour plots. To do that, I need values for each point rather than for each triangle vertex. So, we convert from the solution DOFs to point values using the triangle array as the mapping.

Also, to easily plot both the Tectosaur and Okada solutions, I extend the Okada solution with zeros so that it has the same length as the Tectosaur solution (the total number of points in the full fault + surface mesh, not just the number of points in the surface mesh).

A quick function to plot all three components of the displacement.

The match is very nice!

Let's plot the difference.

In general, the error is small. It is largest near the fault where a finer discretization would help resolve the large displacement gradients. You can also see there's some larger errors near the edge of the surface in the $u_z$ component. That's because we're trying to match Okada dislocations, which assume an infinite planar free surface, with a finite planar free surface. We could minimize that component of the error by increasing the size of the surface mesh. However, for practical Earth Science problems where 1\% error is acceptable, using a free surface that is 2-5x the fault length seems sufficient. In addition, the infinite planar free surface is itself an approximation to a spherical Earth. Using Tectosaur, it's better to just embrace that spherical geometry!

Alright, we're all done here. We've replicated the surface displacement results of using Okada dislocation by using a very different type of BEM based on numerical integration and linear basis functions!