Qualitative properties of mathematical model for data flow

1.   Introduction

The increasing number and diversity of processing units in modern supercomputers presents a significant challenge to the understanding of how data is globally distributed in these systems. This issue is critically important in emerging computer architectures for which the cost of data movement is becoming a critical driver for system design, affecting hardware and software choices and even basic algorithms [12, 17, 18]. At the same time, these supercomputers are becoming increasingly difficult to model. Indeed, simulating the detailed state of a supercomputer requires an even larger computer, thereby making such approaches prohibitive for the largest systems.

With the help of new mathematical models, a more efficient utilization of available processing units might be possible. Such models are necessarily coarse-grained approximations, used to understand the gross behavior of the system in order to make design, control, and optimization more tractable. In this spirit, a coarse-grained mathematical model has been recently derived in [5] to describe the movement of processed data in an extreme-scale computer system. This model focuses on a idealized setting, in which a lattice of processors perform operations in parallel and communicate with nearest neighbors. Although very simple, the model allows for local variations in processor speed, due to differences in hardware performance or local problem complexity, as well as asynchronously operation, meaning that processors manipulate and share data with their neighbors as they are able, as opposed to a staged process in which processors communicate only after they have all completed a step of assigned work.

The model proposed in [5] takes the form of a PDE for the density of data at a particular processor location and stage in the global computation. It is derived formally from a microscopic ODE description. The notion of coarse-graining discrete and semi-discrete model of networks and related systems to derive PDE models is common mathematical strategy that is used in many applications, including vehicular [4, 11, 16] and pedestrian traffic flow [1, 6, 8, 15], supply chains [2, 3], bacterial movement [14] and machine learning applications [7].

In the current work, we establish some qualitative properties of the coarse-grained PDE model in [5]. We focus specifically on the phenomena of front propagation and present simulation results that highlight the theoretical findings. In the context of the current application, these fronts describe the movement of data through the different stages of a computational task. The reason for this analysis is two-fold. First, it was observed in [5] that, after a transient phase, many initial conditions eventually relax to a front whose shape depends on local variations in the processor speed. Second, these fronts can be used to characterize three important features for application: the first time at which some portion of the data reaches a state of completion, the rate at which the rest of the data reaches this state, and the time at which all data reaches this stage. Moreover, the relatively simple form of these fronts enables the investigation of optimization and control strategies based on some application-specific metrics. As a first step in this direction, we explore how variations in the processor speed affect front behavior.

2.   Mathematical description of the model

In this section, we briefly summary the microscopic ODE model from [5] and the macroscopic model that is obtained from it. The microscopic model describes a computer system as a lattice of processors, the dimension of which can be arbitrary, but finite. For application purposes, only one-, two-, and three-dimensional lattices are realistic, and here we focus the one-dimensional case for simplicity. It is assumed in [5] that each processor performs the same task, which is broken into discrete stages. The rate at which a given processor moves data from stage to stage depends not only on its intrinsic processing rate, but also on the availability of usable data in the processor and its neighbors.

Mathematically, the amount of data at time $t$ that sits in stage $k \in \{1,\dots,k_{\max}\}$ of processor $i \in \{1,\dots,i_{\max}\}$ is given by $q_{i,k}(t)$. In the absence of any throttling, the flow of data between stages is given by a (processor dependent) rate $a_{i}(t) \geq 0$. However the actual flow may be reduced due to a lack of available data in a given processor or its nearest neighbors. Throttling due to a lack of usable data in the nearest neighbors is characterized by the parameter $\eta: = k_{\max} / i_{\max} >0$. As $\eta$ increases the global effect of local slowdowns increases.

For each $i \in \{1,\dots,i_{\max}\}$ and $k \in \{1,\dots,k_{\max}\}$, the evolution of $q_{i,k}$ is modeled by the ordinary differential equation

where ${q}^0_{i,k}$ is the (known) initial amount of data present in each processor $i$ at each stage $k$ and $f^{\rm{in}}_{i}(t)$ is the (known) inflow of unprocessed (or raw) data at processor $i$ and time $t$. For $k \in \{1,\dots,k_{\max}\}$, $f_{i,k}$ is the flow of data in processor $i$ from stage $k$ to $k+1$. It is given by the product of the maximum processing rate $a_i$ and the composition of two throttling functions $v_1$ and $v_2$; that is,

where $Q_{i,k}(t)$ is the total amount of data in processor $i$ that by time $t$ has traversed the first $k-1$ stages of the computation:

and $Q_{i\pm1,k}-Q_{i,k}+q_{i,k}$ is the amount of usable data available from processor $i \pm1$. The throttling functions are given by

The function $v_1$ is a linear ramp, sometimes referred to as roof-line model [19]. The steepness of the ramp is determined by the parameter $q_\ast>0$, the value at which the processor is fully on. The function $v_2$ determines the amount of data ready to process based on what is available locally and from nearest neighbors. In particular, if not enough data is available from neighboring processors, then $f_{i,k} < a_i$.

¹ In [5], there is another parameter $\beta \in [0,1]$ in the definition of $v_2$ that measures the effective strength of coupling between processors. Here we assume $\beta = 1$ and do not consider it any further.¹

The macroscopic model in [5] is a continuum approximation for (1) that is derived in the limit of infinitely many processors ($i_{\max} \to \infty$) and stages ($k_{\max} \to \infty$), subject to that constraint that $\eta$ is a fixed, positive, finite constant. This continuum model is a PDE for a function $\rho$, where $\rho(t,x,z)$ is interpreted as the local density of data at time $t$, processor location $x$, and stage completion variable $z$. It takes the form

where $\mathbb{R}^+ = (0, \infty)$; $\mathbb{T}$ is the torus parameterized by $x \in [0,1)$; $F^{\rm{in}}$ is a scaled, continuous version of $f^{\rm{in}}$; $P$ and $\Phi$ are given by

and $w_1$ and $w_2$ are scaled version of $v_1$ and $v_2$:

The model (6) is a conservation law that describes the movement of data entering the system at $z = 0$ and exiting at $z = 1$. The flux $\Phi$ is the composition of the two scaled throttling functions: First $w_2$ determines the amount of data available based on the local density and the density of left and right neighbors. Once $w_2$ is specified, $w_1$ computes the effective processing speed assuming a linear ramp, where the parameter $\rho_* > 0$ is the minimum data density needed to operate at full capacity. The parameter $\eta$ controls the degree of throttling due to the lack of available data in neighboring processors.

At first glance, (6) appears to include a recursive definition of the flux function $\Phi$ due to equation (7a). However, after integrating (6a) over $(s,z) \in [0,t] \times [0,1]$ and applying (6c) and (6b), the result can substituted into (7a) to find the following formula for $P$:

Alternatively, since $\partial_z P(t,x,z) = - \rho(t,x,z)$, (6a) can be rewritten as a closed Hamilton-Jacobi equation for $P$:

The reformulation of (6a) in terms of (10) is the basis for the analysis and simulation performed in [5]. In the current work, we proceed by analyzing (6) directly in terms of $\rho$.

3.   Qualitative properties of the mathematical model

In this section, we investigate qualitative properties of (6). We begin by simplifying the expression for $\Phi$ in (7b) using the auxiliary function

where

(see Figure 1). In terms of $w_1$ and $W_2$, $\Phi$ takes the form

Next for simplicity, we artificially extend the domain in $z$ to $[0,\infty)$ and assume initial data with compact support. This assumption allows us to ignore the contribution of the outflow in (7a). Under this assumption, the definition of $P$ is modified as follows:

Hence, in the following we consider the model (6) on the extended domain

and introduce the auxiliary variable

The resulting system for $(\rho,\sigma) \colon \bf{D}^2 \to \mathbb{R}^2$ is then

The relation of the quantities $(\rho,\sigma)$ to the original model (6), modified by the previous Assumption, is as follows:

Definition 3.1 [Weak solution] Given initial data $\rho_0 \in L^\infty( \mathbb{T} \times \mathbb{R}^+ )$ and boundary data $(\rho_b,\sigma_b)\in L^\infty( \mathbb{R}^+ \times \mathbb{T})^2$, we call $(\rho,\sigma) \in [L^\infty({\bf D})]^2$ a weak solution of (17) if for all smooth compactly supported functions $\varphi:{\bf D}\to \mathbb{R}$ and $\psi:{\bf D}\to \mathbb{R}$,

We do not know whether a solution in the sense of Definition 3.1 exists for (17) with general initial and boundary data or whether such a solution is unique and depends continuously on the data.² However, for the particular initial and boundary data given in (28), we establish the existence of special solutions below.

²In [5], a continuous vanishing viscosity solution is established using known results from the mathematical literature on Hamilton-Jacobi equations. However, translating such a result to an $L^1$-theory for $\rho$ would require $P$ to be absolutely continuous. While possible, this additional regularity has not yet been determined.

The weak formulation (19) gives rise to a Rankine–Hugoniot jump condition [10]. Consider the surface

for a differentiable function $\zeta \colon \mathbb{R}^+_0 \times \mathbb{T} \to \mathbb{R}^+_0$ and assume that the functions $\rho$ and $\sigma$ satisfy (17) point-wise in the interior of ${\bf D} \backslash S$. Then standard arguments (see for example [13, Section 11.1.1]) can be used to establish that

where

3.1.   Special solutions for front propagation

In this section, we investigate solutions to (17) that take the form of fronts; see Figure 2. As discussed in the introduction, such solutions are important from the point of view of applications.

3.1.1.   Existence of Solutions

In the case of a single front, the initial and boundary conditions for $\rho$ take the form

where $H$ denotes the Heaviside function and $r>0$ is a positive constant. In addition, we enforce the consistency of $\sigma_b$ as prescribed in (18a):

Theorem 3.2. Let $\rho_*,\eta>0$ be positive constants and let $\alpha = \alpha(t,x) \in C^1( \mathbb{R}^+ \times \mathbb{T} )$ be a strictly positive function. Given initial and boundary conditions of the form (24) and (25), where

and $\zeta_0 \in C^1(\mathbb{T})$ is non-negative, assume that there is a $C^1$ solution $\zeta: \mathbb{R}^+_0\times \mathbb{T} \to \mathbb{R}^+$ that satisfies

Then there pair ($\rho$, $\sigma$) given by

is a weak solution for (17) in the sense of Definition 3.1 with initial and boundary conditions given by (24) and (25).

Before proving this result, some remarks are in order.

1. The condition in (26) guarantees that $\frac{r}{\rho_*} W_2( \frac{\sigma}{r} ) < 1$. In this case the flux function $\Phi$ (c.f. (13)) is given by the simplified formula

This formula and the Rankine–Hugoniot condition in (21) together form the key components of the existence proof.

2. While the condition $r<\rho_*$ is sufficient to obtain (29), it is not necessary. Indeed $W_2(\sigma/\rho) = 0$ for $\sigma$ sufficiently large; in particular, (29) holds at any $(x,t)$ such that $\partial_x \zeta(t,x)> \eta^{-1}$.

3. It is currently open as to whether equation (27) has a $C^1$ solution. However, since $\alpha, \rho_*$, and $W_2$ are non-negative, it is clear that any such solution will be non-decreasing. This fact is consistent with the notion that data is always processed toward completion.

4. A result similar to Theorem 3.2 can be obtained for initial data separated by a finite number of non-intersecting fronts. Consider for example the case of two fronts (see Figure 2-b):

with positive constants $r_1 \not = r_2$ and $r_1<\rho_*$ and $r_2<\rho_*$. If $0<\zeta_{1,0}(x)<\zeta_{2,0}(x)$ for all $x \in \mathbb{T}$, then there is a weak solution $(\rho,\sigma)$ in the sense of Definition 3.1 of the form

provided $\zeta_1(0,x) = \zeta_{1,0}(x)$, $\zeta_2(0,x) = \zeta_{2,0}(x)$, and

Proof of Theorem 3.2. With the $C^1$ assumption on $\zeta$, $(\rho,\sigma)$ given by equation (28) clearly belongs to $L^\infty({\bf{D}})^2$. By construction the given solution (28) fulfills boundary and initial condition pointwise and hence also in weak form. It therefore remains to verify that $(\rho,\sigma)$ is a weak solution in the interior of $\bf{D}.$

We first verify that $(\rho,\sigma)$ fulfills (19b). For $\psi$ smooth and compactly supported in the interior of $\bf{D}$, we need to show that:

By definition of the Heaviside function,

where the integral in the last line vanishes due to periodicity of the domain with respect to $x.$

We next verify that $(\rho,\sigma)$ fulfills (19a). Let $S = \{ (t,x,z)\subset \bf{D}: z = \zeta(t,x) \} \subset \bf{D}$ be the surface of the front described by $\zeta$, and let $\bf{S}_+: = \{ (t,x,z) \subset \bf{D}: z>\zeta(t,x) \}$ and $\bf{S}_-: = \{ (t,x,z) \subset \bf{D}: z<\zeta(t,x) \}.$ For $z<\zeta(t,x)$ the explicit form of $(\rho,\sigma)$ implies the following relation

Hence according to (29),

Let $\varphi$ be any smooth, compactly supported function in the interior of $\bf{D}.$ In $\bf{S}_{\pm}$, $\Phi(\rho,\sigma)$ and $\rho$ are constant, and in $\bf{S}_{+}$, they both vanish. Hence,

where $\mathrm{d} A(t,z)$ denotes the surface measure on the curve $\{ (t,z): z = \zeta(t,x) \}$. The final expression above vanishes due to the definition of $\zeta$ given in (27). This finishes the proof.

3.1.2.   Properties of fronts and their associated weak solutions

It is the equation for the front $\zeta$ in (27) that essentially determines the behavior of the special solutions (28). We assume again that $0 < r<\rho_*$. With the definition of $W_2$ in (11), this equation takes the form

In general, we do not know if there exists a solution to (38), due to the discontinuous flux. However, we may consider particular solutions related to special initial conditions $\zeta_0(x).$ Those special solutions are used as test cases for numerical simulations in the next section.

● If for all $(t,x) \in \mathbb{R}^+ \times \mathbb{T}$, $|\partial_x \zeta(t,x)| < \eta^{-1}$ and $\alpha(t,x) = \alpha$ is constant, then formally $\partial_x \zeta(t,x)$ fulfills the conservation law

with flux function $f(u) = C |u|$ and $C = \frac{\alpha\eta}{\rho_*}.$ It is known that there exists a weak (entropic) solution to (39) as long as $\partial_x \zeta_0(x)$ is in $L^\infty$, has locally bounded variation, and is continuous from the left [9, Proposition 3.1]. Moreover, if $\partial_x \zeta_0(x)$ is piecewise constant, this solution is given by piecewise constant states separated by a finite number of traveling discontinuities [9, Lemma 3.2]. We investigate such piece-wise solutions in Examples 2 and 3 in the numerical results of Section 4.

● If $\alpha$ is constant, and if $0 \leq \partial_x \zeta_0(x) < \eta^{-1}$ locally in $x$, then for sufficiently small $t$,

Similarly if $-\eta^{-1} < \partial_x \zeta_0(x) \leq 0$ locally, then for sufficiently small $t$,

We use these formulas to make comparisons with numerical results in Examples 2-4 in Section 4.

● If locally $|\partial_x \zeta_0(x)| \geq \eta^{-1}$, then $\zeta(t,x) = \zeta_0(x)$ for $t$ sufficiently small. In this case the model predicts that all processors are stalled. However, the condition on $\zeta_0$ cannot hold globally due to the periodicity assumption in $x$. We explore the behavior of the model for this scenario in Example 5 in Section 4.

4.   Numerical scheme and computational simulations

In this section, we construct a numerical discretization for the system. Results of numerical simulations using this discretization are presented to illustrate the behavior of special front-type solutions discussed in the previous section.

4.1.   Numerical discretization

Our strategy for approximating (17) is based on the relaxed system of equations

where $\epsilon > 0$ is a small parameter. Formally, (17a) and (17b) are obtained from (42) in the limit as $\epsilon \to 0$. This limit is considered later after a suitable discretization of (42).

Before discretizing (42), we investigate hyperbolicity in the simplified case that $\alpha(t,x) > 0$ is constant.³ In this case, (42) takes the non-conservative form

³ Hyperbolicity will not change if $\alpha$ is space or time–dependent, provided it is sufficiently smooth.

where $U = (\rho,\sigma)$ and

Under the assumption (26) and provided that $\rho>0$, the simplified form of $\Phi$ in (29) implies that

where $s = \sigma/\rho$.

The system (43) is hyperbolic if for any $\xi = (\xi_1, \xi_2) \in \mathbb{R}^2$ the matrix $B(U, \xi) = \xi_1 B^x(U) + \xi_2 B^z(U)$ is diagonalizable in the field of real numbers [13]. Because the function $W_2$ is not smooth, $B(U,\xi)$ is not differentiable when $s = 0$ or $s = \pm 1/\eta$. However, away from these points,

and the eigenvalues of $B(U,\xi)$ are

where $(\lambda_1^*,\lambda_2^*)$ are the roots of the characteristic polynomial

For $|s| > \eta^{-1}$, the system is strictly hyperbolic. For $0 < |s| < \eta^{-1}$, the polynomial $p$ in (48) has discriminant

which, for $\xi_2 \ne 0$, is negative (so that $\lambda_1^*$ and $\lambda_2^*$ have non-zero imaginary component) if and only if

For any $\epsilon >0$, the non-hyperbolic region $\{ (\xi_1,\xi_2) : d_\epsilon(\xi_1,\xi_2) < 0 \}$ is nontrivial; see Figure 3. However, for any fixed $(\xi_1,\xi_2)$, there exists an $\epsilon^*(\xi_1,\xi_2)$ small enough that $d_\epsilon(\xi_1,\xi_2) \geq 0$ for all $\epsilon < \epsilon^*(\xi_1,\xi_2)$. Hence as $\epsilon \to 0^+$, the non-hyperbolic region vanishes. The numerical method below is constructed by discretizing the relaxation system in (42) and then formally setting $\epsilon = 0$. This fact partially justifies the use of (42) even though it is not everywhere hyperbolic when $\epsilon > 0$.

To derive a numerical scheme for (42), we first discretize in time using a first-order method. We are eventually interested in the $\epsilon \to 0$ limit of the fully discretized scheme; hence the equation for $\rho$ in (42a) is treated explicitly, but the equation for $\sigma$ in (42b) is treated implicitly. Let $\rho^n(x,z) \simeq \rho(n \Delta t,x,z)$ and $\sigma^n(x,z) \simeq \sigma(n \Delta t,x,z)$ be the discrete approximations of $\rho$ and $\sigma$, respectively, at time $n \Delta t.$ We set

To discretize (51) in $x$ and $z$, we introduce a bounded computational domain $(x,z) \in [0,1)\times [0,1]$, which is divided into $N^x\times N^z$ uniform cells of size $\Delta x \times \Delta z$. The cell centers are denoted by $(x_i,z_j)$, for $i \in \{1,\dots,N^x\}$ and $j \in \{1,\dots,N^z\}$; and the cell edges are denoted by $x_{i+\frac{1}2} = (i + \frac{1}2) \Delta x$ and $z_{j+\frac{1}2} = (j + \frac{1}2) \Delta z$, respectively, for $i \in \{0,\dots,N^x\}$ and $j \in \{0,\dots,N^z\}$. The approximate cell averages of $\rho^n$ and $\sigma^n$ on the cell centered at $(x_i,z_j)$ are denoted by

and

respectively.

We use a Lax-Friedrichs approximation of $\partial_z \Phi$ in (51a) and in (51b) a centered discretization for $\partial_x \rho$ and a one-sided discretization of $\partial_z \sigma$. The resulting fully discrete scheme is

where

The monotonicity of the numerical flux $F_{i,j+\frac{1}2}^n$ is guaranteed provided that

This condition motivates a dynamic choice of $\Delta t$: At each time level $t^n$ we set $\Delta t = \frac{\Delta z}{a^n}$ for

The right-biased stencil in the discretization of $\partial_z \sigma$ in equation (54b) warrants some discussion. Indeed, given that the boundary condition for $\sigma$ is given at $z = 0$, it seems more natural to use a left-biased stencil. However, such an approach does not allow the boundary condition to be computed implicitly via the consistency relation in (18a). However, if $\rho$ has compact support, then (18c) holds. As long as the support of the solution does not reach the boundary of the computational domain at $z = 1$, this condition can be enforced there, in an implicit fashion. This approach is consistent with the fact that $\sigma$ at a given $z_*$ is determined entirely by $\rho$ at values of $z > z_*$. Moreover, numerical calculations confirm that enforcing the boundary condition in this way yields stable results.

The final scheme used in Section 4.2 is given by the formal limit of (54) in the limit $\epsilon \to 0$:

This scheme above must be accompanied by boundary and initial conditions. The initial values $\overline{R}_{ij}^0$ are obtained by integration of the initial data $\rho_0$:

The cell averages $\overline{S}_{ij}^0$ are obtained by applying a midpoint rule to the consistency relation (18a):

For $\overline{S}_{ij}^n$, zero boundary conditions are prescribed at $j = N^z$ consistent with the right-based stencil in the discretization of $\partial_z \sigma$:

For $\overline{R}_{ij}^n$, boundary conditions at $z = 0$ are computed from $\rho_{b}(t).$ For the special solutions discussed in the previous section $\rho_{b} = r$ and hence

Due to the central discretization of $\partial_z \rho$, boundary conditions for $\overline{R}^n_{ij}$ at $j = N^z$ need to be prescribed. Since in the simulation no data reaches this boundary we implement Neumann boundary conditions

4.2.   Simulation results

Example 1: validation with the microscopic model. We consider a test problem that has been used to to validate the macroscopic model (17) and the numerical method (58) by comparing with a simulation of the microscopic model introduced in [5]. For this problem, $\rho_* = 1.0$ and $\eta = 1.0$. The initial and boundary conditions are given by

and the processor speed is

Both models are simulated up to a final time $T_{\rm{fin}} = 0.5.$

The $(x,z)$ domain for the macroscopic model is discretized using $N^x = N^z = 800$ cells, and simulated with the algorithm in (58). In order to be consistent with the choice of $\eta$, the microscopic model uses $i_{\max} = 800$ processors and $k_{\max} = 800$ stages. The microscopic model is integrated in time using the explicit Euler method.

In Figure 4 we show the initial density and the final density profiles for both models at time $T_{\rm{fin}}$. We observe qualitative agreement between the microscopic and the macroscopic solutions. In both cases, the data around $x = 0.5$ is processed more slowly than the data in the rest of the domain, due to a slower processing rate $\alpha$ there.

In the remaining examples of this subsection, we explore the behavior of front-type solutions. Unless otherwise stated the following parameters are used in all simulations:

The value of $T_{\rm{fin}}$ ensures that data does not reach $z = 1$; in particular the assumption in (14) holds. These parameter choices, along with the initial and boundary conditions ensure that

in which case the flux $\Phi$ is given by the simplified formula in (29).

Example 2: front with constant profile. We consider the case of an initial constant front, namely

which provides a trivial example for a piece-wise constant solution $\partial_x \zeta$ to (39). Moreover, the evolution of this front is given explicitly in (40) and (41):

In this example $\eta$ is set to $0.5$ although the solution in (69) is independent of this choice.

Numerical results are shown in Figure 5. Since the processing rate $\alpha$ is constant the front moves with the same speed for each $x\in[0,1)$, as expected. Moreover, the analytical and the numerical solutions provide the same position of the front at final time.

Example 3: $\textsf{V}$-shaped front with small $\eta$ We consider an initial condition characterized by a $\textsf{V}$-shape front under the condition that $\eta | \partial_x \zeta | <1$; cf. (38). This condition models the physical situation in which processors are not initially at the same stage of completion. Such a situation may be realized if a local group of processors is, at some previous time, slower than the rest [5]. Specifically, we let

where $z_0 = \zeta_0(0.5) < 0.5$ is a fixed parameter.

Based on the theoretical findings in [9] for (39), we expect that the front at later times will be piecewise linear. Moreover, the formulas in (40) and (41) can be used to generate the local solution analytically. Indeed, under the assumption $\eta | \partial_x \zeta | <1$, a small calculation shows that the local, short time solution away from $x = 0.5$ is given by

where $|\partial_x \zeta_0| = (1-2z_0)$. However at $x = 0.5$, Theorem 3.2 does not apply, since the front is not in differentiable there; hence (40) and (41) are not valid. We conjecture that at this point, the front moves forward at full speed $\alpha / \rho_*$ and, as it catches up with neighboring points in the front, these points also move forward at full speed. This means that the global front solution is given by

In Figure 6 we show numerical results for this test problem using the algorithm in (58).. We use a grid with $N^x = N^z = 800$ computational cells. We set $z_0 = 0.2$ (so that $|\partial_x \zeta _0(x)| = 0.6)$ and let $\eta = 1/1.1$. It is clear from these plots that the analytic solution (72) and the numerical solution provide the same position of the front.

Example 4: Smooth front with small $\eta$. We repeat the test problem above, except now the initial front is given by a smooth function:

In order to fulfill the condition $\eta |\partial_x \zeta|<1$ we choose $\eta^{-1} = \max\limits_{x \in\mathbb{T}} \zeta_0(x) + 1 = {\pi}/2 +1$.

In Figure 7 we show numerical results for this test problem using the algorithm in (58). We use a grid with $N^x = N^z = 800$ computational cells and simulate the solution to a time horizon $T_{\rm{fin}} = 2$. As in the previous example, (72) can be used to evaluate the location of the front analytically. The results in Figure 7 demonstrate that the numerical and analytical solutions agree even when $\partial_x \zeta_0$ is not piece-wise constant.

Example 5: $\textsf{V}$-shaped front with large $\eta$. We consider again the initial condition (70) characterized by a $\textsf{V}$-shape front. The model parameters are the same as in Example 3, except that $\eta$ is now large enough to ensure that $\eta | \partial_x \zeta | >1$. In particular $|\partial_x \zeta | = 0.6$, but $\eta = 10$. In this case, we expect stalling for those processors $x\neq 0.5$; see (38). Again at $x = 0.5$ the solution is not described by Theorem 3.2, since there the front is not $C^1$. We expect that the front at this point will proceed at the maximum speed $\alpha/\rho_*$ and that the stalled points away front $x = 0.5$ will also begin to move at this rate after they have been overtaken. In other words, we conjecture that the global formula for the solution front is

In Figure 8 we show numerical results for a grid with $N^x = N^z = 800$ computational cells and a time horizon $T_{\rm{fin}} = 2$. These results agree with the analytical solution in (74).

4.3.   Variations in the processor rate.

The discrete model (1) can be used to predict how variations in the processor rate affect computer performance, and to understand whether these rates can be controlled to achieve a prescribed objective, such as faster throughput. For the continuum model in (17), such a control is specified via the function $\alpha$. In this section, we briefly explore the effects of adapting $\alpha$ with the intent of a more extensive study in later work.

In order to normalize the control action, we assume that for all $t$,

Since $\alpha$ is non-negative, so too is $\bar{\alpha}$; we assume further that $\bar{\alpha}>0$. Then several choices for $\alpha$ can be envisioned. The simplest choice is to set $\alpha$ equal to a constant. In this case, the normalization in (75) implies that $\alpha = \bar{\alpha}$. For a front-type solution, this simple choice will be compared with a policy that depends on the initial data

Since the slope $\partial_x \zeta(t,x)$ controls the throttling due to neighboring processors, the condition $\partial_x \zeta(t,x) = 0$ is desirable. We therefore propose a choice for $\alpha$ that gives priority to processors that are trailing in the initial configuration. Specifically, we consider

where $\zeta_{\rm{max}}>\max_{x \in \mathbb{T}} \zeta_0(x)$ is a fixed constant and $C_\alpha$, which depends on $\zeta_0$, is chosen so that (75) holds. This choice of $\alpha$ will increase the speed of trailing processors

In order to compare the temporal and spatial evolution of $\rho$ based on the policy (77) versus the constant $\alpha = \bar{\alpha}$ policy, we compute some quantities of interest

The quantity $\omega_1$ measures the cumulative outflow at stage of completion $z$ up to time $t$; the quantity $\omega_2$ measures the current outflow at stage of completion $z$; and the quantity $\omega_2$ measures the cumulative part density at stage of completion $z$. The first two quantities are indicators of the processed data whereas the last indicator shows the load of the processors, which may be related to the consumption of energy during computation. In each case, larger values of $\omega_i$ are preferable.

For our numerical tests, we set $T_{\rm{fin}} = 1$ and $\bar{\alpha} = 0.5$. The initial front $\zeta_0(x)$ is given as in (70) with $z_0 = 0.1$, in which case $C_\alpha \approx 0.5/0.575$. The quantities of interest are evaluated for $z = 0.5$ and $z = 0.75$. Figure 9 shows the initial density profile and the profiles at $T_{\rm{fin}} = 1$ for each of the two policies. For the policy based on equation (77), the control has been been active long enough to push the processors that were initially trailing ahead of the others. This means the control has been on for too long and should have been modified or turned off at an some earlier time. While the results of the two policies appear similar, the quantities of interest $w_i(t)$, $i = 1,2,3$ displayed in Figure 10 demonstrate the differences. In particular, with respect to $w_1$, the policy in (77) outperforms the constant $\alpha$ policy.

5.   Summary

In this paper, we have analyzed a recently derived mathematical model for the evolution of processed data in large-scale, asynchronous computers [5]. After the introduction of an auxiliary variable, the model is expressed as a system of partial differential equations. It is possible to prove existence of discontinuous solutions to this system for a particular class of initial and boundary conditions. These solution describes the flow of information as a series of propagating fronts. A numerical scheme for the reformulated system has also been designed based on a relaxation approximation. Numerical simulations based on this scheme demonstrate qualitative agreement with theoretical findings.

We have also briefly explored the effects of local processor speed on quantities of interest predicted by the model. In future work, we intend to investigate more extensively control mechanisms for optimizing important objectives related to performance of the large-scale computers. Further it is planned to validate the model based on experimental data.

Acknowledgments

MH and GV thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for the financial support through 20021702/ GRK2326, 333849990/IRTG-2379, HE5386/14, 15, 18-1, 19-1, 22-1 and under Germany's Excellence Strategy EXC-2023 Internet of Production 390621612. The funding through HIDSS-004 is acknowledged. This work has been supported also by the US National Science Foundation, RNMS (KI-Net) grant 11-07444, and the U.S. Department of Energy, Office of Advanced Scientific Computing Research. The work of CDH was performed at the Oak Ridge National Laboratory, which is managed by UT-Battelle, LLC under Contract No. De-AC05-00OR22725.