Introduction

SpecFem3D Cartesian is an open source software for seismic wave propagation in 3D. It is based upon the Spectral Elements Method (SEM, we will discuss what it is all about later in this article), encompasses its own mesher and can also perform Full Waveform Inversion among other things. It is flexible and powerful enough to run seismic simulations at a local and regional scale. The software was developed by researchers from the CNRS, Princeton and ETH Zürich since the late 90’s. Note that if you feel like a regional scale is not enough for you, SpecFem3D Globe also exists.

Use cases

High performance seismic simulations have a lot of real-world applications and have been on the forefront of the race for High Performance Computation (HPC) for a long time now. Those simulations can allow for a better understanding of seismic risks and help design human infrastructures resilient to the worst cases. Particularly in regions of the world with high seismic activity like California, Mexico or even Chile and even more for infrastructures whose structural integrity is critical like nuclear power plants (i.e. Fukushima or the ITER project near Cadarache in France).

7.0 magnitude earthquake simulation on the Hayward fault near San Fransisco in California. source.

But the biggest commercial application may be the use of numerical simulations for hydrocarbon exploration. One classical strategy to find whether we have oil or gas under our feet without drilling expensive holes randomly into the ground, is to send “seismic” waves into the ground (using dynamite, air guns, vibrators…) and to listen for the reflection of those waves when they come back to the surface similarly to how sonars or medical imaging works. This data is then fed to the software that will compute a subsurface model that can orient drilling campaigns. Seismic imaging can also be used to gain knowledge on fault areas to further the knowledge of seismic risks and improve security measures.

Classical hydrocarbon exploration data collection

Each geophone (or seismometer in the case of earthquakes) will record ground movement. This data is usually displayed with each recording side by side as pictured below :

Example of Seismic/Exploration records.

This is is the data that is then processed by software to deduce ground composition.

Velocity models of marine soils using Full Waveform Inversion (FWI) and First-Arrival Tomography (FAT) from data collected by ocean bottom seismometers of a crossection of the Nankai Trough : source

Computational Physics and Numerical methods

Classical FEM method

In the world of computational physics and numerical simulations, the Finite Elements Method (FEM) is one of the standard methods to approximate the solutions to Partial Differential Equations (PDEs) which are present in a lot of problems such as the spread of heat through a material, fluid flow, electromagnetic waves and in our case, seismic waves. PDEs are equations that “imposes the relations between the various partial derivatives of a multivariable function” . In the case of our seismic simulation for example, we are going to look for a function representing the displacement of the soil for each axis in 3D and through time, giving us a 4 variables function. Those PDEs are usually impossible or too cumbersome to solve mathematically, therefore we are going to build approximations through numerical means.

To solve equations numerically, we have to transcribe our continuous problem/equation into a finite model/problem that can be solved by our computer with finite memory space and computing power. Time is discretized by cutting the total time window in regular time steps.

The strategy of the FEM is to divide our large problem through discretization onto smaller parts called finite elements on which we are going to approximate our unknown solution function. Those finite elements build a mesh representing the total domain of the problem. The definition of a mesh is a critical part of the process as it needs to be well defined to avoid numerical errors, optimize computation power etc…

Once our domain has been correctly defined and our mesh built The FEM needs us to rewrite our equation in another form. The workflow looks something like this :

We take a classical Dirichlet problem as an example where we want to find a function that satisfies :

Building the weak formulation by multiplying both sides of the equation by a test function phi and using Green’s Formula :

And defining :

Gives us the variational form :

Under the right circumstances, the Lax-Milgram theorem proves that such function exists. To approximate our solution we will be searching for solutions on a subspace of finite dimension :

And particularly in the FEM, this subspace that we design is going to be the space of piecewise continuous polynomial functions with each of those pieces being one finite element. Now that our solution space is finite and our equation is linear, testing u for every v in that space just becomes testing u for each element of the basis of the space. That gives us N equations and N coefficients to find u as a linear combination of the basis functions. This set of equations can now be numerically solved.

Here on a 1D problem with the “real” function in blue, the FEM approximation in red and the basis functions at the bottom i.e. the red dotted line is a linear combination of the basis functions.

Of course we are approximating the solution and we cannot be sure of how close we are going to be eventually from the real solution. But Céa’s lemma gives us an upper limit of the error that we make.

With our mesh and our variational form, we are ready to give our problem to a FEM software that will be approximating a solution of our problem.

SEM method

Now that we have seen what the FEM is, we can talk a little bit about the Spectral Elements Method (SEM) which is used in SpecFem3D. The method was originally developed in the Computational Fluid Dynamics domain before being eventually implemented on seismic simulations.

Whereas FEM usually uses low-degree polynomials as basis functions to represent functions on an element, the SEM used in SpecFem3d uses Lagrange polynomials of degree 4 to 10 and inside each elements, the function is interpolated at specific points called Gauss-Lobatto-Legendre (GLL) where the Lagrange polynomials reach either 0 or 1 (see below). Using those Lagrange polynomials as base functions has the neat property of being orthogonal with one another and greatly simplifying the numerical algorithms needed to solve our problem numerically. This “version” of the FEM gives some computational easiness and works really well on parallel computers.

Left : representation of Lagrange interpolants of degree N=4. Right : representaion of the N+1 Gauss-Lobatto-Legendre (GLL) interpolation points in an 2D element. Note that on each of these points, the Lagrange interpolants are valued at either 0 or 1.

If you wish to learn more about the SEM and how it is implemented in SpecFem3D and you are not afraid of some equations, you can refer to Tromp, J., Komatitsch, D., & Liu, Q. (2008). Spectral-element and adjoint methods in seismology. Communications in Computational Physics, 3(1), 1-32. or these lecture notes.

Versions

Versions available on Qarnot's cloud platform.

If you are interested in another version, please send us an email at qlab@qarnot.com.

Prerequisites

Before launching the case, please ensure that the following prerequisites have been met.

• Create an account (here)
• Retrieve the authentication token (here)
• Install one of Qarnot’s SDK (Python) (Node.js) (C#) (Commande Line)

Test cases

We will be running one of the examples that can be found from the official SpecFem3d git repository. This simulation demonstrates how SpecFem3D works with coupled acoustic/elastic domains by simulating a seismic source inside a coffee cup. The mesh is already provided in this case.

The example folder contains :

• A DATA folder containing :
  • Par_file : input file configuration
  • STATIONS : coordinates of the receivers that will output time recordings
  • CMTSOLUTION : source definition
• MESH-default : mesh files
• run.sh : batch script that will be running SpecFem3D on Qarnot.

You can download scripts to run these examples from the archive just below. You need to unzip this folder to be able to use it.

Launching a case

Copy the following code in a Python script and save it next to the input folder you unzipped before. Be sure you have copied your authentication token in the script (instead of <<<MY_SECRET_TOKEN>>>) to be able to launch the task on Qarnot.

The script will run the simulation on the Qarnot platform and output the result in the ‘sf3d_results’ on your computer once it is completed.

Results

At any given time, you can monitor the status of your task on Tasq.

Once the task is deployed, it should take around 5 minutes to run, and a few more to download the results to your computer. You can visualize the AVS_movie_*.inp files using Paraview or any other viewer that supports the format. You can also use our Paraview Web payload to view the results online if needed. It must look something like this :

Advanced usage

If you wish to run Specfem3D using GPUs or a cluster from Qarnot, please contact us at qlab@qarnot.com.

Wrapping up

That’s it! If you have any questions, please contact qlab@qarnot.com and we will help you with pleasure !

Did you like this article ? If so you might be interested in those other articles :

• Freefem on Qarnot
• Paraview web on Qarnot
• Gmsh on Qarnot Cloud