Performance Analysis of Speculative Parallel Adaptive Local Timestepping for Conservation Laws

1.1 Previous Work

The idea of allowing different subsystems to advance with different timesteps has been extensively examined for Runge-Kutta methods [38, 60]. Extending these methods to conservation laws has seen considerable work [16, 23, 26, 42, 48, 58, 61, 66]. In addition, other fields have arrived at related time integration schemes. Quantized state systems [19, 31, 47] have been used to model systems of ordinary differential equations (particularly for power electronics modeling) by approximating state values with time-dependent polynomials and updating estimates when they deviate more than a given threshold. Those systems focus on enforcing bounds on state variable deviations, as opposed to enforcing particular mathematical constraints, such as a CFL condition. Asynchronous variational integration is another locally timestepped method designed for linear elastodynamics [50]. Whereas local timestepping for conservation laws is focused on maintaining non-linear stability, asynchronous variational integration focuses on preserving the total energy of the system.

Beyond the various applications, parallelizations share common motifs based on the characterization of the algorithm. By assuming constant timesteps or variations in timestepping sizes, local timestepping methods—also known as multi-rate integrators—can be parallelized using task-graph-like approaches (e.g., [13, 43, 67]). However, the inability to account for variations in the wave speed limits the practicality of these approaches for non-linear problems. Methods such as quantized state systems and the solver in the work of Omelchenko and Karimabadi [55, 56, 57] rely on a discrete event simulation framework to account for changes in timestep and correctly resolve dependencies in a parallel context. Our ALTS implementation differs from these other discrete event simulation based approaches in that (1) ALTS strictly satisfies the total variation diminishing CFL condition leading to a different control flow and (2) previous approaches have not utilized and examined the performance of ALTS using optimistic parallel discrete event simulation (PDES).

To arrive at a performant implementation of ALTS, we employ two types of optimizations: throttling of optimistic execution and dynamic load balancing. The limiting of the progress logical processes are able to make to improve the performance of optimistic PDES simulators is an extensively studied subject. Many previous approaches throttle optimistic PDES by inspecting either simulator state (e.g., [30, 65]) or hardware resource utilization [25]. We refer the reader to the review of Das [24] and therein cited works for a comprehensive overview. A common limitation of these methods is that they fail to fully exploit application specific knowledge but rather can only infer it through statistical methods. In the work of Bauer et al. [3] and Lindén et al. ??, application-specific knowledge is used to provide upper bounds by which neighboring communication must be received. These dynamic local time windows block logical process execution once further progress is guaranteed to result in rollback. This technique is used to improve the performance of PDES-based simulation of the reaction-diffusion master equation. Simulation of conservation laws shares two important characteristics with reaction-diffusion master equation simulation: (1) it is not possible to place a lower bound on event arrival times (precluding the use of lookahead), and (2) it is possible to easily predict scheduling times of future events. Our optimizations rely on a similar mechanism to identify and eliminate misspeculation for the local timestepping algorithm.

Dynamic load balancing for PDES has also been examined extensively in previous work. Earlier work focused mostly on measured compute load of the cores to determine when and how to rebalance. For example, in the work of Boukerche and Das [7], a process migration strategy was presented that balanced workload based on CPU-queue length. Glazer and Tropper [35] present a method to dynamically load balance by migrating processes based on measurements of local clock progress and a bin packing algorithm. The approach presented by Avril and Tropper [2] uses the number of events processed but added estimates of the change in communication as an imbalance metric within the context of digital circuit simulation. In the work of Lindén et al. [52], a voxel-based, load balancing protocol was introduced using an imbalance metric based on rollbacks. More recently, Mikida et al. ?? examined the use of load balancing capabilities in the Charm++ [45] runtime to provide dynamic load balancing capabilities to PDES simulation models in ROSS [17]. They examined the efficacy of both a centralized greedy heuristic as well as a distributed strategy that utilizes a gossip protocol to avoid the synchronization cost of more centralized methods [54]. Our technique has some similarities to these previous methods but offers a more tailored approach given the specific characteristics of our application (e.g., handling wet and dry cells separately) and takes advantage of the 1D spatial domain topology to compute partitions extremely efficiently.

A key outstanding issue has been the adoption of adaptive locally timestepped methods by the scientific computing community. PDES as well as the mathematical foundations of ALTS have been around for more than 30 years, yet neither have managed to gain substantial traction in the scientific computing community [33]. Key barriers to entry remain the productive and maintainable introduction of ALTS into an existing code base while still meeting the performance needs of the higher-performance computing community [40]. Adaptive mesh refinement packages such as Clawpack [21] and ExaHyPE [59] solve conservation laws by requiring the user only to specify key characteristics of the conservation law. Other packages try to provide isolated functionality that can be integrated into a larger software stack—for example, SUNDIALS [41] provides optimized routines for timestep integration of initial value problems. Our proposed API most closely follows the second approach. The objective being to rapidly introduce ALTS into existing applications to demonstrate advantageous speedup and to contain the complexity of ALTS within a timestepping layer.

2.1 Numerical Solutions to Conservation Laws

Conservation laws model the evolution of conserved quantities between regions of the domain. As a partial differential equation, these laws are represented as (1) \( \begin{equation} \partial _t u + \partial _x f(u) = S, \end{equation} \) where u is our conserved quantity (e.g., mass or momentum) and f is the flux, which encodes how those quantities are exchanged with the surrounding region, and S is the source term (e.g., friction). These equations are notoriously difficult to stably simulate, as their solutions may form discontinuities. When discretizing these methods, particular care must be taken to faithfully approximate the exchange of conserved quantities between regions of the domain.

In this work, we only consider first-order methods for 1D domains. Consider a grid, \( x_{lb} = x_0 \lt x_1 \lt \dots \lt x_{n}=x_{ub} \). The first-order finite volume method will approximate the evolution of the average conserved quantity in cell j as (2) \( \begin{equation} U_j(t) - U_j(s) = \Delta x_j \int _s^t S_j(\tau) - \frac{\hat{F}_{j+1/2}(\tau) - \hat{F}_{j-1/2}(\tau)}{\Delta x_j} \,\mathrm{d}\tau . \end{equation} \) Our discretization will assume that conserved quantities are piecewise constant in time. The use of a continuous integral is meant to save on indexing headaches but should be practically understood as a discrete sum that can be evaluated exactly. The numerical flux \( \hat{F}_{j+1/2} \) approximates the exchange of conserved quantities between cells j and \( j+1 \). The maximum speed at which information can pass from one cell into another is denoted as \( |\Lambda | \), which is computed to be the magnitude of the largest eigenvalue of the Jacobian of the flux, \( \partial _u f(u) \).

The purpose of local timestepping is to not require each cell to update at each discrete time. To help identify when cells update, let \( t^m_j \) denote the simulation time of the m-th update on cell j. Note that the state \( U_j \) must satisfy Equation (2) at these update times.

To ensure that the solution is free of oscillations, we require two conditions to be satisfied. These restrictions arise due to mathematical constraints, not limitations of the implementation. First, updates must obey the CFL condition, which limits the allowable time between successive updates [39]. Consider an element j and its neighbors N, and define \( t^* \) as the most recent time that every element has updated—that is, (3) \( \begin{equation} t^* = \min _{n \in N} \max _m \lbrace t_n^m \rbrace . \end{equation} \) A proposed update advancing element j from \( t_j^m \) to \( t_j^{m+1} \) satisfies the CFL condition if (4) \( \begin{equation} t^{m+1}_j - t^* \le C_{CFL} \frac{\Delta x_j}{ \max _{n\in N} |\Lambda _n|(\tau)} \quad \text{ for all } t^* \le \tau \lt t^{m+1}_j, \end{equation} \) where \( C_{CFL} \) is a discretization specific constant, \( |\Lambda | \) is the largest wave speed over the element’s stencil, and \( \Delta x_j \) is the size of the element. This inequality limits the allowable change of the conserved quantities within an element between successive updates by restricting the timestep size relative to flow rates. Second, elements must update in a locally ordered manner. In the ALTS setting, not every element updates at every timestamp. Given two adjacent elements, only one element may update more frequently than the other. To flip which element substeps relative to the other element, the two elements must pass through a synchronization time where both elements update at the same timestamp. We refer to element update traces that obey this behavior as locally ordered.

Listing 1.

View Figure

For performance purposes, we partition elements into submeshes. Submeshes update constituent elements using a single timestep, which is limited by the most stringent local CFL condition. Since elements within the submesh step synchronously, local ordering is trivially satisfied within the submesh. By changing the number of elements per submesh, we gain the ability to tune the computational cost of submesh updates to amortize simulator overheads. As we show later, balancing simulator overheads with computational work is essential to obtaining good performance.

As a discrete event simulation, each actor corresponds to a submesh along with information about neighboring submeshes required to update element densities and compute the CFL condition at the boundary of the submesh. Also referred to as logical processes, actors contain the state of the system, interacting with one another only through messages managed by the PDES engine. Events are scheduled using a schedule routine that accepts a timestamp, actor, and functor as arguments or can be executed within the same event using a schedule_inline routine. The algorithm consists of two types of events: update shown in Listing 1 and push_flux shown in Listing 2. A submesh executes an update event to advance the submesh state to the event execution time (Listing 1, line 10). The updated halo state needs to be shared with neighboring submeshes, which is done via push_flux events (Listing 1, line 16). Note that updates can be scheduled as part of the push_flux event (Listing 2, line 9), leading to a potential cascade of updates. Showing that this interleaved execution of update and push_flux events utilizes the correct state and terminates has been carefully examined in the work of Bremer et al. [10]. The eps delay in scheduling the event push_flux is a small amount of subtime used to correctly order updates and push fluxes. The advance, accumulate, and make_meta_data functions encapsulate the mathematical discretization of the conservation law. The halo regions computed as part of make_meta_data for discontinuous Galerkin-based approaches consist of the state variables evaluated at shared interfaces. Higher-order finite volume methods would require sending larger halo regions. As presented here, ALTS does assume that the halo region is entirely contained in a single neighboring submesh. The next timestep is computed in compute_t_next to be the largest timestep that satisfies (4). Since the timestep depends on the boundary metadata, the next update time may change as push_flux events are processed and the update event must be rescheduled. We know that the most recent actor state will have access to all relevant metadata, and thus we only execute the update consistent with the current submesh.t_next. All other updates become no-ops, either because their scheduled execution time does not match submesh.t_next or the update has already been executed.

The two functions n_outstanding_msgs and inhibits_progress examine the local representation of the distributed state of the system to determine what events should be scheduled next. The n_outstanding_msgs function provides a lower bound on the number of outstanding messages to be processed on the actor by time t. By knowing that a message will be processed before the next update time, we defer scheduling updates until the waited upon messages are received, as any event executed after n_outstanding_msgs that returns a number greater than zero is guaranteed to be rolled back. The inhibits_progress function checks to see whether the actor is unable to make progress without forcing its neighbor to update. Forcing a neighbor to update may arise due to a local ordering violation, an inability to propose a CFL-stable next timestamp, or performance benefits of having the neighbor prematurely update.

Listing 2.

View Figure

2.2 The Devastator Runtime

To measure the performance of the ALTS algorithm, we implemented it using the Devastator PDES runtime (first described and used in the work of Chan et al. [20]). Devastator is a general-purpose, optimistic PDES engine that coordinates event execution and rollback required by Jefferson’s Timewarp protocol. Devastator is built using the GASNet-EX runtime [6] and utilizes light-weight active messages (RPCs) to invoke events on remote logical processes. Since GASNet-EX active messages are tailored to the specific interconnect hardware, they are efficient at delivering the small asynchronous events required of PDES. It is implemented with modern C++ to improve programmer productivity, uses parallel threads and processes to enable distributed memory execution, and includes several optimizations to improve communication and synchronization efficiency. The modern C++ interface allows clean specification of event classes, their serialization, and their execute and rollback functions. Since conservation laws are not reversible in general, portions of the actor state are logged during execute, and the actor state can be restored to the state prior to the execute call by calling the rollback function.

To illustrate the parallel efficiency of the runtime, Figure 1 shows weak scaling performance of Devastator running a modified version of the well-known PHOLD benchmark [32]. PHOLD consists of a number of rays that bounce between actors. The ray arrival times are distributed exponentially, and a ray terminates after a fixed number of events. To determine a ray’s next destination, the actors are lexicographically sorted based on rank ID and then local actor ID, and the target rank is selected by using a wrapped normal distribution with standard deviation of 16 and a mean of the sending actor ID. This modification allows us to tune the communication pattern of the benchmark to more closely resemble applications with locality, such as the ALTS simulations studied here. In the figure, we show the per core event commit rate as the number of cores is varied (with the number of logical processes per core fixed at 1,000). The different series in the figure vary the number of active events (rays) per logical process in the system (corresponding to the average event queue depth for each logical process). Each configuration is run three times, with each trial represented as a point and the line representing the average commit rates of the trials. As can be seen in Figure 1, parallel efficiency (performance per core) decreases to approximately 55% to 75% as the number of cores increases from 32 cores to 16,384 cores, corresponding to a 270× to 3800× overall speedup in aggregate event commit rate. In certain cases, significant variance is observed. Determining whether the jagged nature of the scaling results is due to specific paths taken by the rays or due to system noise will be addressed in a more thorough future treatment of Devastator.

Fig. 1.

View Figure

2.3 Performance Optimizations

Although satisfying the CFL condition and local ordering leads to stable computation, careful management of PDES overheads is required to achieve expected speedups. In this section, we introduce a set of performance optimizations focused on eliminating bad speculation (Sections 2.3.1—2.3.4) and a dynamic load balancing scheme (Section 2.3.5).

2.3.1 Binning Timestep Sizes.

The need for local ordering of events leads to potentially increased forced synchronizations. For example, given submeshes stepping with timesteps 17 and 18, the submesh stepping with size 17 would step to time 17 and then have to update again at time 18 to maintain local ordering. To address this, we bin future timesteps into powers of 2 times a fixed user-supplied \( \Delta t_{\min } \). In this example, both submeshes would take binned timesteps of size 16. To simplify computation all timestamps are kept as integers with units of \( \Delta t_{\min } \). The Pythonic pseudocode for timestep binning is

View Figure

In practice, we implement this efficiently using bit shifts. The impact of this optimization is that we reduce forced updates by having submeshes with roughly the same timestep step to the same timestamps. This also concentrates updates to certain multiples of \( \Delta t_{\min } \) leading to improved simulator parallelism and reduces the depth of the algorithm. The downside of this approach is that we are unable to take the largest feasible timestep, potentially increasing the total number of updates.

3.1 Computation Required by ALTS

The purpose of the Problem class is to specify the physics required to update the conserved variables and compute local wave speeds required to determine timestep restrictions. We present the Problem class interface in Listing 3.

We begin by focusing on the physics information required during a push_flux event. When a submesh updates, the numerical fluxes on neighboring submeshes must be updated to ensure that Equation (2) continues to hold. To achieve this, the updated submesh calls the make_message function to construct a message of type msg_t, which encapsulates the information needed by its neighbors to update their fluxes accordingly. Note that often there is an opportunity to reduce the required message sizes by computing derived quantities before sending them to neighboring submeshes. For example, in the case of the discontinuous Galerkin method, where conserved quantities are approximated by polynomials, the numerical flux only depends on the value of the conserved quantity at the shared interface. In this case, the make_message function computes the necessary derived quantities and sends only those to the submesh’s neighbors as part of the push_flux event.

When the neighboring submesh executes the push_flux event, it accumulates the flux exchanged based on the previous state up until that point in time and computes a new flux as part of the process_message routine. The representation of this partially computed right-hand side is encapsulated with the bdry_t type. Since the integration of this state is done in the local timestepping application layer, the bdry_t type must support standard vector space operations (multiplication by a scalar and addition of two bdry_t objects). Finally, when the neighboring submesh is updated, the boundary flux values are applied to all dependent degrees of freedom as part of the lift_bdry function. To make the preceding processes concrete, for the methods considered in this work, the make_message routine would simply identify the average conserved quantity associated with the cell neighboring the submesh and send that state as a msg_t. When the neighbor executes the associated push_flux event, the previous flux at the shared interface is updated by scaling it based on the elapsed time, and a new current flux is computed in process_message.

Next, we turn to the physics required by the update event. The build_u_next function takes the state’s current time and the desired timestep as an argument and then returns the internal state at time t + dt. In this case, the implementation is left to the domain scientists. We implement this as a simple forward Euler step. In addition, build_u_next is responsible for integrating boundary conditions. The build_u_next function in conjunction with the lift_bdry function discussed in the previous paragraph computes all of the necessary terms required to satisfy Equation (2) for a particular submesh.

One critique of the build_u_next abstraction is that it requires the application to understand how to perform a simple sequential integration, undermining the physics/timestepping separation. However, we believe that existing applications (e.g., flat MPI code) would already have implemented this functionality, therefore not imposing a significant burden on domain scientists. The advantage of this approach is that it allows for loop fusion. Low-order methods such as the ones considered in this article are memory bound. So requiring the entire computation of the right-hand side followed by scaling the result by the timestep and adding it to the previous state would incur additional memory traffic. Given the emphasis on performance, we have chosen the less clean but more performant option for the abstraction.

The second functionality that the problem layer must supply are the CFL-related constraints. Similar to the state updates, we partition the computation into two calls based on whether the CFL condition depends on solely local or distributed information. We refer to K as the wave speed \( |\Lambda |/\Delta x \). The local CFL bound is computed as part of compute_K_internal. From an implementation perspective, the computation of the CFL condition tends to be quite expensive. Thus, we opt to fuse the computation of the internal K with the build_u_next call and cache the result. When computing the internal K value, the only additional information that requires consideration are the domain’s boundary conditions that may depend on the time argument, t. Note that these evaluations are only at snapshots in time, the control flow associated with properly evaluating the CFL condition, Equation (4) is delegated to the timestepping layer.

When a neighbor updates and the associated push_flux is processed, we need to recheck the CFL condition. Since the neighboring state only affects a limited number of cells, we introduce an additional call compute_K_bdry, which will recompute \( |\Lambda |/\Delta \), based on the neighbor’s updated wave speed. Note that we assume that the relevant neighbor’s \( |\Lambda | \) can be computed based on the received msg_t. In our experiments, we have always found this to be the case.

We conclude with a remark on the extension of this API to higher-order timestepping methods. One approach is based on ADER methods [29]. These are single-step methods that rely on a polynomial space-time representation of the state. To implement ADER timestepping, we would only need to implement a single ADER timestep along with the relevant boundary types at the application layer. The helper function implementation could then be passed to the unmodified timestepping layer.

Alternatively, high-order strong stability preserving Runge-Kutta methods [23, 36] are also commonly used due to their non-linear stability for conservation laws. The key property of these methods is that they can be represented as convex combinations of Euler steps. Implementing this class of Runge-Kutta methods would require modifications to the timestepping layer. Advancing a submesh would require a second small subtime quantity \( \epsilon ^2 \), which would be used to identify which stage is currently being processed. Each update would cause several halo exchanges to occur at time t with the final state being updated at the update time t + dt. All of the control flow details are managed in the timestepping layer, and the forward Euler implementation of build_u_next function is sufficient for supporting this class of Runge-Kutta methods. Introduction of multi-stage Runge-Kutta methods introduces considerable complexity, but the hope is by sequestering the implementation details from users, we hope to promote adoption of these techniques.

3.2 Rolling Back the Application Layer State

The discussion up until this point has been agnostic to PDES. However, one important aspect of the API is how the application deals with rollback. The API has been designed to allow the user to optimize the amount of state required to be captured by separating data members that vary between updates (e.g., state variables) from those that remain constant (e.g., mesh geometry). Since the computations associated with build_u_next are not reversible, Devastator is required to store the entire dynamic state until global virtual time (GVT) passes the event time and the memory can be reclaimed. For this reason, we return the next state by value, which we will swap with the old state. Afterward, the current dynamic state will be stored in the Problem class and the old dynamic state will be managed by the Devastator runtime, which will deallocate it once it is safe to do so. Similarly, for the push flux event, we split apart the dynamic state of the boundary from things that remain unchanged during the simulation. By explicitly requiring these types to be defined, we enable Devastator to ensure that events remain reversible. With these dynamic types, rollback can simply be performed by swapping old dynamic types into the problem class, restoring the state prior to an event.

Last, we require that the msg_t and Problem types be serializable. Since push fluxes are sent between nodes and submeshes are migrated during load rebalancing phases, we need to ensure that the state can be serialized into a buffer and sent across the wire.

To understand the practical impact of this interface, we have collected lines of code for the various modules. Using the Problem API, we have implemented three different sets of conservation laws, which will be outlined in more detail in the next section. In addition, we have also used the code to implement a serial synchronously stepped implementation as well as an MPI-based implementation. In Table 1, we compare the lines of code for each module. Note that problem classes are specified by the equation they solve (Burgers’, Shallow Water, Euler), and the timestepping layer is labeled according to the algorithm (STS or ALTS) and how it is run (Devastator, which we abbreviate as deva, MPI, or serial). We additionally note a block of code as common. This code is used across all denoted modules and thus cannot be assigned merely to a single group. From the data, we remark that the ALTS implementation is significantly larger than the STS implementation. This is in part due to more features including performance instrumentation and load balancing capabilities as well as more complexity when compared with the STS-based implementations. However, we observe that using this abstraction layer, we achieve significant code reuse with a single problem class consisting of approximately 1,000 lines of code being able to leverage three different runtimes. In other words, across the three different conservation laws, no less than 74% of the code is reused.

4.1 Event Overheads

What makes ALTS different from traditional PDES applications is that we can freely manage event granularities by tuning the number of elements assigned to each submesh. We solve the lake at rest problem on a fixed mesh with 20,000 elements for the uniform mesh and 75,000 elements for the polynomial mesh. We selected the simulation time to be sufficiently long so that reliable measurements could be made. We set \( t^{\text{end}}=64,\,(|\Lambda |/\Delta x)_{\max }=20{,}200 \) for the uniform mesh, and \( t^{\text{end}}=1,\,(|\Lambda |/\Delta x)_{\max }=1{,}115{,}208 \) for the polynomial mesh. The simulation is run varying the number of submeshes, ranging from one to eight submeshes per worker thread. As the number of submeshes increases, the concurrency available to the discrete event simulator increases (i.e., overdecomposition), allowing useful work to occur while waiting for messages to arrive. Furthermore, each submesh must satisfy the CFL condition for fewer elements, leading to a more optimal local CFL condition. This is an important performance tuning knob for irregular algorithms, as we are able to weigh the benefits of increased irregularity versus reduced PDES overheads.

We ran the following experiments on 64 cores on Cori. Each run consists of four processes—one per NUMA domain—with 15 worker threads per process. Enumerating the submeshes from left to right, we assign submeshes in contiguous chunks—that is, rank k is assigned submeshes \( C_F k \) through \( C_F (k+1) \), where \( C_F \) is the number of submeshes per rank. The critical path will lie on the rank with the submeshes with the fewest elements. Since the work performed per submesh is balanced for the lake at rest problem, the ALTS-related work completed on this rank is comparable to that of other ranks. However, this rank must execute more updates and therefore incurs more event overheads.

In Figure 4, we plot execution time versus the lowest number of elements assigned to a submesh. Each configuration was run five times, and the standard deviation of each set of runs was within \( 0.99\% \) of the respective means. For either configuration, we see that once the element size per submesh drops below 70 elements, event overheads begin to dominate the execution time. Although there is a slight benefit to oversubscribing compute resources for the uniform mesh, for the polynomial mesh we observe a speedup of \( 1.26 \) going from an oversubscription factor of 1 (407 minimum elements per submesh) to 2 (196 minimum elements per submesh). For the uniform mesh, dependencies are approximately satisfied at the same time, but for the polynomial mesh, a neighboring submesh that is substepping must first make progress before the submesh with the larger timestep is able to progress causing the performance degradation for the largest minimum submesh in Figure 4(b). Oversubscription allows us to update a different submesh while waiting for the substepped submesh to satisfy the stalled submesh’s dependencies.

Fig. 4.

View Figure

Last, we build a cost model that approximates the time elapsed updating a submesh of \( n_{\text{elem}} \). Since the cost of the update kernel is \( \mathcal {O}(n_{\text{elem}}) \), we assume a linear model of the form \( cost(n_{\text{elem}})= C_2n_{\text{elem}}+ C_1 \), where \( C_2 \) corresponds to the cost to update one element and \( C_1 \) to roughly be the event overhead. Using the data for the polynomial mesh runs and omitting the data point for an oversubscription factor of 1 due to the aforementioned discrete effects impeding performance, we estimate the average task execution time by dividing the overall execution time by the largest number of updates committed on any rank, then use the average number of elements per submesh on that rank as the \( n_{\text{elem}} \) input. Due to the smooth nature of the mesh refinement, we assume that all submeshes on the selected rank update the same number of times. Via linear regression, we find \( C_2=0.0438\,{\mu \rm {s}}/\text{element} \) and \( C_1=1.474\,{\mu \rm {s}} \) with \( R^2=0.999 \). We note that we do not expect this cost model to generalize to larger problems. The selected problem sizes fit in the core’s local L1 cache. As elements spill out of the local caches, we expect the cost to update an element, \( C_2 \), to increase, reflecting the reduced memory bandwidth. Using this cost model as a lower bound, we approximate that roughly 100 elements are needed to achieve a good resource utilization. At this task granularity, \( 74.8\% \) of the time is spent doing useful work. This corresponds to an event size of \( 5.9\,{\mu \rm {s}} \).

If we think of the lake at rest problem as an ordinary 1D stencil dependency graph, we can compare Devastator’s event execution overheads with the large number of existing task-based runtimes using the analysis described in the work of Slaughter et al. [62]. The figure of merit in this analysis is the minimum effective task granularity (METG). Specifically, METG(50%) is the smallest task size, in our case cells per submesh, for which half the time is spent doing useful work. In our case, the METG(50%) is extrapolated to \( 3\,{\mu \rm {s}} \) (i.e., \( 2C_1 \)), making it competitive with the fastest reported implementations in the work of Slaughter et al. [62]. It is in some sense unsurprising that Devastator would do well in this metric since, like many other PDES engines, it trades off some task scheduling flexibility (e.g., local work stealing) for the ability to efficiently execute small events, where execution time is largely driven by \( C_1 \). This poses interesting performance trade-offs for the PDES and task-graph runtime communities. If the irregularity of applications is determined to be limited, optimizing for small event execution may lead to better strong scaling properties. The other interesting question is how this efficiency changes as we scale out. As shown in Figure 1, the commit rate drops as we scale out to larger node counts. Quantifying the impact of the GVT reduction on METG(50%) is an important research question we will pursue in future work.

4.2 Description of Misspeculation

To explore when and where misspeculation is occurring, we consider a single-node (world=threads) problem with a 300,000 element mesh with a 6x oversubscription factor. We use a trivial actor partitioning, leading to imbalanced work-to-rank assignments. However, as a study of speculation, this serves as a stress test since underworked ranks have the possibility to do more work that could be rolled back while waiting on messages from overworked neighbors. We run multiple configurations of the Burgers’, Euler, and shallow water equations with several initial conditions and mesh types. We report the number of element updates committed and rolled back for each configuration as well as the ratio of the two statistics in Table 3. At a high level, we remark that the waiting heuristics outlined in Section 2.3 work remarkably well. For the SWE-lake at rest problem, no rollback of updates occurred during the simulation. For this configuration, the only source of variation in timestep sizes is due to variation in element sizes, not the local wave speed, which equals 1 throughout the entire domain. Insofar, the information sent via the upper bound messages anticipates when messages will arrive and allows the submesh to avoid performing updates that would be rolled back. Otherwise, the highest observed rate of rollback was for the Euler-blast wave problem with 1.5% and 0.2% of element updates rolled back for the uniform and polynomial meshes, respectively. Beyond the blast wave problem, all other configurations had a rollback rate of less than 0.05%.

To help understand where and why rollback is occurring, Figure 5 shows the number of rollbacks per second occurring on each submesh (actor) versus wall clock time compared to the timestep sizes taken by the submeshes versus simulated time for three of the problem configurations. Since we are considering problems where the critical path goes through a single actor, the simulated and wallclock times are roughly proportional. We claim that the cause for misspeculation arises from the unanticipated changes in timestep size and the resulting rollback cascade. When a submesh is stepping with a larger timestep than its neighbor, it will try to speculate ahead of its neighbors. While the neighbor’s timestep remains stable, the submesh can use its knowledge of the neighbor’s timestep along with the upper bound mechanism (see Section 2.3.4) to avoid misspeculation by waiting on an appropriate number of messages. However, due to the evolution of the equations, that neighbor may ultimately require a smaller or larger timestep, leading to unanticipated push flux events and corresponding rollbacks.

Fig. 5.

View Figure

Consider the case of a timestep restriction where two neighbors are updating alongside one another, and one suddenly requires a smaller timestep. This is the case for the SWE dam break problem in Figure 5(a) and (b). In this problem, there are two waves emanating from the initial discontinuity for which the emergent interior region requires a more strict timestep than either side of the initial conditions. As seen in Figure 5(a), rollback occurs in elements in front of the shockwave. With the more strict timestep moving through the domain, coarser timesteps need to be rolled back and re-executed with finer timesteps. Since all cells in a submesh are updated as a unit, these rollbacks only occur when the timestep change lines up precisely with a submesh boundary.

At the coarse-fine timestep boundary, when the finely timestepped region updates, it relaxes its neighbor’s upper bound, eventually permitting a coarse timestep update in the neighboring submesh. The process of each submesh updating and relaxing its neighbor’s upper bound then propagates away from the boundary, creating a work cascade, whereby a set of tasks along the coarsely stepped light cone become executable. The speed of the work cascade propagation \( v_{w} \) is determined by the wall clock time required to execute update events and send push-flux events to neighboring submeshes. So long as the submeshes’ timesteps do not change, the work done in each work cascade is useful work that does not need to be rolled back. However, when a submesh’s timestep must be refined due to the simulated shockwave’s propagation, it induces a rollback cascade along another light cone with speed \( v_{rb} \). Since executing an update is more computationally expensive than processing a rollback, \( v_{w} \lt v_{rb} \), allowing the rollback cascade to catch up with and terminate the work cascade. We can observe this effect in Figure 5(a), where rollbacks emanate from the timestep differences into the less frequently stepped regions.

But the interesting question is this: why does the cascade emanate further from the shock front in the more coarsely stepped region of the mesh (i.e., the actors with higher IDs)? The work associated with either work cascade is the same, although the cascade with the larger timestep will advance further in simulation time. Letting \( \delta _{rb} \) be the delay between the misspeculated work cascade and the subsequent rollback cascade, the rollback cascade would terminate at a distance from the shock of \( (v_{rb} v_w \delta _{rb}) / (v_{rb} - v_w) \). We hypothesize that this delay varies based on the ratio of timestepping groups. When there is a larger difference in timesteps, there is more work to be performed at the shockwave interface that must now advance GVT at a slower rate. This leads to a longer delay before the rollback is issued, causing more rollback for the high ID actors than the low ID actors in Figure 5(a).

This speculation distance is particularly important for the Euler-blast wave problem, shown in Figure 5(e) and (f). This problem consists of two shockwaves meeting in the middle of the domain and interacting. The high actor IDs are simulating the lower energy blast, whereas the low actor IDs are simulating the higher energy blast. The two blasts are separated by a low pressure region, which is able to take a large timestep. In this case, the timestepping ratio is so large that the work cascade is able to make it across the entire vacuum state and cause erroneous events to be scheduled for the lower energy blast. Although not visible in Figure 5(e), the vacuum region does experience a rollback (\( \mathcal {O}(10^{-2})-\mathcal {O}(10^{2}) \) rollbacks per second). Since the work associated with this lower energy blast is non-trivial, most rollback is observable in the lower energy blast. We also see substantially more rollback when compared to the other two configurations in Figure 5(a) and (c). We determine that there is relatively limited benefit to speculating ahead of the high energy blast (i.e., independently advancing the lower energy blast). Rollback cascades from the misspeculation of the more computationally intensive higher energy blast undo speculative work of the lower energy blast. We also observe that the misspeculation substantially reduces once the two shocks meet and the discrepancy in timestep reduces. From the mathematical perspective, it is interesting to note that these rollback cascades are unphysical. The work done for the lower energy blast should be independent of what is happening far away. We believe that for most of the rollback cascade, the re-executed push flux has an identical state to the unexecuted push flux. One area of future work is implementing lazy cancellation [34] in Devastator to disrupt these rollback cascades.

These experiments showcase the latent irregularity in the application. We see that timestep changes are relatively infrequent and localized, although they may move throughout the domain. The benefits of speculation are less about doing useful work while being uncertain if it will need to be rolled back, but rather avoiding paying the expensive synchronization overheads required to prove that no causality violations occur.

The occurrence of speculation is fairly limited in all of these configurations. When the local prediction of incoming messages is accurate, we are able to prevent bad speculation by waiting rather than executing updates that will inevitably be rolled back. The execution of the local timestepping algorithm is equivalent to that of the task-graph-based implementations of local timestepping for linear problems for which \( |\Lambda | \) is independent of the submesh state. Furthermore, in this section, we have artificially aggravated misspeculation by running an imbalanced configuration. However, most often the misspeculation happens off the critical path. In Figure 6, we show the cumulative element updates and rollbacks for each actor for the uniform Euler-blast wave problem. Rolled back element updates are concentrated on actors that are idly waiting for updates. If we identify the rank with the most work (consisting of actors 18 through 23), we note that only \( 3.1\cdot 10^{-3}\% \) of element updates are rolled back.

Fig. 6.

View Figure

4.3 Strong Scaling Comparison to STS-MPI

We turn to distributed performance with a strong scaling study of the Carrier-Greenspan problem on a polynomial mesh with 1.2 million elements. We compare the performance of the ALTS Devastator (ALTS-deva) to a synchronously stepped MPI implementation (STS-MPI). At 1.2 million elements, the smallest submesh size is 99 elements, limiting Devastator overheads per Section 4.1. For the MPI version, each submesh steps at each timestep (i.e., synchronously). To ensure a fair comparison over the existing baseline parallelization, STS-MPI implementation utilizes point-to-point non-blocking messages and overlaps communication with available local work, consistent with previous work [11, 14]. Furthermore, both implementations use the same actor data structures and application layer helper functions.

Since ALTS improves performance via work reduction, we use a throughput measure that approximates the throughput as work done for a synchronous problem divided by the execution time—that is, \( \begin{equation*} TP = \frac{ t^{\text{end}}n_{\text{elem}}}{ \Delta t^{\min }T^{n}_X }, \end{equation*} \) where \( t^{\text{end}}=2\pi \) is the simulation time, \( n_{\text{elem}}=1.2 \cdot 10^6 \) is the number of elements in the mesh, \( \Delta t_{\min }=7.8125\cdot 10^{-8}\,((|\Lambda |/\Delta x)_{\max }=6.4\cdot 10^6) \) is the smallest timestep taken, and \( T^n_X \) is the execution time measured in seconds with n cores using parallelization X—either Devastator or MPI. Figure 7 compares the output of the MPI-synchronous timestepping (STS-MPI) implementation to the Devastator-ALTS (ALTS-deva) implementation. For the MPI implementation, the 32 core run has been omitted since it was unable to finish within queue time constraints. The core counts are chosen to highlight the scaling behavior of the ALTS-deva implementation. At 1,024 cores, the oversubscription factor is three submeshes per rank. Scaling out further without modifying submesh sizes would lead to discrete load balancing issues, degrading performance. Since ALTS achieves its speedup through work reduction, we observe an improvement in per core throughput when run on 128 to 1,024 cores rather than an improvement in the maximum number of cores to which we are able to scale out. At 128 cores, the ALTS-deva implementation has a speedup of 3.75x over the STS-MPI implementation. This speedup degrades to 2.57x at 768 cores. Nonetheless, at 768 cores, the ALTS-deva implementation achieves a time to solution comparable to the 2,048-core STS-MPI configuration.

Fig. 7.

View Figure

This strong scaling study is designed to show tangible benefits of ALTS to accelerate effective per core throughput. However, the more regular workload of the STS-MPI application remains more scalable. We have shown that the 2,048-core MPI run matches the best throughput of the ALTS-deva implementation. At 6,144 cores, the STS-MPI implementation has a parallel efficiency of 80% and a speedup over the best ALTS-deva time to solution of 2.42. This result is important to emphasize. If fastest time to solution remains the primary objective, synchronous timestepping remains the best option. However, if the problem exhibits large spatial variance in the wave speed \( |\Lambda | \), can be solved within provided time constraints (e.g., supercomputer queue time limits), and can efficiently utilize compute resources—by breaking the mesh into sufficiently large submeshes while retaining sufficiently many submeshes per core—we expect ALTS to substantially reduce the requisite computational effort to simulate the given configuration.

The speedups obtained are ultimately determined by the distribution of element sizes \( \lbrace \Delta x_j\rbrace \) and problem dynamics. Taking the ratio of element updates for the synchronous timestepping method divided by the number of element updates for ALTS, we find a theoretical speedup based on work estimates of \( S^{work}=5.38 \). The remainder of this section will focus on explaining the discrepancy between the observed speedups between the STS-MPI implementation and the ALTS-deva implementation. There are five contributors worth examining when comparing the two implementations:

Based on these analyses, we determine that the primary drivers of performance loss are computing the CFL condition and worsening load balance. Nonetheless, at low core counts, the problem can be sufficiently balanced, and the ALTS-deva implementation achieves throughput improvements of up to 4× over the STS-MPI implementation, which is 74.5% of the peak theoretical speedup.

4.4 Weak Scaling Performance

To show that ALTS can be scaled to larger distributed problems, we perform a weak scaling study using the Carrier-Greenspan problem on the uniform mesh. Since the performance of the ALTS method strongly depends on the distribution of the wave speed, \( |\Lambda | \), it is important to simulate the full \( 2 \pi \) time to compare performance across different problem sizes. Due to the fact that the amount of computational work scales quadratically with the number of elements, and the fact that the polynomial mesh is computationally too expensive, using the uniform mesh allows us to consider a larger range of problem sizes.

We consider four submeshes with 500,000, 1 million, 2 million, and 4 million elements. The problem is partitioned into 2,880 actors per 500,000 elements. For reference, storm surge simulation meshes have \( \mathcal {O}(10^6) \) elements [27]. This leads to a minimum submesh size of 174 elements, ensuring that simulator overheads are reasonably amortized. The rebalance frequency is taken to be once every 8,192 \( \Delta t_{\min } \). As the number of elements doubles, \( \Delta t_{\min } \) is halved, and the total number of load balancing phases doubles. With the work per core per epoch remaining constant across weak scaling configurations, we observed that this rebalance frequency provided the optimal performance.

Each configuration is run over a few core counts to demonstrate strong scaling characteristics. The parallel efficiency versus core count is shown in Figure 8. The simulation weak scales well with the 4 million element mesh obtaining 90% of ideal per core throughput at 1,024 cores and 68% of ideal per core throughput at 3,072 cores. Looking at the imbalance across the runs, we see that the average imbalance for the 32-core 500,000 element mesh was 0.03. For the 4 million element mesh, the imbalance increases to 0.14 for the 1,024-core configuration to 0.33 for the 3,072-core configuration, leading to a slight degradation in scalability. The choice of a uniform mesh will lead to higher fractions of elements updating each timestep than a polynomial mesh. Thus, the average number of elements per core may not be a good indicator of work required to achieve good scaling behavior, and nonetheless we conclude that if the problem is load balanced and sufficiently large, the ALTS method can be scaled out.

Fig. 8.

View Figure

4.5 Impact of Performance Optimizations

To help quantify the impact of the performance optimization strategies described in Section 2.3, we consider a few problems and then enable the strategies to observe their impact. This section is split into two sections highlighting the benefits of waiting on misspeculation and load balancing.

4.5.1 Waiting Heuristics.

To determine the impact of the performance optimizations outlined in Sections 2.3.1 through 2.3.4, we incrementally enable them and report the impact on speedups relative to a local timestepping baseline with no waiting heuristics enabled. We consider a 512-core lake at rest problem (\( n_{\text{elem}}=6\cdot 10^5,\,(|\Lambda |/\Delta x)_{\max }=1,\!370,\!368,\!000,\,t^{\text{end}}=0.1,\, n_{\text{actors}}=1,\!440 \)) and 768-core lake at rest problem (\( n_{\text{elem}}=1.2\cdot 10^6,\,(|\Lambda |/\Delta x)_{\max }=2,\!740,\!736,\!000,\,t^{\text{end}}=0.05,\,n_{\text{actors}}= 1,\!440 \)). The \( (|\Lambda |/\Delta x)_{\max } \) parameters have been chosen to be intentionally large—that is, \( (|\Lambda |/\Delta x)_{\max }\approx 256/\Delta x_{\min } \). To avoid exaggerating the impact of the waiting heuristics, these problems have been sized to be sufficiently large to sustain a reasonable parallel efficiency. The lake at rest configurations achieve 90% of the throughput rate of lake at rest problems with twice the number of elements. We also consider two 512-core Carrier-Greenspan configurations: one that exactly mimics the strong scaling configuration and the second the same as the first except load balancing is disabled. For the Carrier-Greenspan configurations, \( (|\Lambda |/\Delta x)_{\max } \) is chosen to be close to the optimal value. The cumulative speedups for each configuration relative to a baseline without any optimizations enabled are presented in Figure 9. For these bar graphs, cumulative implies that the preceding optimizations remain enabled (e.g., the speedup for Section 2.3.3 has also enabled Sections 2.3.1 and 2.3.2).

Fig. 9.

View Figure

For the lake at rest problems, we observe the largest benefits occur when binning timesteps to powers of 2 (Section 2.3.1) and forcing neighbors to update (Section 2.3.3). The large \( (|\Lambda |/\Delta x)_{\max } \) values cause a substantial number of local ordering updates to occur as described in Section 2.3.1 and impede the submesh partitioner, as the expected number of element updates is not well approximated by the local CFL condition. Introducing timestep binning (Section 2.3.1) elides local ordering-related updates with CFL-condition-related updates and causes the number of element updates committed per rank to become roughly equal. For the 512-core lake at rest problem, this corresponds to the most overworked rank committing 1.67x the average number of element updates before timestep binning is introduced and 1.01x the average number of element updates after timestep binning is introduced. Forcing neighbors to update (Section 2.3.2) and waiting on those updates (Section 2.3.2) eliminates most misspeculation, causing a second increase in improvement. For the Carrier-Greenspan problem, the close-to-optimal \( (|\Lambda |/\Delta x)_{\max } \) reduces the extra updates associated with enforcing local ordering. In that case without Section 2.3.1, timesteps are binned into multiples of \( \Delta t_{\min } \) (rather than \( 2^8 \Delta t_{\min } \)). In conjunction with a more optimal timestep taken at the critical path, we observe that the time to solution slightly increases after binning timesteps into powers of 2. Similar to the lake at rest problem, a large performance improvement is only observed after the first three optimizations have been enabled.

Another takeaway from these results is that misspeculation impacts time to solution only when it delays the critical path. Since rollback is triggered by arrival of delayed push fluxes from neighboring submeshes, when the simulation is very imbalanced, the rank which is on the critical path has the slowest rate of advancing the simulation time, and as this rank makes progress, messages that would incur rollback have most likely already arrived. The unbalanced Carrier-Greenspan configuration with no performance optimizations enabled executed 1.3x more element updates than the configuration with all optimizations enabled yet ran 10% faster. Only when considering load balanced problems, either through construction (lake at rest) or via load balancing (Carrier-Greenspan), do we observe performance improvements by reducing misspeculation. In these cases, the rollback-causing messages have not necessarily arrived, and anticipating and waiting on these messages, when feasible, leads to improved time to solution. For the 768-core lake at rest problem, the percentage of the cell updates rolled back on the critical path goes down from \( 46\% \) down for the Section 2.3.1 configuration to zero for the fully optimized case.

4.5.2 Load Balancing.

For the 512-core strong scaling problem, load balancing plays a significant role, delivering a speedup of 1.8x over a run without dynamic load balancing. In Figure 10, we compare the impact of rebalancing frequency on time to solution. We note that the problem has a relatively wide minimum. With only 614 epochs, the time to solution is only 13% larger than the optimal time to solution, which occurs with 9,818 epochs. To understand this behavior, note that the worst case imbalance corresponds to when all submeshes on a rank are wet. In that sense, once the rebalance frequency begins to significantly exceed the rate of flooding, the problem becomes effectively balanced.

Fig. 10.

View Figure

Thus, it is important to emphasize that while the algorithm is and must be stable under very rapidly changing values of \( |\Lambda | \), these scenarios are rarely seen in practice. Rather, we see gradual shifts in timestepping group membership, which if not accounted for will lead to instabilities but are nevertheless gradual enough that these simulations can be load balanced effectively. Although we have only considered a 1D problem, the simulated dynamics—local refinement in the middle of the domain and flooding moving into the refined region—represent a worst case scenario for hurricane storm surge modeling. Rarely would all inundation be concentrated in a single submesh; rather, inundation would be spread out over a large front, lowering the overall rate of imbalance generation.