A Lagrangian meshfree method applied to linear and nonlinear elasticity

The repeated replacement method (RRM) is a Lagrangian meshfree method which we have previously applied to the Euler equations for compressible fluid flow. In this paper we present new enhancements to RRM, and we apply the enhanced method to both linear and nonlinear elasticity. We compare the results of ten test problems to those of analytic solvers, to demonstrate that RRM can successfully simulate these elastic systems without many of the requirements of traditional numerical methods such as numerical derivatives, equation system solvers, or Riemann solvers. We also show the relationship between error and computational effort for RRM on these systems, and compare RRM to other methods to highlight its strengths and weaknesses. And to further explain the two elastic equations used in the paper, we demonstrate the mathematical procedure used to create Riemann and Sedov-Taylor solvers for them, and detail the numerical techniques needed to embody those solvers in code.


Introduction
The repeated replacement method (RRM) [1] is a Lagrangian meshfree method for the simulation of time-dependent systems of conservation laws. We previously used RRM to simulate the one-dimensional Euler equations for compressible fluid flow to establish the basic functionality of the method, and we compared the results to an exact Riemann solver for the Euler equations given by Toro [2]. In this paper, we enhance RRM to increase its speed and accuracy, and apply it to the more challenging problems of one-dimensional linear and nonlinear elasticity. This allows us to demonstrate that RRM works for a range of constitutive equations, while maintaining good accuracy and error scaling behavior.
In this paper, we first motivate and derive our linear and nonlinear elastic constitutive equations, to define the terms and symbols used in subsequent sections. Second, we give anoverview of the repeated replacement method with a detailed comparison to earlier work, highlighting how the strengths and weaknesses of RRM differ from those of previous methods. Third, we explain the improvements made to RRM in the course of adapting it to handle elastic systems. Fourth, we present the derivation and explanation of our Riemann and Sedov-Taylor solvers. Then we show the results of RRM for eleven test problems, validating the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 results against Riemann and Sedov-Taylor solvers, and show that RRM's error convergence behavior is somewhat super-linear. Finally, we present a summary and conclusion.

Elastic equations
For both linear and nonlinear cases, we assume a homogeneous, one-dimensional, elastic continuum material that is subject to finite strain, which is a strain large enough that we cannot simplify our mathematical treatment by assuming it to be infinitesimal. The material's reference (or undeformed) density is ρ 0 , and its spatial (or deformed) density at each point x is ρ(x). We will express the deformation in terms of the stretch λ, which is written in terms of density as λ(x) = ρ 0 /ρ(x). For a finite-sized piece of material, stretch is defined as λ l/l 0 , where l 0 and l are the undeformed and deformed lengths of the piece, respectively. Strain is defined as (l − l 0 )/l 0 = λ − 1.
We assume that the internal energy density of the material due to its state of strain and temperature, per unit of reference length, can be expressed as some function C(λ, T) = C s (λ) + C t (λ, T), where C s is the strain energy density, and C t is the thermal energy density.
The Cauchy stress in such a material has both an elastic and a thermal component, and can be derived from the energy density as sðl; TÞ ¼ s e ðlÞ þ s t ðl; TÞ ¼ @C s ðlÞ @l þ @C t ðl; TÞ @l ð1Þ A material whose stress is derived from an energy density function in this way is called a Green-elastic or hyperelastic material.

Linear elasticity
Linear elasticity is defined by the direct proportionality of stress to strain. The elastic part of the stress is where E is the elastic modulus, which has the units (kg Á m)/s 2 in one dimension. To obtain the strain energy density, we note that the integral of stress over distance is strain energy, and strain energy density with respect to the reference length is simply strain energy/l 0 . So, making use of the fact that λ = l/l 0 , the strain energy density is To complete this simple model, we choose a thermal energy density C t (T) = C v T where C v is the heat capacity in J/(m Á K). This C t has no dependence on λ, since it is difficult to form such an expression that has the correct signs for both energy and stress equations without making the stress nonlinear in λ. This results in a total stress of Note that this stress does not model thermal expansion, and also allows the material to be compressed down to zero length with finite work C s (λ)| λ = 0 = E/2, so we have sacrificed physicality for the sake of linearity in this simple model. The total energy density in spatial coordinates, including kinetic energy, is C tot ðl; u; TÞ ¼ ru 2 where u is velocity, and the division of the C term by λ is needed to convert the energy density from reference to spatial coordinates. Finally, making use of the fact that λ = ρ 0 /ρ, the speed of sound in the material in terms of the density is aðrÞ ¼ ffiffiffiffiffi ffi @s @r Nonlinear elasticity There are many existing models of nonlinear elasticity. Some, such as that of Ogden [3], are complex enough to form an energy density from any desired polynomial in λ to better match the behavior of real materials. In our case, we are merely trying to test a numerical scheme, not model a specific material, so instead of adopting a general model, we modify the linear elastic equation only as much as needed to prevent manifestly unphysical behavior. Specifically, if we change the −1 in the linear elastic stress to −1/λ 2 , we obtain This new nonlinear stress approaches −1 as the stretch approaches zero, which prevents the material from being compressed to zero length with a finite amount of energy. Integrating this stress the same way we did for linear elasticity, we obtain the strain energy density C s ðlÞ ¼ E 2l ðl À 1Þ 2 ðl þ 2Þ.
To model the fact that real materials expand when heated, we add a thermal energy density term C t (λ, T) = C v T/λ. The division by λ is because we want the sign of this term to become negative in the stress equation to reflect this expansion.
The total stress σ is again obtained from the energy density function by sðl; TÞ ¼ @C s ðlÞ @l þ @C t ðl; TÞ @l and the total energy density in spatial coordinates is again given by Finally, the speed of sound in the material in terms of the density and temperature is aðr; TÞ ¼ ffiffiffiffiffi ffi @s @r Note that the dependence of the thermal stress σ t on λ is what causes the T term to appear in the speed of sound, without which the material would only support shock waves and not rarefactions. We demonstrate this in the course of deriving our Riemann solvers.

Overview of the repeated replacement method
RRM is based on the observation that when we model conservation laws using a field of piecewise-constant cells, the primitive values in the cells' interiors do not change directly. Instead, wavefronts of change expand out from the edges between adjacent, dissimilar cells, where the spatial derivative is nonzero. These expanding wavefronts carry primitive value changes into the cells' interiors.
To motivate this statement, consider the following. As we will see in the section on the derivation of Riemann solvers, we can write a system of conservation laws in the form W t + A(W) W x = 0, where W denotes a vector of primitive variables, A(W) is a Jacobian matrix derived from the constitutive equations, and the subscripts represent partial differentiation with respect to t and x. From this equation, we can see that if the spatial derivative W x of the field is zero at a point, then the temporal derivative W t of the field must also be zero at that point. Since our cells are piecewise constant, the spatial derivative of the field is zero everywhere except at the edges between cells, so time evolution can only begin at these edges.
In RRM, we track the expansion of these wavefronts over time. Once a wavefront reaches a certain width, we chop out the cells and parts of cells it encompasses, and replace them with a new piecewise-constant cell containing the same mass, momentum, and energy. The more different the adjacent cells are, and the lower the user has set the maximum allowable error, the less a wavefront is allowed to expand before it is replaced by a new cell. Each new cell gives rise to two new wavefronts, one at each edge.
To process these wavefronts, RRM uses an event-based simulator. Wavefronts are stored in a queue, ordered by the wavefronts' replacement time. Each simulation event consists of updating the simulation time to the replacement time of the soonest wavefront in the queue, removing that wavefront from the queue, using it to create a new cell in the field, and adding the new cell's two new wavefronts back into the queue. Continuing this process drives the time evolution of the field. The following sections describe this process in greater detail.

Cells, wavefronts and tracer particles
RRM subdivides the domain into cells c i , each of which contains constant values of the primitive variables W = [ρ u T] T . A simple two-cell set of initial conditions is shown in Fig 1. At the start of the simulation, the cells are typically chosen to have no gaps or overlaps between them. But since this is a Lagrangian meshfree method, the cells are what moves, not a mesh, so gaps and overlaps between cells will appear and disappear over the course of the simulation.
If we write a system of conservation laws as W t + A(W)W x = 0, we can see that @W @t ¼ 0 if @W @x ¼ 0, so no evolution takes place inside constant-valued cells. Instead, all evolution originates at the cell edges where the spatial derivative is nonzero. For the edges between c 1 and c 2 , we trace out an expanding wavefront w 12 using a pair of tracer particles p 1r and p 2l moving at the speed of sound in c 1 and c 2 , respectively, as shown in Fig 2. We call these "tracer" particles because they are just an aid to simulation, rather than a representation of real physical particles.
For the tracer particles, we define an error metric which gives a measure of the distance-weighted total variation in primitive values which the tracer particle encounters as it moves from cell c 1 to cell c n , where cell c 0 is the cell on the opposite side of the wavefront where the tracer particle was created, d i is the distance traveled by the tracer particle inside cell c i , and W i is the vector of primitive values for cell c i . Fig 3 shows what that same wavefront expansion looks like as a spacetime diagram. In this format, we make the cells implicit as white areas, and show wavefronts as colored triangles Expanding wavefront. Wavefront w 12 and tracer particles p 1r and p 2l expanding from the edges between cells c 1 and c 2 . Tracer particle p 1r moves at speed a(W 1 ), and tracer particle p 2l moves at speed a (W 2 ). The left and right bounds of the domain are inactive in this example, otherwise we would show wavefronts expanding from them as well.
https://doi.org/10.1371/journal.pone.0186345.g002 expanding in space as the time increases going downward. This format will allow us to represent many hundreds of events succinctly on one diagram.
Once a wavefront gets so wide that the error metric of one of its tracer particles exceeds some user-defined limit Δ max , we chop out the cells and parts of cells under the wavefront, accumulate all their conserved quantities, and replace that area with a new constant-valued cell. This "flattening" process is shown in Fig 4. If our simulation includes multiple material types, the flattening process may produce multiple new constant-valued cells, one for each connected area of a single material type. The same flattening process is shown as a spacetime diagram in Fig 5. Conversion between momentum, kinetic energy, and potential energy is modeled using the stress momentum P σ = σΔt and stress energy E σ = σuΔt, which we store in each cell in addition to the primitive values. Each cell c i initially contains to the right, where w i is the width of the cell and a i = a(W i ) is the speed of sound in the cell. These quantities sum to zero for each cell, so they do not affect the overall amount of momentum and energy in the simulation. Tracking P σ and E σ in this way allows conservation of momentum and energy to be checked at every simulation event, even though the entire field is never updated simultaneously. Note that P σ and E σ are not precisely fluxes, since they do not cross cell boundaries, but they do serve a similar purpose in that they drive the time evolution of the field.
During flattening, each cell edge inside the wavefront contributes P σ = ±σΔt and E σ = ±σuΔt to the flattening process, where Δt is the amount of time since the cell was created, or since that side of the cell was last chopped. The total mass, momentum and energy M, P, E accumulated by one wavefront chop of cells c 1 through c n by wavefront w 1n are where M c 1c and M c nc are the partial masses chopped from c 1 and c n , respectively, by the wavefront, and M c 2 through M c n − 1 are the entire masses of cells c 2 through c n − 1 which are subsumed by the wavefront. A similar notation holds for the momenta P and energies E, using the energy density C tot from the elastic model for the chopped energy. So for example P σc 1 r + P σc 2 l is the net stress momentum across the edge between cells c 1 and c 2 , and E σc 1 r + E σc 2 l is the net stress energy across the edge between cells c 1 and c 2 . Once we have accumulated the conserved quantities, we can determine the primitive values of the new cell algebraically by where w w is the width of the wavefront at replacement time, KE and PE are the kinetic and potential energies of the new cell, respectively, and C s and C t are taken from the appropriate elastic model.
In the case where more than one type of material is present inside a wavefront during flattening, we create a new cell for each contiguous area of each material type, and flatten those areas separately, with the exception of stress momentum and stress energy, which are divided among the new cells in proportion to their size.

Early replacement and wavefront merging
When two adjacent cells have very different primitive values, the wavefront created across their shared edges will be replaced quickly, resulting in the creation of a relatively small new cell and an overall increase in the number of cells in the simulation. This is called early replacement, and is the case shown in Fig 4. In areas where cells' primitive values are similar, adjacent wavefronts will intersect before the error metric grows large enough to require replacement. To reduce the number of tracer particles we must track during simulation, we do not allow wavefronts to overlap, so we merge intersecting wavefronts together unless their summed error metrics would exceed the userdefined limit Δ max . This wavefront merging results in a decrease in the number of particles and cells in the simulation, and is illustrated in Fig 6. Any number of wavefronts may be merged together in this way, so long as the total error metric remains lower than Δ max . Fig 7 shows an entire cycle of wavefront expansion, flattening, new wavefront creation, and wavefront merging in one spacetime diagram. We will use this more succinct format when we show long sequences of events.
The interaction between early replacement and wavefront merging is what makes RRM adapt to local conditions across the field. And since replacement and merging are scheduled by event-based simulation instead of a global time tick, greater activity in one area need not affect other areas. The current implementation of RRM uses a single event queue, but events are only dependent if their wavefronts are spatially adjacent, so this event queue could be broken up into separate queues for different regions of the field to parallelize the algorithm.
An intentional result of wavefront merging is that it allows RRM to preserve a constantvalued field exactly. In other words, it allows RRM to satisfy the geometric conservation law (GCL), so named by Thomas and Lombard [4] and more recently discussed by Guillard and Farhat [5]. We can see that in a constant-valued field, stress momentum P σ and stress energy E σ across any pair of cell edges will sum to zero, since the cells' primitive values are equal. So then Eqs (12) and (13) tell us that if we chop out any wavefront in a constant-valued field, it will encompass an amount of material which, upon flattening, will yield a new cell with the same primitive values. However, this is only true if wavefronts are not allowed to overlap. If they do overlap, it is possible for a wavefront in a constant-valued field to encompass unbalanced stress momentum and stress energy, unless all the overlapping wavefronts are unioned and replaced as one, as mentioned in the section on improvements to RRM.

Positivity preservation
We say that a scheme is positivity-preserving if it is guaranteed never to produce a negative value for a quantity which should not be negative, such as density, pressure, or temperature. For density, RRM is always positivity-preserving, since new cells' masses are added up from chopped parts of other cells. And for the Euler equations, though we occasionally see a negative pressure upon flattening, those negative pressures are the result of numerical error and are on the order of the machine precision. But for nonlinear elasticity, negative temperatures can occur because of the form of the energy density Eq (9), which is dependent upon both stretch λ and temperature T. In the Euler equations, once a cell's size, mass, and velocity have been determined, any remaining energy can be assigned to potential energy simply by setting the pressure p to the necessary value. But in nonlinear elasticity, for a specific cell stretch and mass, a specific amount of strain potential energy is required, which may not leave enough energy for the thermal potential energy at that same stretch. The result is a negative temperature.
A positivity analysis shows that this effect is most pronounced at low temperatures, when neighboring cells are very different in density and velocity. Interestingly, the effect does not depend on how long a wavefront has been allowed to expand, so it cannot be remedied simply by increasing the temporal precision. Instead, the effect seems to be inherent in the constitutive equation itself, which may point to an additional constraint we should impose during the construction of such equations.
In the test cases we examined, negative temperatures only occur near contacts of dissimilar materials, since the lack of heat diffusion across the contact restricts the energy available to each flattened cell. The solution is to add a small amount of thermal energy to bring the flattened cell up to T = 0, and then to maintain energy conservation by adding a corresponding negative strain energy to the flattened cell. This negative energy is carried forward a short time into the future, until sufficient energy is encompassed by nearby wavefronts to cancel it out. We use the same solution for the smaller numerical errors. Fig 8 shows the processing of events during simulation. The initial events are those enqueued for the wavefronts spanning the edges of the initial condition cells. The event queue is arranged in order of increasing wavefront replacement time, so the soonest event is always the next one in the queue, though new events may be inserted at any position.

Flowchart
To keep the flowchart readable, we omit a few optimizations that are present in the current implementation of RRM. For example, if two wavefronts intersect, but their combined error would be over the limit, they are chopped out into two adjacent new cells instead of being merged. And at the edge between the two new cells, only one new wavefronts is created instead of two, to avoid two identical new wavefronts on top of each other. Comparison to previous work RRM was designed for use on systems of conservation laws that are mathematically inconvenient, so it requires only that the system have conserved quantities that can be numerically integrated, and that some expression for the speed of sound can be obtained. RRM does not require the evaluation of a Riemann solver, so it may be applied to conservation equations for which the Riemann solver derivation procedure is intractable. RRM also does not require the evaluation of spatial or temporal derivatives, so it may be applied to equations where such derivatives are expensive or not well defined at every point.
RRM was also designed to be adaptive. Temporally, it is unconditionally stable for time steps of any length, though it loses accuracy as the time step length is increased. Spatially, it does not require a mesh, and it can add or remove cells wherever they are needed to maintain accuracy, without the requirement to stitch together or nest meshes of different resolutions.
Of course, advantages generally come paired with disadvantages. RRM's main disadvantages are its high computational intensity, high programming complexity, and low-order reconstruction of solutions.
The largest difference between RRM and methods such as the finite difference (FD) method, finite volume (FV) method, and the finite element method (FEM) is that RRM does not assemble a system of explicit or implicit simultaneous equations for some area of the field and then apply a solver to find the vector of unknown primitive values for that whole area at each time step. Instead, RRM is a purely local method, where wavefronts are added and removed only as required to meet the user's accuracy goals, and the temporal and spatial steps can be different for every wavefront. Another major difference is that RRM does not explicitly bound its time steps using the Courant-Friedrichs-Lewy (CFL) condition, and such a condition is not required for stability. RRM does track the intersection of wavefronts, but if the intersecting wavefronts are sufficiently similar, they may simply be merged together, which causes a loss of accuracy, but does not affect stability. Since wavefronts travel at the speed of sound, their intersection generally happens at a time later than the "signal time" used in codes like Whitehurst's FLAME [6], which is based on the time taken for opposing shocks to traverse a cell and intersect within it. RRM cannot use the signal time because it does not solve the Riemann problem across the cell edges, so it does not know the shock speeds. However, RRM does track the motion of cell edges, and in cases where two cells interpenetrate at a speed faster than the speed of sound, the wavefront emanating from that pair of edges is broadened so that it always encompasses them.
RRM is probably most similar to cell-centered Lagrangian finite volume schemes such as those described by Godunov [7], then later by Després and Mazeran [8], Maire et al. [9], Maire [10], Carré et al. [11], Kluth and Després [12], Burton et al. [13] [14], and Vilar et al. [15], among others. The main difference between these schemes and RRM is that instead of changing the primitive values and shapes of existing cells based on momentum and energy fluxes across the cell faces, in RRM cells are unchanged after their creation, except for having parts chopped off to form new cells. In RRM the stress momentum P σ and stress energy E σ play a role similar to momentum and energy flux in Lagrangian finite volume schemes, but strictly speaking they are not fluxes, since they do not cross from one existing cell to another. Instead, they are one source of the conserved quantities used to form new cells, the other source being the chopped-off parts of other cells.
As a Lagrangian method, RRM automatically achieves Galilean invariance, whereas in some Eulerian methods the bulk velocity of the material may alter the simulation results, as noted by Wadsley et al. [16]. And unlike other adaptive methods such as adaptive mesh refinement (AMR) [17], RRM does not require cell sizes or time steps at different refinement levels to be multiples of each other, and does not require grids of varying degrees of refinement to be nested in any specific way.
Some variants of the finite element method such as those of Peraire et al. [18] and Radovitzky and Ortiz [19] support fully Lagrangian operation. These methods use adaptive and frequent remeshing to prevent the mesh from becoming tangled in areas of high deformation. The remeshing methods used, such as the advancing front algorithm, are paired with refinement indicators that attempt to evenly distribute quantities such as deformation among the elements of the mesh. This permits simultaneous coarsening and refinement of the mesh in different locations. However, remeshing the entire field is an expensive operation which must be carefully designed to minimize impact on simulation speed, and the transfer of conserved quantities from the old to the new mesh can introduce numerical diffusion. Such diffusion may be minimized by arbitrary Lagrangian-Eulerian (ALE) methods (originated by Hirt et al. [20] and recently reviewed by Barlow et al. [21]), which allow the mesh to move independently of the material only where required to fix mesh tangling.
Moving-mesh codes which are based on Delaunay or Voronoi tessellations of a set of moving points (such as FLAME [6], Springel's AREPO [22], Duffell and MacFadyen's TESS [23], and that of Gaburov et al. [24]) are similar to RRM in that a set of freely-moving objects serves as the basis of the field. However, these tessellation-based methods may require some care to insure that the point positions do not imply a degenerate mesh, and they may need to add points, fuse points, or steer the trajectories of the points to keep them near the centers of mass of their cells, thereby keeping the cells rounder and more tractable. Since RRM does not compute any derivatives, it does not require any constraints on cell motion or position to avoid degeneracy. However, RRM does need to prevent cells from being subdivided down to a size too close to machine precision, to insure that the edges of the cells remain distinguishable from each other when subjected to numerical error.
RRM is similar to AREPO [22] in that it can incrementally add or remove cells in specific areas of the field. In AREPO, these processes are called refinement (when a cell is split in two) and de-refinement (where a cell is removed and its conserved quantities distributed among its neighbors). The difference is that in RRM, cells are never enlarged, they are either nibbled away by neighbors over multiple simulation events, or are removed all at once after their left and right wavefronts intersect inside the cell.
RRM ensures the "manifestness" of conserved quantities by storing two extra quantities in each cell, the stress momentum P σ and the stress energy E σ . As noted in the section on wavefronts and particles, these quantities insure that after each event, global momentum and energy are still exactly conserved, even if, for example, the event introduces some change in momentum or energy which has not yet been balanced by a corresponding event elsewhere in the field. Other codes such as AREPO [22] and Hopkins' GIZMO [25] solve the same problem by requiring time steps to conform to a power-of-two hierarchy so that time steps at a given level of the hierarchy may be synchronized without synchronizing the entire field, and by requiring fluxes on both sides of a face to be updated in the same time step.
RRM is similar to discontinuous Galerkin (DG) methods [26] in that there is no global enforcement of continuity across cell boundaries as there is in the finite element method. Discontinuous Galerkin methods do assemble a matrix of equations which must be solved, but this matrix is typically block diagonal, so it may often be solved more simply than in the case of FEM.
RRM is somewhat similar to the original Eulerian version of the Godunov method [27], in that they both represent the field as piecewise-constant values, and both explicitly model time evolution at cell edges. RRM is simpler in one way, because instead of solving an exact or approximate Riemann problem at each set of edges, it only models a wavefront of two characteristics expanding outward at the speed of sound. For shock waves, Godunov's method directly models the shock speed by recovering it from the underlying Riemann solver. As we will see in the results for a nonlinear shock tube, RRM achieves the correct shock speed as a dynamic phenomenon whereby a thin velocity and stress spike forms at the shock, but without explicitly calculating the shock speed. RRM is also Lagrangian instead of Eulerian, which makes it more complex to code than the original Godunov method since cells can move and overlap.
Compared to Riemann-solver-based high-order reconstruction schemes such as MPWENO (as used for nonlinear elasticity for example by Barton et al. [28]), RRM does not require slope limiters, multiple weighted stencils, or other special care to avoid oscillation near steep gradients, since its solution reconstruction is currently only piecewise constant. RRM can also be used in cases where solution of the Riemann problem is prohibitive, either mathematically or computationally. However RRM does suffer from a one-way overshoot at shocks, as we will see in the results for a nonlinear shock tube, so it could be argued that RRM could also benefit from some special treatment to avoid creating new extrema of the solution in that case. In principle, the cells of RRM could be made higher-order in a way similar to finite volume schemes, but we have not yet investigated this possibility.
RRM bears many similarities to other meshfree methods such as smoothed-particle hydrodynamics (SPH) [29] [30] and more recent kernel-based methods such as those of Lanson and Vila [31], Gaburov and Nitadori [32], and Hopkins [25]. The cells of RRM could be compared to the particles of SPH, with the wavefronts of RRM playing a role similar to SPH's kernel functions. In SPH, for a locally-supported kernel function, the value of the field at a given point is determined by the weighted contributions of some number of nearby particles. In RRM, the "support" of a wavefront is the set of cells and parts of cells that it encompasses, but wavefronts never overlap, and the primitive values within a given wavefront depend only on that wavefront. RRM's maintenance of connections between wavefronts and cells is similar to certain acceleration techniques used in SPH implementations to reduce the number of particles that must be examined to produce the field values at each point.
The Finite Mass Method (FMM) published by Gauger et al. [33] and Klingler et al. [34] is a meshfree method that uses extended pieces of material in a way somewhat similar to the cells of RRM. FMM divides the material being simulated into finite-sized "mass packets", which are modeled by differentiable functions with compact support such as cubic B-splines. FMM packets can interpenetrate, as cells can in RRM, but FMM packets can also deform and change in size, whereas in RRM the cells are of fixed size and are chosen to be piecewise-constant specifically so their shape does not change over time. In FMM, packets are coupled by frictional forces, where in RRM the cells are independent of each other. The solution procedure of FMM is similar to that of FEM, where we assemble a system of differential equations containing the packet locations, deformation matrices, and entropies, and solve the system with standard techniques, as opposed to the event-based simulation of RRM. Finally, since FMM packets can change shape, they can become degenerate, in which case they must be recreated and the simulation continued, as noted by Klingler et al. [35].
RRM was partly inspired by work from computational geometry, specifically Chaikin's corner-cutting algorithm [36] and Catmull's subdivision surfaces [37]. Both algorithms start with simple curves and progressively refine them by performing some repeated, local operations. RRM adds the notion of conserving mass, momentum, and energy during the process, and can de-refine as well as refine, but shares a similar basic idea.

Improvements to the repeated replacement method
There are several differences between the original version of RRM and the one described in this paper. Originally, RRM was demonstrated only for the Euler equations. In the process of adapting RRM to handle elasticity, we made a number of improvements.
The change that improved accuracy the most was the addition of support for multi-material fields. In its original form, RRM is diffusive across contacts, since the chopping and replacement process mixes together material which was initially at different temperatures. An example of this behavior can be seen in the double rarefaction test in the results section. Arguably this is physically accurate, since it models heat diffusion, but it is at odds with the analytic solutions we obtain from Riemann solvers, which are adiabatic. So to make RRM adiabatic, we added the option to mark cells as being of different material types. For initial conditions with two different temperatures, if we use a different material type for each temperature, it allows RRM to avoid mixing them together. This requires cell flattening to create as many new cells as there are contiguous areas of each material type inside the wavefront. So flattening across a contact will result in two adjacent cells, which keeps contacts sharp and removes the major source of error versus a Riemann solver. Most of the other tests in the results section make use of this feature. The error analysis section contrasts the error behavior with and without this feature to quantify the improvement it brings.
Another change from the original version of RRM is in the improved handling of wavefront intersection. Originally, tracer particles were allowed to cross each other and move from one cell to another. Then at replacement time, all the wavefronts that overlapped the replacing wavefront were unioned together and replaced as one. This unioning was needed to prevent wavefronts from sometimes chopping out an unbalanced amount of stress momentum and energy, which could cause large negative temperatures or violations of the geometric conservation law, as mentioned in the sections on wavefront merging and positivity preservation. But such unioning presented problems for shock-only systems like linear elasticity, where indiscriminate unioning of wavefronts could destroy the shocks. It was also difficult to implement particle motion between cells properly for cases where multiple cells overlapped or interpenetrated at high speeds. So we changed to a system where wavefronts never overlap, and must be explicitly merged together when they intersect. This gives us the option of wavefront replacement before merge to preserve a shock or to keep error below the user threshold. It also removes the requirement to track particle motion between cells or to support more than two particles in a single cell.
The change that improved the speed of RRM the most was the reduction of the computational order of the wavefront-cell intersection operation. The previous implementation of RRM was OðN 2 Þ for this operation, because at replacement time it would intersect the wavefront with every other wavefront to check for overlaps, and with every cell to determine which cells were being removed from the field. The current implementation is OðNÞ, because a replacing wavefront never overlaps another, and because we keep track of which cells each wavefront has crossed so we can remove them from the field without intersecting wavefronts against cells.
To demonstrate the speed improvement due to this change, we profiled the performance of the same test problem for the old and new versions of the simulator. Since the old simulator only supported the Euler equations, we chose a shock-tube problem with initial conditions (ρ l , u l , p l ) = (1.0, 0.0, 100000.0), (ρ r , u r , p r ) = (0.01, 0.0, 1000.0), with the ratio of specific heats set to γ = 1.4. The simulation domain is x 2 [−1.0, 1.0], and the simulation time runs from t = 0.0 to t = 1.5 seconds. The precision settings were chosen such that the number of cells would exceed 1000, and were set the same in both versions of the simulator. The solution consists of a left rarefaction, a contact, and a right shock. Both tests were run on a single core of an Intel Core i7-3820 running at 3.6 GHz.
First we check the number of cells used over the course of both simulations, to be sure our performance numbers will be directly comparable across simulator versions.  that follow should reflect the true differences in the simulation algorithm, rather than some other change that alters the amount of work the simulator performs. Fig 10 shows the increase in simulation rate that the new version of the simulator achieves over the old version. We define simulation rate as the amount of simulation time we can process per unit of wall-clock time. The speedup is the ratio of the simulation rate of the new version to the simulation rate of the old version. Early in the simulation, when the number of cells is still small, both simulators are fast, and the new version shows less than a 2× speedup. But by the time the number of cells exceeds approximately 600, the new simulator is 6× as fast as the old one, and the speedup continues to increase over the course of the simulation as the number of cells increases.
A more human-centric measure of simulation performance is how the rate of wall-clock time increase changes over the course of the simulation. Fig 11 shows that as the simulation progresses, the amount of wall-clock time required to process a unit of simulation time goes up superlinearly for the old version of the simulator, resulting in the perception that the simulation is running more and more slowly. Some slowdown is to be expected of course, since the number of cells is increasing, which requires more computational effort. But the superlinear nature of it makes this expected behavior feel worse. In contrast, this quantity goes up only linearly for the new version of the simulator.
The final performance measure we explore is how much wall-clock time is required to process each simulation event.

Derivation of Riemann solvers
To validate our simulation results for both linear and nonlinear elasticity, we derive solvers for the Riemann problem of left and right constant states W L and W R separated by a discontinuity at x = 0. The primitive variables we use are either W = [w 1 w 2 w 3 ] T = [ρ u T] T or [ρ u σ] T , where ρ is density in spatial coordinates, u is velocity, T is temperature, and σ is Cauchy stress, as discussed in the elastic equations section. We will use one set or the other of primitive variables depending on which is more mathematically convenient for a given part of the derivation, but we will always express initial conditions in terms of temperature instead of stress, because we find the temperature values more intuitive.
In this section we use the theory of quasi-linear systems of hyperbolic partial differential equations, and generally follow the notation of Toro [2] and more recently Titarev et al. [38] But this section is meant only to outline the application of this procedure to these specific elastic systems, not as a rigorous derivation of the procedure itself, so the reader is advised to consult the cited references for more detailed information.  Introduction to the Riemann solver derivation procedure For both the linear and nonlinear elastic systems described above, we will initially write the conservation equations in differential form where the conserved variables U and the fluxes F(U) of mass, momentum, and energy are U ¼ where definitions of C tot (λ, u, T) and σ(λ, T) are given in the sections on linear and nonlinear elasticity, and we use λ = ρ 0 /ρ to transform C tot (λ, u, T) to C tot (ρ, u, T) and σ(λ, T) to σ(ρ, T). To convert the system to conservative quasi-linear form we use the Jacobian matrix To convert the system into non-conservative quasi-linear form we first write down Eq (14), then evaluate the partial derivatives of the primitive variables with respect to t and x and rearrange terms to form the matrix A(W). A(W), we can use the eigenvalues λ i and eigenvectors K i = [k i1 k i2 k i3 ] T of these matrices, together with the generalized Riemann invariants across rarefactions and contacts, and the Rankine-Hugoniot conditions across shocks, to derive a complete solution f(x, t,W L ,W R ) to the Riemann problem.

Once we have A(U) and
The solution will consist of four constant states W L , W ÃL , W ÃR , W R separated by three waves which propagate at speeds given by the eigenvalues λ 1 , λ 2 , λ 3 . The left (1) and right (3) waves can be either shocks or rarefactions; the center (2) wave is always a contact. The star left state W ÃL and star right state W ÃR surround the contact; finding their unknown components in terms of the initial conditions W L and W R is the first part of solving a Riemann problem.
To determine if a wave i is a shock, a contact, or a rarefaction, we check if its characteristics are converging, parallel, or diverging, respectively, on the left (l) and right (r) sides of the wave: Note that contacts and degenerate shocks can both have parallel characteristics. The difference is that velocity changes across shocks, but not across contacts. We can derive generalized Riemann invariants from the eigenvectors of the three waves: These ordinary differential equations allow us to form two invariants that hold across each wave i. These invariants hold only for rarefactions and contacts. A component k ij of an eigenvector K i will be zero where the corresponding variable does not change across wave i. Note that for eigenvectors derived from non-conservative forms, this may only be true for rarefactions, depending on which set of primitive variables is chosen.
The Rankine-Hugoniot conditions express the fact that mass, momentum and energy must be conserved across shocks. They are written where S i is the speed of the shock at characteristic i, and the fluxes and conserved variables are those on the left and right sides of the shock. In a coordinate system moving with the shock, the Rankine-Hugoniot relations reduce to F(U l ) = F(U r ), since the local shock speed is zero.

Riemann solver for linear elasticity
Evaluating the Jacobian A(U) = @F/@U for our linear elastic system, we obtain The eigenvalues of this system are where the speed of sound a ¼ ffiffiffiffiffiffiffi ffi Er 0 p =r. This corresponds to a system with a shock or rarefaction traveling left at speed u − a, a contact moving with the material at speed u, and a shock or rarefaction traveling right at speed u + a.
In non-conservative form, the matrix is where a 0 = a| ρ = ρ 0 . The eigenvalues of this system are the same as those derived from the conservative formulation, which is a useful check. The eigenvectors of the non-conservative form are Looking at the zero components of the eigenvectors, we can see that the temperature is constant across the left and right waves, and the density and velocity are constant across the contact, which means the solution consists of the four states separated by the three waves traveling at speeds λ 1 , λ 2 , and λ 3 . There are two ways of deriving expressions for ρ Ã and u Ã . The first is to use the generalized Riemann invariants; the second is to use the Rankine-Hugoniot conditions. For this linear elastic system, both methods give the same result for ρ Ã and u Ã , and our analysis of the characteristic speeds will show that the waves are always degenerate shocks, so we will only show the derivation of the Rankine-Hugoniot conditions here.
First, we transform the speeds u L , u Ã , and u R , to a coordinate system that moves with the shock i by usingû In a coordinate system moving with the shock, the Rankine-Hugoniot relations reduce to F(U l ) = F(U r ), since the local shock speed is zero. So for each wave, we have two coordinate transformation equations and three Rankine-Hugoniot equations, one each for mass, momentum, and energy. Eliminating the hatted velocities and shock speeds and solving for u Ã for waves 1 and 3, we obtain two expressions for u Ã Solving this for ρ Ã , we obtain Alternatively, using the same five equations per wave, but eliminating u Ã and solving for shock speeds gives which describes a left shock traveling at the speed of sound, and a right shock traveling at the speed of sound.
We can see that both waves are shocks because Note that contacts can also have parallel characteristics, but we know these are shocks because velocity changes across them, since u L 6 ¼ u Ã and u Ã 6 ¼ u R in general. We call these shocks "degenerate" because they are discontinuous but travel only at the speed of sound, whereas shocks in more complex systems can travel faster than sound. Now that we can calculate ρ Ã , u Ã , S 1 , and S 3 , this completes the solution of the Riemann problem for linear elasticity. Note that the temperature does not appear in any of these expressions, and since this solution procedure does not model diffusion across the contact, the left and right temperature states never mix. Sampling the solution requires simply finding the shock positions at the desired time t, then outputting the correct one of the four constant states W L , W ÃL , W ÃR , or W R that corresponds to the desired spatial coordinate x.

Riemann solver for nonlinear elasticity
Evaluating the Jacobian A(U) = @F/@U for our nonlinear elastic system, we obtain The eigenvalues of this system are where the speed of sound a ¼ . This corresponds to a system with a shock or rarefaction traveling left at speed u − a, a contact moving with the material at speed u, and a shock or rarefaction traveling right at speed u + a. Note that for nonlinear elasticity the speed of sound depends on T, which will complicate the solution process.
In non-conservative form, the matrix is The eigenvalues of this system are the same as those derived from the conservative formulation, which is a useful check. The eigenvectors are Looking at the zero components of the eigenvectors, it appears that the temperature is constant across the left and right waves, and the velocity is constant across the contact, which means that the solution would consist of the four states  This is the more general solution structure, since it allows for differing temperatures T ÃL and T ÃR inside shocks, and it conveniently incorporates the continuity of stress across the contact, so we will use this structure below.
In the linear elastic case, we had only two unknowns ρ Ã and u Ã . In the nonlinear elastic case we have four: ρ ÃL , ρ ÃR , u Ã and σ Ã . Our solution strategy will be to create two functions u Ã = f L (σ Ã ,W L ) = f R (σ Ã ,W R ) which we will equate and solve iteratively for σ Ã . There will be one form of u Ã = f(σ Ã ,W) for rarefactions, which we will derive from the generalized Riemann invariants, and another for shocks, which we will derive from the Rankine-Hugoniot conditions.
Star region velocity function for rarefactions. Writing down the generalized Riemann invariants using the eigenvectors from the first non-conservative form and integrating, we have R k 12 which yield the invariants Invariants I 12 , I 21 , I 22 , and I 32 are obvious by inspection of zeros in the eigenvectors; the generalized Riemann invariants just confirm them. Invariants I 11 and I 31 can be solved in terms of elliptic functions, but we leave them unevaluated for brevity. In our Riemann solver, we evaluate them by converting the integral to an initial value problem of an ordinary differential equation and integrating numerically. Applying invariants I 11 and I 31 across the left and right waves and solving for u Ã , we obtain Star region velocity function for shocks. Applying the same velocity-transforming Rankine-Hugoniot procedure as we did for linear elasticity, we get with a similar expression for u ÃR which merely substitutes R for L and reverses the sign on the right-hand term. Root finder for star region velocity and stress. We now have the ingredients to write an implicit equation u ÃL = u ÃR across the star region for the four possible solution cases: left shock/right shock, left shock/right rarefaction, left rarefaction/right shock, and left rarefaction/right rarefaction. We solve this implicit equation using a numerical root finder, whose result is a value of σ Ã which makes u Ã invariant across the contact. During the root finding process, we choose the shock or rarefaction u Ã function on each side depending on whether the characteristics converge or diverge, respectively. Analysis of the shape of the u Ã function shows that it has cusps at σ Ã = σ L and σ Ã = σ R which must not be stepped across during root finding, but the function is otherwise tractable, sloping upward for increasing σ Ã with no local maxima or minima.
Finding star region densities and shock speeds. The root finder gives us values of u Ã and σ Ã . To find ρ ÃL and ρ ÃR for rarefactions, we use the fact that the temperature is constant across rarefactions, as seen from the zeros of the eigenvectors. If we solve for temperature in the stress equation, and equate it across a left rarefaction, we obtain This can be solved explicitly for ρ ÃL as a cubic in σ Ã , but this solution is numerically problematic because it involves complex-valued intermediate calculations, though the result is always real. So instead, we solve it numerically, which is simplified by the fact that T is a smooth function of σ. The solution for a right rarefaction is the same, but with R substituted for L.
For shocks, we use the Rankine-Hugoniot equations again, but this time solving for ρ ÃL and ρ ÃR since we have the values of u Ã and σ Ã available. The result is for the left side, with ρ ÃR the same but with R substituted for L.
Once ρ Ã is known for a shock, we can use the Rankine-Hugoniot equations one final time, this time solving for the left shock speed S 1 and the right shock speed S 3 : Solutions inside rarefactions. Once we have values for ρ ÃL , ρ ÃR , u Ã , and σ Ã , the last step is to derive expressions for ρ and u inside the rarefactions, since they change smoothly over some distance instead of discontinuously as in shocks. We use the fact that all characteristics start at point (x, t) = (0, 0) in the x-t plane. This means the slope of any characteristic on the left side is x t ¼ u À a, and on the right side is x t ¼ u þ a. Using these equations and invariants I 11 and I 31 from Eq (45), we obtain for the left and right sides We solve these using a root finder to determine ρ for a given point x t . Then we substitute that ρ into these equations derived from the same invariants to get u at the same point With this, the Riemann solver for nonlinear elasticity is complete. As in the case of linear elasticity, for any given point x t , we first determine which characteristics it falls between, then sample the appropriate part of the solution. Note that this solver is somewhat numerically intensive, since it nests root finders and integrators inside of root finders. This is acceptable here, since our goal is merely to validate another numerical scheme with it. But this solver would likely be too inefficient to incorporate inside another scheme such as a Godunov method.

Derivation of a Sedov-Taylor solver
A Sedov-Taylor blast wave [39] [40] [41] is a blast wave formed by placing a large amount of energy in a very small area of a stationary field whose ambient density and temperature are small by comparison. Later in the results section we will present a test problem of this type as a supplement to the other Riemann test problems.
The solution of this type of blast wave is symmetric around the blast center, and consists of two very strong diverging shocks separated by two rarefactions. We can solve for the shock positions and post-shock primitive values by using the strong-shock assumption that the speed of sound a u = a L = a R in the unshocked material is much less than the shock speed S L = S R . The solution for the rarefactions is more involved, and is not presented here.

ð58Þ
Substituting these into the Rankine-Hugoniot conditions with initial velocity u R = 0, eliminating the velocity and stress variables, and taking the limit as a R ! 0, we obtain an expression for the right post-shock density ρ ÃR as a function of the right initial density ρ R and right shock speed S R : Similar manipulations give us the right post-shock velocity u RÃ and right post-shock stress σ RÃ By symmetry, the left post-shock values are the same as on the right, but with the sign of the velocity reversed.
To get the shock speed S R , we first find the radius of the blast wave R(t) by dimensional analysis. Our initial values are the density of the unshocked field ρ u = ρ L = ρ R and the total energy E c deposited at the blast center. The unshocked velocity u u = u L = u R is zero, and the unshocked stress σ u = σ L = σ R is assumed to be small compared to that of the blast. So the dimensional quantities available to us are We can then construct a quantity with the dimension of length The final expression for the radius of the blast wave as a function of time also requires a dimensionless factor η s % 1 to set the scale of the result And since we can substitute this shock speed into equations Eqs (59), (60) and (61) to obtain the postshock primitive values.

Results
Here we show the results of eleven numerical test cases. The first three tests are straightforward shock tube problems, and are used to show the general nature of RRM's solutions and to illustrate RRM's operation using spacetime diagrams. The fourth through tenth tests are of more specialized initial conditions which were chosen because they may cause difficulties for various other types of numerical methods. The eleventh test case shows RRM's handling of free boundary conditions. Where available, an analytic solution from an exact Riemann solver or a Sedov-Taylor solver is shown for comparison to the simulated solution. The derivation of these solvers is presented in a separate section. After the test cases, we analyze how the simulation error decreases as we increase the computational effort by decreasing the user-specified maximum tracer particle error metric Δ max . Fig 13 shows a shock-tube-like evolution of left rarefaction, center contact, and right shock for nonlinear elasticity. Normally the contact would be s-shaped, because RRM is not naturally adiabatic across contacts, since the replacement process mimics heat diffusion. However, in this test case we have specified in the initial conditions that the left and right cell are two different materials with identical properties, which causes flattening across contacts to produce two cells, which prevents heat diffusion and results in a sharp contact.

Nonlinear elastic shock tube
Note that there is a velocity and stress spike at the shock front which we do not see in the linear elastic results. This is because when we model shocks in a nonlinear elastic medium using RRM, the shock itself is a dynamic phenomenon. We do not explicitly calculate the shock speed and use it in the simulation as Godunov methods do. The tracer particles in the wavefronts spanning the shock front travel at the local speed of sound, the same as all the other tracer particles in the system. The spike is analogous to the finite thickness that real physical shocks possess, though its overshot profile differs from shocks in a real material, which are approximately sigmoidal as shown by Alsmeyer [42]. Fig 14 shows a spacetime diagram of the shock tube problem, zoomed in to show the solution structure. Each expanding wavefront is represented by a colored triangle. When wavefronts collide and merge, the left wavefront in the merger is drawn under the right wavefront that it merges with. Any number of wavefronts may be merged together, either all at once if multiple wavefronts collide simultaneously, or successively if wavefront collisions are spread out in time.
The inverted white triangles are cells, which start at their maximum width and are progressively encompassed by the wavefronts expanding out from their left and right edges. Usually only one cell is created from a single wavefront, but along the contact, we can see single wavefronts flattening into two cells to prevent heat diffusion across the contact.
Linear elastic double shock Fig 15 shows an asymmetric double shock for linear elasticity. We can see that RRM keeps perfectly sharp shocks, and this is the case even for very long simulations or cases where the shocks reflect off the boundaries (which are not shown here). To prevent shocks from being erroneously merged together, we disable wavefront merging in systems like linear elasticity where we know that the solutions are always shocks, except in the case of wavefronts with very small differences between them.
One interesting peculiarity shows up on the right side of the graph, where we can see that at t = 0.5, the Riemann solver shows the shock slightly ahead of the simulation. This is not simply an error in shock placement. Instead, it is a reflection of the fact that RRM does not update the entire field at every time step, so the value of t may be slightly different for every cell. The shocks are only perfectly placed at the exact moment a replacement occurs for the wavefront Fig 14. Nonlinear elastic shock tube as a spactime diagram. Shock-tube evolution from t = 0.0 to t = 0.005, with initial conditions the same as in Fig 13, but zoomed in on the center of the simulation domain. Each colored triangle is an expanding wavefront, and each white triangle is a cell. We can see the closely-spaced small triangles at left capturing the rarefaction, the line of small triangles down the middle following the contact, and the shock front using the smallest triangles at the right side.
https://doi.org/10.1371/journal.pone.0186345.g014 Fig 15. Linear elastic asymmetric double shock. Two linear elastic cells converging at different speeds to form an asymmetric double shock at t = 0.5. The time value is chosen to highlight (at the right shock) how the shock placement can be slightly behind the analytic value. Because the time steps in RRM only affect individual cells rather than the entire field, the shock placement will be exact only at the moments when the wavefront spanning the shock is replaced. This placement error does not accumulate over time, and can be reduced to any desired degree by reducing the user-defined error metric limit Δ max . Initial conditions are (ρ l , u l , spanning the shock. If desired, we can increase this replacement frequency by decreasing the user-defined error metric limit Δ max of the wavefronts. Or if we wish to take a synchronized snapshot of the entire field, we could simply force all pending events to occur at the snapshot time. This placement error does not accumulate over time, so it does not affect the overall course of the simulation. Figs 16 and 17 show spacetime diagrams of the same initial conditions, zoomed in to show the solution structure at two different levels. Note how the pattern of triangles is asymmetric, reflecting the asymmetric nature of the initial conditions. Note also the horizontally-striped "rainbow triangles" on the left and right. These illustrate an optimization called "wavefront spanning", where we avoid creating adjacent cells with the same primitive values by expanding existing wavefronts to cover them. The expansion process creates a new wavefront and merges the old one onto it, which creates a new stripe of color since wavefronts are colored by modulo indexing into a color table using the wavefront serial number. However, a single rainbow triangle may be considered to be one expanding wavefront, albeit one where the left and right tracer particles are changed over time. solver data closely. These initial conditions include a larger-than-normal value of C v to make the test harder by increasing the step size across the contact. Fig 19 shows a spacetime diagram of the same initial conditions, zoomed in to show the solution structure. Note how the cell sizes vary smoothly across both space and time in accord with the local demands of the simulation, without any sharp demarcations for submeshes or other stepwise adaptive techniques.

Stationary contact
Fig 20 shows a stationary contact in a nonlinear elastic material. These contacts are formed by finding two (ρ, T) pairs that result in identical stresses, and using those for the initial conditions. In Eulerian methods, contacts typically broaden over time due to numerical diffusion. In some Lagrangian methods such as smoothed particle hydrodynamics, contacts can change shape due to a numerical surface tension effect, depending on how those methods are implemented [43] [25] [44]. Contacts in RRM maintain their sharpness indefinitely, as long as we make the simulation adiabatic. Moving-grid methods can also maintain sharp contacts, at the cost of occasional remeshing operations when the grid becomes too deformed. RRM has the same property, but with continuous rather than occasional remeshing.   [45] or impact test in a nonlinear elastic material. This is a test of shock placement, symmetry, and the treatment of the trivial contact at x = 0.0. Even some recent and sophisticated Lagrangian techniques such as those described by Hopkins [25] and Frontiere et al. [44] can show a phenomenon known as wall-heating at the trivial contact, which causes spurious dips or peaks in the primitive variables. RRM shows accurate shock placement and no evidence of wall-heating, though it also displays a typical impulse-like overshoot at the shocks where the primitive variables transition from their pre-to post-shock values.  [38], but at higher separation velocity. RRM shows a flat solution across the trivial contact at x = 0.0, where some other methods may show dips in the primitive values at this location [38]. This test also stresses our nonlinear elastic Riemann solver, which required a more accurate numerical integrator for Eq 46 to handle these initial conditions. These initial conditions are likely beyond the range of physical validity of our constitutive equations, so we use them here merely to demonstrate numerical robustness.   in numerical methods for fluid dynamics [46], and the same can be true for nonlinear elasticity [38]. In this test case, we see that at the location of the sonic point at x % −0.47, RRM does not produce any spurious jumps in the primitive values.

Woodward-Colella-like blast wave
Fig 25 shows a blast wave in a nonlinear elastic material with a temperature ratio of 10,000:1 from left to right, which is similar to the pressure ratio across the left side of the Woodward-Colella blast wave test [47]. The solution is a tall and narrow right shock, the trailing side of which is a contact moving at the same speed. This tests a numerical method's ability to resolve parallel discontinuities of differing types, as well as general numerical robustness. RRM does well on this test, with the exception of a small overshoot at the shock front.  [41], but in a nonlinear elastic material instead of a gas. A large amount of energy E c is placed in a narrow center region of an otherwise stationary and uniform material. These initial conditions evolve into two diverging stong shocks, separated by two rarefactions. This is a test of numerical robustness, because a numerical method must be able to convert a very concentrated energy source into a symmetrically expanding blast wave.

Sedov-Taylor-like blast wave
Comparison of the shock placement and post-shock values to an analytic solution shows good agreement. An analytic solution is not currently available for the interior of this problem, but the qualitative features of the solution are similar to those of the analytic solution for the Euler equations [44], namely the convex-upward density curve, approximately linear diverging velocity, and stress which is flat in the center. Unlike for the Euler equations, these nonlinear elastic equations do not pull the density and stress very low in the center, due to the quadratic term in the stress equation Eq (7).
Nonlinear elastic double shock with free boundaries Fig 27 shows an asymmetric double shock in a nonlinear elastic rod with two free boundaries (the two ends of the rod). The rod initially consists of two converging cells. We can see that the length of the rod decreases initially as the two shocks propagate outwards from the center, then the rod rebounds out to its maximum stretched length at approximately t = 8.0, then it contracts and starts the cycle over at approximately t = 13.0. As the rod oscillates in and out, the shocks gradually lose their sharpness as they cross each other and reflect off the ends of the rod. Eventually the primitive values form smooth curves across the length of the rod as it continues to oscillate indefinitely.
The usual free boundary condition of an elastic substance is that σ = 0 at the boundaries. RRM does not need to enforce this condition explicitly. Instead, it happens naturally because at the ends of the rod, the stress momentum and stress energy are not balanced out by those of an adjacent cell, so any stress at the ends is quickly converted to kinetic energy.

Error analysis
If we define the point error between a simulation and a Riemann solver as e(x, t) = W RRM (x, t) − W Riemann (x, t), then we can define a useful measure of the error of an entire simulation run as  where e max is the maximum value that the space integral of the absolute value of the error takes for any time t during the simulation run. If we plot this quantity versus the maximum number of cells n used at any time during the simulation run, this will show us how the error decreases as the user decreases the tracer particle error metric limit Δ max , and will give us a measure of the computational requirements of the simulation. The result, shown in Fig  28, is that the error is approximately proportional to n −1.4 . And since the computational effort of RRM scales as OðNÞ, where N is the number of cells, the error in the simulation decreases slightly faster than linearly as we increase the computational effort. Since RRM uses constant-valued cells, we might expect the error to decrease only linearly as we add cells. However, since the cells are not evenly spaced, but instead are concentrated in areas of higher primitive variable slope, we obtain this slightly superlinear behavior. We have not yet performed an analysis of how alternative error metrics to Eq (11) might affect this behavior.
To see what RRM's error behavior looked like before its improvements for this paper, we redo the error analysis on the same test case, but with a diffusive (non-adiabatic) contact instead of a sharp contact. This gives the results shown in Fig 29. Note that e max , shown by the three dashed lines, shows approximately 10x less reduction at 400 cells than it did in the sharp contact case. In particular, the density component of the error e maxρ is much worse than before. This is because in RRM, diffusive contacts converge to a sigmoidal density curve, not a step, so RRM is unable to reduce the difference between its solution and that of the Riemann solver. And since all the error reduction must come from the non-contact parts of the domain, the system must work much harder, which is reflected in a sublinear behavior where the error is approximately proportional to n −0.65 .

Summary and conclusions
We have shown that RRM performs well on problems in linear and nonlinear elasticity, giving results whose shocks and rarefactions match the Riemann and Sedov-Taylor solvers, with a well-behaved error that decreases somewhat superlinearly with increased computational effort. We have also shown how to motivate and construct a simple constitutive equation for nonlinear elasticity, and have demonstrated in the derivation section how to create Riemann solvers for both types of elastic system.
The derivative-free, solver-free nature of RRM was chosen for robustness when simulating mathematically inconvenient systems. Now that RRM has been verified to work correctly for the Euler equations, linear elasticity, and nonlinear elasticity, future work on RRM will concentrate on systems for which no analytic solution is known.
In addition to application to new systems, future work is also possible in extending the generality and speed of RRM within existing systems. This work is described below.

Extension to 2D and 3D
The current implementation of RRM models cells as finite elements. As simulation progresses, these elements are chopped up and replaced with new cells. This approach could be directly extended to higher dimensions, but it would likely prove cumbersome to implement, since it Error versus computational effort, diffusive contact. The three components e maxρ , e maxu , and e maxσ of the maximum integral error e max versus the maximum number of cells, for the nonlinear elastic shock tube test case from t = 0.0 to t = 0.4, with a diffusive contact. The error in all three primitive values is approximately proportional to n −0.65 where n is the maximum number of cells used in the simulation run. This is significantly sublinear because without sharp contacts turned on, RRM is unable to reduce the difference between its solution and that of the Riemann solver around the contact. Note that the density component of the error e maxρ is the worst of the three, since the density is what changes across a contact. The value α is an arbitrary constant chosen to put the proportion line below the others for ease of comparison. Not shown on this graph, all three components of Δ max , corresponding to the user-specified error metric limit for ρ, u, and σ, are decreased from 10 −1 to 10 −5 from left to right. would require cells to increase in geometric complexity as they are chopped by other cells at any angle or position. Complex cells of this type would require more work to test for intersection, since we could no longer assume convexity, and such cells could become degenerate in a variety of ways.
To simplify the implementation, one possibility is "sampled RRM" (SRRM), which replaces the cells with finite-mass points, and replaces numerical integration with discrete summation over those points. This would take RRM's implementation closer to that of point-based Lagrangian methods like SPH. A kernel would also be needed in order to calculate the local speed of sound used by the tracer particles. SRRM would still differ from SPH and similar Lagrangian methods in that the velocities of the mass points would be unchanged from creation until chopping. The chopping out of wavefronts would also be performed in the same way as in RRM, albeit against mass points instead of cells.

Extension to high-performance computing architectures
The current implementation of RRM is a single-threaded event-based simulator. To run on HPC (high-performance computing) systems, which are typically made up of many thousands of computing nodes connected by a high-performance network, RRM will require parallelization.
Numerical methods based on vector and matrix operations are often amenable to intranode parallelization via vectorization libraries such as OpenMP (Open Multi-Processing) [48]. However, event-based simulators like RRM are harder to parallelize in this way, due to the data-dependent nature of the branches in the code. So parallelization of RRM will likely take place at the inter-core level using threads, and at the inter-node level using a library such as MPI (Message Passing Interface) [49], which is widely used in the HPC community.
Previous work in the field of parallel discrete event simulators (PDES), including the seminal works of Jefferson [50] and Fujimoto [51] as well as more recent work by Barnes et al. [52], suggests how RRM could be parallelized at scale. We would divide the domain into regions, each of which would have its own event queue and run in its own thread on a dedicated core. Event queues of neighboring regions would maintain loose time synchronization, requiring a tighter synchronization handshake only at the substitution of a wavefront that spans two or more regions. The PDES literature also suggests optimizations such as the use of speculation and rollback to permit looser coupling of event queues [50].
Additional optimizations are possible in the specific case of PDES applied to RRM. The simulation domain of a physical system is simply-connected, unlike systems such as electronic circuits which may be multiply connected. This simple connectivity should make it easier to split the domain at low-activity areas to reduce the frequency of event synchronization operations. And since RRM simulates a physical system with conserved quantities, occasional nonmonotonicity in event times near region boundaries may be acceptable as long as conservation is strictly maintained.