Parallel discrete-event simulation (PDES) is a large-scale simulation technique with various applications in engineering, computer science, economics, transportation, the military, etc. For example, it is used to simulate modern computer networks [1], air traffic [2], disease spread [3], the evolution of biological systems [4], catalytic kinetics [5], and metal welding [6]. PDES also provides an environment for war-gaming exercises [7] and intelligent combat communications simulations [8].

The discrete-event simulation (DES) model assumes that the system being simulated changes state only at discrete moments in the simulated time. These changes are called events. Events are not synchronized by global time but occur at irregular intervals [9]. Each event contains a timestamp and usually represents some change in the state of the system being modeled. The simulation process can be described as a loop in which the event with the smallest timestamp is repeatedly removed from the event list and processed. Handling an event means changing some state variables and scheduling new events for the future.

Serial DES typically consumes a large amount of time/memory resources, but due to its design, it can also be easily parallelized. In parallel DES, the system being modeled is split into subsystems, and each subsystem, in turn, is modeled by a separate process called a logical process (LP). LPs are performed by actual processing elements such as threads, cores, or nodes. To preserve the causality of event processing, PDES uses special synchronization algorithms: the conservative algorithm [10, 11], optimistic algorithm [12], and freeze-and-shift algorithm [13].

A key concept of PDES used in synchronization algorithms is the concept of virtual time [12]. Each LP has a local variable denoting its local virtual time (LVT). The profile of local virtual times for parallel processes grows during the simulation. It is known that the evolution of the LVT profile in the conservative algorithm belongs [14] to the Kardar–Parizi–Zhang (KPZ) universality class [15] and that the evolution of the LVT profile in the optimistic algorithm belongs [16] to the directed percolation (PD) universality class [17, 18].

In this paper, we use the model of LVT evolution in optimistically synchronized PDES [16] for studying the properties of the personal communication service (PCS) model [19]. For this purpose, we associate the parameters of a theoretical LVT evolution model [20] with those of a realistic PCS simulation using an optimistic PDES engine called Rensselaer’s Optimistic Simulation System (ROSS) [21]. We find that the critical behavior of both models near the roughening transition is quite similar, suggesting that they belong to the same universality class.

The development of parallel discrete-event simulation began in the 1960s and 1970s of the last century, when the first commercial multiprocessor computers became available for researchers and engineers. The development of the method proceeded in two directions. First, the PDES method is of interest to scientists because it allows solving problems that cannot be solved using sequential algorithms. Second, PDES is interesting in its own right because it represents a particular problem area in parallel computing in general.

It seems that the first application of the PDES method in physics was carried out to simulate spin dynamics [22]. An important property of the parallel modeling method is that it is event-driven: instead of updating sites at regular intervals Δt, an event-driven simulation stores a history of updates, and the sequence of events (spin-flip attempts) does not change in favor of parallelization. In the following sections, we will refer to the method proposed in [22] as the basic conservative PDES scheme. This scheme is ideal for parallel implementation of the Markov Chain Monte Carlo algorithms [23] and any stochastic cellular automata with local dynamics, where the discrete events are Poisson arrivals.

In real systems, events have temporal order and causality relations in real time. Lamport was the first to show that it was possible to determine the order of events in distributed systems using a system of special artificial clocks. The clock is used to mark events with a unique value [24] in each process. He also formulated the so-called “Lamport clock conditions,” which are crucial for preserving causality when computing events in parallel. Jefferson later formalized the concepts of local and global virtual times [12] and introduced the famous Time Warp synchronization algorithm.

Local virtual time refers to the current time of logical processes, which is usually equal to the arrival time of the last processed events. Global virtual time (GVT) is the smallest LVT in the system. GVT is used to measure the progress of computation. The LVTs of all logical processes form the so-called LVT profile, which grows as the modeling process progresses. The instantaneous LVT profile is shown as an example in Figure 1. The concept of virtual time underlies all PDES synchronization protocols and, as will be discussed later in this article, is a useful tool for analyzing the efficiency and scalability of PDES algorithms.

It is quite interesting that the analogy of the particular PDES algorithm, the conservative algorithm, with the KPZ model [14] can be used to classify possible PDES algorithms [13]. The reasoning is as follows: the KPZ equation is a partial differential equation. The solution to the partial differential equation is defined by the choice of specific boundary conditions. [13] expressed the idea that the type of synchronization of LPs with each other is associated with the boundary conditions of the PDES algorithms. Thus, the continuous boundary conditions correspond to the conservative scheme, the free boundary conditions correspond to the optimistic class of algorithms, and the fixed type of boundary conditions should correspond to the new class of PDES synchronization algorithms. This class is called the “freeze-and-shift” algorithm [13].

In the following subsections, we briefly describe two LVT profile simulation models that reflect the general features of conservative and optimistic algorithms and have similarities with the phase transitions in the KPZ model [14] and directed percolation model [16].

A two-dimensional mobile network, also known as a personal communication service (PCS) network, is a wireless communication network with distributed communication ports, each with its own coverage area or cell zone, set of communication channels, and users [19, 29, 30]. Users send and receive phone calls using these communication channels. When a user makes a phone call, the system tries to allocate a channel for the connection. If all channels are busy, the phone call terminates. The network is represented by square cells, assuming that each cell is assigned a fixed number of radio channels. The goal of PCS modeling is to minimize the probability of dropped calls.

The PCS model is implemented in the ROSS simulator [21]. ROSS serves as a simulation engine that provides a simulation loop for generating, dispatching, and processing events by parallel processes. It is implemented in the C++ programming language and uses the message-passing interface (MPI) for parallel operation. Parallel processes are optimistically synchronized.

Let us describe the system in the language of logical processes and events. In this model, ports or cell zones are LP objects. LPs are located on a square lattice (Figure 2). The behavior of a phone call is modeled by different types of events. The system implements four different types of events: 1) NextCall—receipt of a call to a cell; 2) CompletionCall—completed cell call; 3) MoveOut—calling exit from the current cell (remote event); and 4) Moveln—indicates the arrival of a switching call in the cell.

In PCS, the time between events is distributed exponentially. The following parameters set the distribution means:

• MOVE_CALL_MEAN—average time between MoveOut events;

Thus, a call can be scheduled to be sent to the current cell or to one of the neighboring cells. In the ROSS implementation, such events passing between LPs are called remote events. The proportion of deleted events is available in the output statistics. Let us denote this fraction as p. In other words, LPs are located on a square lattice with periodic boundary conditions, and the interaction between nearby LPs occurs with probability p. In addition to the number of remote events p, ROSS statistics provides data such as the total number of events processed, the number of events rolled back, and the rate of events. LVTs are available during the simulation.

Let us construct a model of the evolution of the LVT in the described PCS model, the LVT-PCS model. An important difference between the LVT–PCS model and the general optimistic LVT model [25] is that the parameters p and q are now interdependent and cannot be chosen independently. However, we expect that in the vicinity of the roughening transition, the behavior of the LVT–PCS model will be essentially similar to that of the general LVT model and belong to the DP universality class.

In the LVT–PCS profile evolution model, we, as usual, use the assumption that events are Poisson arrivals, so the time between them is exponentially distributed. Simulation parameters are

• N—total number of LPs,

The values of q and p are set equal to the values of q and p, obtained from the PCS simulation on ROSS.

We first create the LP topology. Since the LPs in ROSS are implemented on a two-dimensional square lattice, in LVT–PCS, we also place the LPs on a two-dimensional square lattice (Figure 2). Each LP can exchange messages with its four nearest neighbors.

The simulation starts with a flat profile τ _i = 0 and i = 0. N − 1, where τ _i is the local virtual time of the ith LP. Next, at each modeling step t, the free growth of the profile is first simulated, and then, rollbacks are simulated.

Modeling free growth means increasing the LVT of each LP by an exponentially distributed random variable η _i with the unit mean: τ i t + 1 = τ i t + η i , i = 1 , … , N .

At this optimistic stage, we model the profile growth, assuming that no violation of the causality has occurred. After this step, we simulate rollbacks. The total number of LPs that need to rollback their LVTs is regulated by the growth rate parameter q (see Eq. 2). Let us denote the average number of rolled back LPs as bN. Assuming that all N LPs proceed at the first step and bN LPs rollback on average, from Eq. 2, we have q = 1 − b N N + b N = 1 1 + b . (3)

From Eq. 3, we obtain b = 1 q − 1 .

Therefore, to simulate rollbacks, we randomly select, on average, bN LPs, which local time will rollback. If the LVT of LP is greater than that of one of its neighbors, chosen with probability p, then the LVT is reduced to the neighbor’s time.

After each time step t, we calculate the average profile velocity v(t) (1) and the average profile width w ²(t) (3). In computer science parlance, the average speed of a profile v is related to the efficiency of the simulation, which is the average load of the parallel processes. In addition, the average profile width w ² is associated with desynchronization between processing elements. The larger the profile width, the more desynchronized the simulation.

The final results are calculated as the average of R independent runs of the random process: τ t = ∑ R τ t R , v t = ∑ R v t R , w 2 t = ∑ R w 2 t R .

The main result of our simulation is that the behavior of the LVT profiles in both models is qualitatively similar near the roughening transition q _c . Near the transition, there are more remote events, and the LVT growth speed becomes lower, as does the event rate (i.e., efficiency), and the average profile width stops growing [31, 20] (i.e., synchronization becomes better), and at the transition value q _c , there is a deadlock in the PCS model.

A comparison of the average profile velocity in the LVT–PCS and PCS models and the event rate in the PCS model is presented in Figures 4, 5. The average velocity in the LVT model corresponds to event utilization in the PCS models (i.e., event rate).

The dependence of the event rate ER on the value of the parameter q is shown in Figure 4. It can be seen that the symbols calculated for system sizes N = 64, 144, and 256 approach each other at a low event rate. We approximate the event rate data using the expression E R = E R 0 ( q − q c ) ν + const and estimated values q _c = 0.064 and ν = 1.345(9) in the case of the largest system size N = 256. Fitting smaller values of N is unreliable as it is well-known that for small system values, scaling corrections can be large [32]. The value ν ≈ 1.3 [17] apparently corresponds to the DP 2 + 1 universality class. The value q _c is a deadlock critical value, i.e., deadlocks occur at point q = q _c , and the profile speed becomes zero. The conditions mentioned at the end of Section 2.2 are satisfied—there is a roughening transition at q = q _c , the event rate is non-negative, the rules are local, and there is no quenched randomness. To summarize, we can conclude that the behavior of the PCS model near the deadlock value belongs to the DP 2 + 1 universality class.

We use the p and q values estimated in the PCS model to simulate the LVT–PVS model. Figure 5 shows the resulting dependence of the average profile velocity in the LVT–PCS model on the parameter q. There are three ranges of q with different velocity dependences. First, for the case of N = 256, this is zero velocity below the value q _c = 0.185(1), which is the roughening critical point [18]. Second, we use the profile approximation v = v 0 ( q − q c ) ν + const of the data in the range q _c < q < 0.3. We estimate the value of the exponent ν = 1.3(1), which is close to the value ν _DP of the directed percolation model in the 2 + 1 dimension, ν _DP ≈ 1.3 [18, 17]. Third, it is a random deposition mode, with a value close to 0.8, where the profile grows randomly on each LP and the average speed reaches a value of unity, the maximum possible profile speed.

Finite-size effects are visible in Figure 5 as a shift in the position q* of velocity that vanishes with the size of the system. General theory (for example, see [32]) predicts that the value of the parameter q* depends on the size of the system and scales with some exponent. We cannot correctly estimate q*(N) for systems of sizes 64 and 144 and provide only a qualitative explanation of the shift. The LVT–PCS model captures important properties of the PCS model that were not known prior to our analysis, exploiting the profile’s similarity to those of the statistical physics of surface growth.

Figure 6 shows the growth of the squared profile width in the LVT–PCS model for two system sizes N = 64 and 256 and for several values of parameter pairs (q, p). For values of q below the corresponding value of q*(N), the width does not increase after some time of transition to the absorbing state. At large values of q, saturation occurs as a function of the values of (q, p), and before this, the system has a value of the exponent β close to 1/2. We summarize estimates of the effective exponent β in Figure 7, which shows a smooth change in 2β, as q increases. This is another illustration of the cross-over from the DP fixed point to the random deposition fixed point.

In this study, we applied the previously proposed optimistic LVT evolution model [16, 20] to analyze the simplest mobile network model, the portable communication service (PCS) model. We propose a way to relate the parameters of the PCS model and our LVT model. We call the resulting model the LVT–PCS model, emphasizing the fact that it uses parameter values specific to the PCS model. This is a way to model the specific behavior of a PDES model using the LVT growth model.

We first simulate the PCS model using the ROSS simulator [21] and extract the values of q and p. Next, we use q and p as input parameters when simulating the LVT–PCS model. In the LVT–PCS model, the parameter q is related to the number of rollbacks in the optimistic algorithm and the parameter p is the probability of time synchronization with one of the cell neighbors. By matching the parameters, our model was used to analyze deadlocks and synchronizations in the PCS model.

We found a roughening transition at q _c in the LVT–PCS model, which apparently belongs to the directed percolation universality class of [17, 26] in dimension 2 + 1. This corresponds to a deadlock in the PCS model. For values of q > q _c , the velocity profile behaves as a power of q with an exponent close to 1.3, which is an estimate of the correlation length exponent in the DP 2 + 1 universality class. Thus, we know how the event rate behaves in the PCS model near the deadlock. An increase in the width of the LVT profile w was observed over time, with an exponent also corresponding to the DP 2 + 1 universality class.

There is a second extremum with a value of q close to unity, which corresponds to the universality class of random ballistic deposition. Again, the behavior of the PCS and LVT–PVS models is similar at this limit.

Between these two extremes, one should not expect universal behavior since it represents a transition from one fixed point of the DP universality class at q = q _c to a second fixed point of the ballistic deposition universality class at q = 1. The crossover region is known to reflect specific details of the models, and it is not surprising that slightly different behavior is observed in the intermediate q region at approximately 0.5.

We hope that the proposed approach can be generalized and applied to other optimistic PDES models. Future research may consider extending this approach to other models and assessing its effectiveness. Additionally, it would be useful to study the practical implications of this approach, such as its potential to improve system performance, reduce resource usage, and improve scalability. These results may have important implications for the computational community and ultimately lead to more efficient and effective modeling methods.

The work was carried out at the Institute of Solid State Physics of the Russian Academy of Sciences within the framework of the state assignment of the Russian Ministry of Education and Science.