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, 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. 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. That being said, every effort in the direction of accelerating the computations on the ﬂow side are a necessary step forward in rendering a fast coupled geodynamics simulator. In this document, we present a framework in which we study the porting of the OpenMP parallelism of the ﬂow simulator to a GPGPU 1 Introduction The earthquake cycle, from slow deformation associated with interseismic behavior to rapid de- formation 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]. Porting OpenMP to GPGPU to accelerate Krylov space computations in coupled geodynamics simulator 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, 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 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 [13] 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. Algorithm 1: Time-marching /* Initialization */ 1whilet<T do /* time marching till final time */ 2 while Not converged do /* Staggered solution algorithm loop */ 3 while Not converged do /* Flow loop */ 4 Solve ﬂow problem for pressures using GMRES; // Krylov subspace solver in OpenMP parallel flow simulator 5 Port information from ﬂow to geodynamics; 6 Solve geodynamics problem for displacements using GMRES; // Krylov subspace solver in MPI parallel geodynamics simulator 7 Port information from geodynamics to ﬂow; 8t t+ t; The simulator is a mixed programming language hybrid MPI and OpenMP framework in which the open source geodynamics code [14, 15] written in Python and C++ is MPI-based and the in-house multiphase/multicomponent ﬂow code written in C++ is OpenMP-based [16]. As shown in Algorithm 1, the framework follows the popular staggered solution algorithm which allows for different simulators written in separate languages with respective parallelism frameworks are coupled together to render a coupled simulator [16 –26]. 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 [27] 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 [28] XThe SWIG [29] code with all the .i ﬁles, and SWIG-generated python code and MODU- LENAME_wrap.cxx which wraps the Python code around the C++ code to create python modules. 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 2 Porting OpenMP to GPGPU to accelerate Krylov space computations in coupled geodynamics simulator ﬂ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 [30] on the ﬂow side, and the second option is to ofﬂoad the OpenMP parallel computations onto a GPGPU [31]. 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 [32] 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