Understanding the causality between the events leading upto and post fault slip and the earthquake recording is important for seismic design and monitoring of underground struc- tures, bridges and reinforced concrete buildings as well as climate mitigation projects like carbon sequestration and energy technologies like enhanced geothermal systems or oilﬁeld wastewater disposal. While the events leading upto fault slip are typically governed by poroelastostatics, the events post fault slip can easily transition into poroelastodynamics territory due to runaway fault slip velocities. There are marked differences in the numerics of poroelastostatics and poroelastodynamics, and a simple switch from one algorithm to another based on fault slip velocities is not trivial. In fact, an understanding of expected fault slip velocities is critical apriori, as an algorithm which can seamlessly transition from time marching in poroelastostatics realm to poroelastodynamics realm and vice-versa is extremely difﬁcult to achieve. We present the numerics of both physics and point out the differences between the two in this work. 1 Introduction The earthquake cycle, from slow deformation associated with interseismic behavior to rapid deformation associated with earthquake rupture, spans spatial scales ranging from fractions of a meter associated with the size of contact asperities on faults and individual grains to hundreds of kilometers associated with plate boundaries [1]. Similarly, temporal scales range from fractions of a second associated with slip at a point during earthquake rupture to hundreds of years of strain accumulation between earthquakes. In many cases, earthquakes are triggered after pore pressure perturbations activate critically stressed seismogenic faults [1], not just due to natural causes like earth tides [2], rainfall [3], snowfall [4], typhoons [5], but also due to human activity [6]. As faults slip, the induced stress ﬁeld spawns seismic waves, which travel through and around the earth and can be recorded with seismometers. There are several different kinds of seismic waves, and they all move in different ways. The two main types of waves are body waves (P- and S-) and surface waves. Body waves can travel through the Earth’s inner layers, but surface waves can only move along the surface of the planet like ripples on water. Understanding the causality between the events leading upto and post fault slip and the earthquake recording is important for seismic design and monitoring of underground structures [7], bridges [8] and reinforced concrete buildings [9] as well as climate mitigation projects like carbon sequestration [10] and energy technologies like enhanced geothermal systems [11] or oilﬁeld wastewater disposal [12]. The complexity of the many physical processes operating over this vast range of scales leads most researchers to focus on a narrow space-time window to isolate just one or a few processes; the limited spatial and temporal coverage of observations also often justiﬁes this narrow focus. For quasi-static simulations [13 –24], the inertial term is ignored and time dependence only enters through the constitutive models and the loading conditions. As a result, a quasi-static simulation is the solution to a series of static problems with potentially time-varying physical properties and boundary conditions. In dynamic simulations, we shall include the Comparing poroelastostatics and poroelastodynamics: Numerics, solvers and algorithms inertial term to resolve the propagation of seismic waves, with an intended focus on applications for earthquake physics and ground-motion simulations. The only place where things can accelerate fast is on the fault and the only time it happens is when fault slips seismically ( >1 cm/sec). If the fault slip is aseismic or creeping (<0.1 mm/sec), then accelerations are negligible and quasi-static and quasi-dynamic [25] approximations are valid. In our models, we know the fault is the location where things will happen, but we do not know a priori is whether it will be seismic, aseismic or in the transition region. While back-of-the-envelope estimations can be done from the strain energy accumulated in the rock surrounding the faults, the fault’s friction properties, and the rock’s elastic moduli, the more optimal and accurate option is to use the simulator to do this estimation by equipping it with the physics of elastodynamics. 1.1 Code development The simulator is a mixed programming language hybrid MPI and OpenMP framework in which the open source geodynamics code [26, 27] written in Python and C++ is MPI-based and the in-house multiphase/mul- ticomponent ﬂow code written in C++ is OpenMP-based [13, 14]. There are three layers in which the code framework operates out of separate directories XThe high-level Python code in the geodynamics module makes use of Pyre [28] to link the sub- modules together at runtime XThe low-level C++ code performs operations on matrices and vectors in parallel using the Standard Template Library [29] XThe SWIG [30] code with all the .i ﬁles, and SWIG-generated python code and MODULE- NAME_wrap.cxx which wraps the Python code around the C++ code to create python modules. The numpy library [31] allows direct manipulation in Python of multi-dimensional arrays, using memory allocated in C++. In general, Python’s ability to easily expose a library API and the use of numpy arrays for data interchange, make it an excellent tool for combining libraries at a higher level to attack more complex multi-model, multi-physics, multi-scale problems [32]. Some of the code snippets are provided in A. Changes made to the C++ code entail a makeﬁle related compilation of the code base. If there are changes to the C++ header ﬁles that need to be exposed to the input ﬁle, then those changes have to be copied to the corresponding module’s .i ﬁle. This means SWIG has to be run again on that module to update MODULENAME_wrap.cxx and Python ﬁles of the module before building the executable using its Makeﬁle. The central disadvantage for such multi-language code frameworks comes in debugging. There are certainly hurdles introduced into the build system, since different compilation and links steps are needed and many more tests are needed to verify interoperability. Currently, the code base is hybrid MPI and OpenMP, in which the geodynamics code is MPI-parallel and the ﬂow code is OpenMP-parallel. This hybrid framework, although accurate, brings in efﬁciency issues for large size simulations. While the geodynamics elements are distributed among MPI processes with each process running the same tasks on a different sub-grids of the main grid, the ﬂow elements are all stacked together in the root process, and the tasks are distributed among threads in the root process. When the data structures port information from the geodynamics code to the ﬂow code, they are assembled into the root process using the MPI_Gather command, and then transfered to the ﬂow simulator. This porting can happen multiple times in each time step depending on the convergence requirements of the staggered solution algorithm, and this adds inefﬁciency to the overall computational run-time. The ﬁrst option to increase efﬁciency is to port the OpenMP parallel functionality to MPI parallel on the ﬂow side, and the second option is to ofﬂoad the OpenMP parallel computations onto a GPGPU by incorporating CUDA and OpenCL libraries [33]. The second option is deﬁnitely more engaging as it leverages the massive computational power of the GPGPU. That being said, the advent of new, massively parallel architectures [34] has greatly increased the penalty for suboptimal non-vectorized code while providing much more energy efﬁcient performance. They also increase the imbalance between ﬂop rate and memory bandwidth or latency, so that thousands of ﬂops can be needed to cover outstanding memory