Model for the evolution of the time profile in optimistic parallel discrete event simulations

We investigate synchronisation aspects of an optimistic algorithm for parallel discrete event simulations (PDES). We present a model for the time evolution in optimistic PDES. This model evaluates the local virtual time profile of the processing elements. We argue that the evolution of the time profile is reminiscent of the surface profile in the directed percolation problem and in unrestricted surface growth. We present results of the simulation of the model and emphasise predictive features of our approach.


Introduction
The Parallel Discrete Event Simulation (PDES) approach [1] is thought to be very promising for large scale simulations. The main attractive feature of the approach is that it is scalable by construction: i) there are limitations on the problem size [2], and ii) it is proven for some simple cases it is deadlock free [3]. Practical implementations for using a PDES approach in computational physics is limited, primarily due to the folklore that de-synchronisation between processes grows with the system size. We have to mention that this fact is proved for the conservative algorithm of PDES [3], and not adequately investigated for the optimistic algorithm. At the same time, the most powerful implementation of the PDES approach is the optimistic algorithm. It is widely used for simulations in areas outside of computational physics [4]. In this paper we investigate the problem of synchronisation in optimistic PDES using a simple model of the optimistic algorithm. This model has some similarity with the Korniss et al model [3] for the conservative algorithm, and captures essential features of the optimistic algorithm. The dynamical behaviour of the optimistic PDES model is quite different from the conservative PDES mode. The two models belong to different universality classes. The Korniss et al model belongs to the universality class of Kardar-Parisi-Zhang equation [KPZ], while our model demonstrates features of the universality class of a roughening transition [5].
In the second section we formulate our model. In the third section we present results and analysis. In the final section we discuss how predictive properties of our model can be mapped onto real simulations in computational physics.

Model
PDES considers the evolution in time of the ensemble of N objects. One can think on the object as an Ising spin in Glauber dynamics [6] or as a particle that forms clusters [7]. Each object i from (1, 2, ...N ) may change its state A i at some moment of time, this change is named an event.
In the examples, the event is the spin flip or the particle sticking to the cluster. The event happens at some definite moment of time, so the time associated with the event is a discrete variable. The idea of PDES is to replace the interaction of objects, which is usually realised in simulations through the access to the data held in the hardware memory, with the mechanism of message sending [1,2]. Object i while changing its state A i sends a message S(A i ; τ i ) with the time stamp τ i , which is the value of the local virtual time (LVT i ) of the object i. This message is received by all objects which may depend on the state A i . The ensemble of the LVT i forms the LVT profile over the objects.
Let us map the ensemble of objects onto the ensemble of processing elements (PEs). A PE may be the computing node, or CPU, or core, or thread, depending on the appropriate hardware/software combination available for the simulation [8]. Accordingly, the LVT is now associated with the PEs. Therefore, the problem of PE synchronisation is connected with the problem of keeping causality in the object evolution. Possible algorithms can be classified into three groups [8] and named as conservative, optimistic, and FaS. We do not discuss here details of the classification and of the associated algorithms, and refer the reader to a recent paper with a short review on the subject [9]. We concentrate here rather on the optimistic algorithm.
The optimistic algorithm name comes from the idea to give each PE the possibility to proceed, with the assumption that causality is not violated. Each PE i holds messages ordered by the time stamp in the list LR i of received messages. In the process of performing the simulation, the PE uses information written in the message, which is the state A j of the PE j changed at LVT=τ j . The process of simulation and the process of receiving messages are performed in parallel. Each PE i holds also a list of sent messages, LS i , which is also ordered with the time stamps.
If it happens that incoming message A j contains the stamp τ j with the time smaller than the time of the message on the top of the LS i , and if the state A i of that message has been simulated using the old state of processing element j, it means that causality has been violated. The problem of causality violation is solved using the mechanism of a rollback in virtual times. In this case, PE i should return back to the state A i contained in the message on the top of LS i . If it happens that the time stamp of this message is still larger than the time of the message on the top of LR i , than one more rollback is required to be performed, and it is necessary to inform other PEs that this message is not correct. This is done through the mechanism of sending anti-messages. PE i takes message A i from the top of LS i and sends anti-messageÃ i with the same information as in A i . PEs receiving anti-messageÃ i annihilate it with the corresponding A i . It is clear that the rollback process may create an avalanche of rollbacks of some length.
We propose the following model for the evolution of the time horizon in optimistic PDES algorithms. We abbreviate in the following our model as OA, optimistic algorithm model. For simplicity, let us consider that all PEs may be ordered in a one-dimensional chain, and only nearest neighbours depend on the information contained in the messages. These assumptions were used previously in a model of the conservative algorithm [3] and of the FaS algorithm [8]. The first step of the model consists of the forward evaluation of LVT i . We start with the flat profile of local virtual times of L processing elements, τ i = 0 with i = 1, 2, ..., N . We advance LVT i of each PE i with a random time step taken from the exponential distribution with mean F , this is the forward step of the model. Then we have to relax the LVT time profile under the rollback process. We assume that the length of the avalanche follows an exponential distribution with mean B, and each (randomly chosen) PE i may rollback B times to the right or left value of LVT, τ i−1 or τ i+1 .  Simulation shows that both the average LVT τ av (t) and the GVT τ min (t) grow linearly with time, and the slope of the curves define the corresponding speeds of the time profile, u av and u min . We found it useful to introduce a variable q = 1/(1 + B). Figures 1 shows variation of these quantities as a function of q for the parameter value F =1. A fit to the data with the form u = u 0 (q − q c ) ν gives u 0 ≈ 1, q c ≈ .115, and ν ≈ 1.73. The value of the exponent ν is very close to the one of directed percolation [10]. It is interesting to compare our model to the model of unrestricted surface growth (USD) [5]; they have some similarities and some differences. The similarity is in the existence of two main steps in both models: adsorption of atoms in the USD model and the time increment in the OA model, and the desorption of atoms in the USD model and rollbacks in the OA model. The differences are in the probabilities, which is q for adsorption and (1 − q)/2 for desorption in the USD model, while the probabilities in the OA model are independent parameters. The value of q c ≈ 0.233 in the USD model is about two times larger than in the OA model. There is the possibility to compare the time evolution in our OA model with the conservative algorithm model discussed in [3]. We can calculate the number of local minima N (t) of the time profile. This is the number of PEs which may successfully evaluate in the conservative algorithm, and number of PEs which will not rollback on the next step in the OA. Figure 4 shows the density of local minima N (t)/N of the time profile for a number of values of the parameter q, as well as for the conservative algorithm. We notice that formally this picture demonstrates that the optimistic algorithm is more efficient as the number of active sites is larger than in the conservative algorithm. At the same time, the speed of the time profile evolution is growing from the top of figure 4 from the bottom (starting from q = 0.185 and progressing toward q = 0.833). The speed of OA model is equal to the speed of the optimistic algorithm model at the value of q ≈ 0.6.

Discussion
We proposed a model for the evolution of the time horizon in optimistic PDES simulations.
Simulations of the OA model demonstrate that some properties of OA is reminiscent of those of unrestricted surface growth [5], which belongs to the universality class of directed percolation [11]. This observation may be used for the prediction of the behaviour of real optimistic PDES simulations. Predictions may be performed by matching parameters of the OA model to those of the PDES simulations. For example, optimistic Time Warp simulations show indications of self-organised criticality [12]. We plan to perform case studies of optimistic PDES, and analyse possible mappings of PDES onto the OA model. Probably, some extensions of the present OA model should also be investigated.

Acknowledgments
This work by LZ and LNS was supported by the Russian Science Foundation (project nos. 14-21-00158). MAN was supported in part by a grant from the US National Science Foundation (DMR-1151387).