1 Introduction

Pipelines play an important role in the infrastructure that provides human society with crucial energy resources [1]. Tens of thousands of pipes in our cities carry and distribute the essential energy resources for daily life, such water and natural gas, from their origins to the homes of individual users. So, with the ongoing expansion of human cities, an incredibly intricate and vast energy pipeline network (EPN) has been established. Energy providers must invest a significant amount of human, material, and financial resources to ensure the stable operation of the energy supply system, which eventually drives up the costs to the businesses. Fortunately, numerical simulation offers us a tool for cutting costs. Turner et al. [2] believe that applying computer simulation to analyze the hydraulic conditions of extensive pipeline networks is a commercially feasible strategy. In recent years, there have been various computer simulation systems [3,4,5,6] for EPNs, and these simulation tools [7,8,9,10] allow energy supply providers to monitor and troubleshoot the overall network operation status in real-time.

In fact, researches on pipeline simulation have been conducted for many years. The numerical modeling and simulation of pipeline networks have advanced from a single pipe to the entire network. For instance, Taylor et al. [11] completed the validation of a test case consisting of 44.9 miles of pipe in 1962 and offered a computer simulation approach for gas flows in long pipelines. Guo et al. constructed a visual software for graphical and numerical simulations of the transient processes in the water supply network [12], and Wang et al. proposed an efficiently decoupled implicit approach to improve the pipeline flow calculation accuracy in complex natural gas pipeline network simulation [13].

Numerous researchers have proposed methods to improve the simulation model’s correctness, and they have been working to develop an analogous system as a replacement for the actual EPN [14,15,16]. Since the modeling of EPNs in computers involves numerous fluid dynamics computations, these high-precision computational methods require enormous computer resources [17]. Normally, the mass and momentum conservation equations must be solved in all fluid simulation systems [18]. In the simulation of compressible fluid, it is necessary to solve the energy control equations.

Some researchers have found that optimizing the solution process of the pipeline fluid governing equation can reduce the runtime of the entire simulation system [19, 20]; however, the computational acceleration hardware seems to be ignored in their approach, which is another way to increase the execution rate. A computing accelerator, such as Graphics Processing Unit (GPU) or Deep Computing Unit (DCU), is a high-performance calculation device that is specifically designed to perform intensive computing tasks, and most high-performance computers (HPCs) also contain this equipment. With the application of parallel computing technologies, HPC offers an opportunity to enhance the computational efficiency for computationally intensive applications, and it can address the large-scale numerical simulation and data processing issue in academia and industry [21]. Nevertheless, while constructing a parallel strategy, it is important to consider the unique characteristics of the application itself. The HPC concurrency strategy often assigns several similar computing operations to various processors, which can significantly reduce the amount of time that programs take to execute [22].

Task graph is widely applied in load optimization for high-performance parallel computing; it may enhance the parallelism and reduce the complexity of the computational process [23]. So that, EPN is often modeled as a directed graph, and the relationship between vertices is mapped as an adjacency matrix. Theoretically, parallel computation may be accomplished by chunking data such as an image. However, the status of each element in the pipe network has a strong dependency relationship in the traditional pipe network simulation method, which starts from the source node and deduces the state of the whole network progressively in accordance with the network’s boundary conditions. This approach largely limits the parallelization of EPN simulation. However, the Decoupled Implicit Method for Efficient Network Simulation proposed in [13] proves that the solution process of the pipe hydraulic equations is parallelizable. Meanwhile, both Lewandowski [24] and Wang [25] et al. proposed an object-oriented methodology to simulate the energy pipeline system. The authors of these articles agree that the object-oriented programming paradigm is also capable of obtaining high-level accurate EPN simulation models. All of these provide an impetus for us to study the parallelization scheme of pipe network simulation calculation.

In this paper, to overcome the restrictions of graph theory on the parallel hydraulic solution of pipeline networks, we introduce the object-oriented method into the traditional EPN simulation system and focus on the acceleration of urban EPN simulation computation with parallel computing and componentization techniques. The validation and performance of the parallel EPN system are verified with real urban water, natural gas, and thermal transmission and distribution networks. The contributions of this paper mainly include the following:

  • A componentized paradigm for EPN simulation calculation

  • A parallelable simulation computing method for energy pipeline networks which constructs the universal parallel simulation framework (UPSF) for the simulation of energy pipeline networks (EPNs)

  • A practicable parallel computing scheme based on UPSF for the urban EPN simulation

We verify the method’s effectiveness in terms of accuracy and usability by comparing the simulated data with the instrumentation value in real-world working circumstances. The entire experimental process will be presented in the following sections.

2 Related works

Energy supply network in cities is a sophisticated system consisting of a series of pipes and junctions. There are many modeling methods for different application scenarios. For example, a simulation model for studying the dynamic behavior of gas distribution pipeline networks was proposed in [26], and the methods of implicit integration and special crossing-branch suggested in this paper speed up hundreds of equations solved in each time step. Additionally, Kritpiphat et al. designed a simulation model for monitoring and controlling the water pipeline network operations based on the Hardy-Cross and Newton–Raphson methods [27] and completed the experimental validation with a municipal water supply system. A Herrán-González et al. deduced the partial differential equations governing the dynamics of the process, modified the models by adding an inclination term, and then offered two simplified models. Finally, they tested the effectiveness of their approach by simulating a gas distribution pipeline network in MATLAB-Simulink. In addition, Ayad et al. [28] developed a hydraulic simulation model for water pipeline networks based on the extended linear graph theory technique [29]. Their simulation model is equivalent to the EPANET which is a water-distribution-system modeling package developed by the U.S. Environmental Protection Agency Water Supply and Water Resources Division.

All of the aforementioned EPN simulation solutions are based on graph theory which is close to the actual physical network topology [30], and many pipeline analysis tools are also based on this classical method [31, 32]. The intricate connectivity between thousands of components is presented in a connected graph. Researchers firstly build an adjacency matrix and circle-back matrix with the breadth-first and depth-first graph search algorithms; and then, a hydraulic simulation model based on Kirchhoff’s Law is constructed, once the boundary parameters (flow rate, pressure, etc.) of the pipeline network are determined, the hydraulic equations are solved iteratively until the junction pressure and pipe flow errors satisfy the convergence requirements. These kinds of models, which were previously mentioned, are equal to the actual pipe network system since they have strong theoretical backing, and their parameters have undergone extensive verification. However, the hydraulic model evolves from the sources of flow to the terminal nodes in EPN simulation system based on graph theory. Therefore, this type of simulation system usually utilizes a serial programming architecture, and it is difficult to apply the parallel technique to address this issue. This is why most researchers pay attention to EPN modeling and the optimization of fluid simulation solvers and neglect the improvement of the computational performance of the simulation system.

Coincidentally, Lewandowski et al. proposed an object-oriented methodology (OOM) for modeling a natural gas transmission network, which can make it easier to convert the graph of the network structure into an equivalent set of ordinary differential equations [24]. Unfortunately, the author’s descriptions are limited to technical implementation, and it seems to leave out details regarding the system’s performance assessment. Wang et al. applied a hydraulic calculation method based on object orientation to the simulation of a looped-pipeline network with multiple heat sources [25]. They abstracted the heat-supply networks as a collection of various components and regarded each element in the collection as a nondetachable hydraulic calculation unit. On the basis of the fundamental pipe hydraulic calculation technique, this approach changed the computational strategy of the simulation system, which could lessen the hydraulic model’s reliance on the data interdependence among the pipes. Moreover, the effectiveness of OOM was proven again in Wang’s later work [33]. In the following years, Wang et al. also have been optimizing their model for thermal networks [20, 34, 35]; however, they also seem to have overlooked the advantages of OOM in accelerating EPN simulations for large-scale pipe networks which need heavy computing resource to solve more hydraulic equations. In fact, the paradigm based on OOM may be interpreted as a decoupling of several associated computational solution processes. This means that the calculation process-based OOM is parallelizable on a multicore computer. P. Wang has demonstrated its potential on a natural gas pipeline network simulation [13]. The DIMENS method proposed by them could assist in the parallel solving of pipeline processes, and they also plan to apply it to study the GPU-accelerated simulation for large-scale natural gas pipeline networks on computers with acceleration cards.

Calculation acceleration cards might improve the computers’ capability and offer powerful assistance when tackling challenging computational engineering problems. Currently, an growing number of experts from various fields are now working on application speedups with parallel techniques. For example, Zou et al. summarized the processing of genetic data into three approaches, one of which is parallel computing [36], to handle the extensive data processing required for gene sequencing. Yang et al. proposed a parallel computing model to speed up data processing, and they also implemented a service-oriented parallel application in a high-quality network environment to address the challenges of processing large GNSS (Global Navigation Satellite System) datasets [37]. A fine-grained parallel strategy based on the shared memory in [38] balanced the computation and extended the application of the direct simulation Monte Carlo method with open multiprocessing. All of the papers show that parallel computing can indeed facilitate many scientific calculations, and few researches study parallel simulation for urban EPNs. However, it is critical to observe the specific programming mode and data properties when we formulate a parallel scheme for a real application. Task-based parallel strategies proposed in [39] submitted the set of CFD computing tasks to the runtime scheduler that maps each task dynamically to resources; moreover, the authors also applied a directed graph to detach the application program.

To improve the computational performance of the existing simulation of EPN, in this paper, we propose a universal parallel computational framework for EPN simulation on the basic pipeline hydraulic model, drawing inspiration from OOM and the distributed heterogeneous parallel computing framework based on component flow [40]. The simulation process of pipe network components is encapsulated in several objects whose hydraulic solving processes are separated from the simulation program. Through dividing the calculation process of network components into different objects, it is easy to separate the simulation program’s calculation logic from the calculation process, reduce coupling of the pipe network component object’s calculation tasks, and facilitate calculation process parallelization. We also verify the usability of UPSF in different energy supply network scenarios and compare the acceleration ratio for different versions of programs in the water supply network to demonstrate the effectiveness of parallel solutions proposed by us.

3 Methodology and parallel model

As an energy supply network consists of various components, component modeling has been a crucial task in the EPN numerical simulations. This section illustrates the basic computational methods of the EPN components and a parallel model for these components in a simulation system.

3.1 EPN numerical models

In a general pipeline network, pipes are connected by non-pipe components, such as pumps, valves, and two-way, tee, and cross-junctions [41]; the fluid is exchanged between different pipes. A pipeline network has been usually modeled as a directed graph \(G = (V,E)\), where non-pipe components are denoted by set V of vertices in G, and the pipelines constitute the edge-set E, which connect all non-pipe elements. In the proposed approach, the definition of the EPN is simplified to a set constructed with a pipeline or junction of energy supply networks.

According to the functionality, components can be roughly divided into two categories: pipes and non-pipes, and their model conceptions and mathematical equations are discussed in the following. To facilitate the presentation of formulas and models, the symbols are defined in the nomenclature table at the end of the article.

3.1.1 Pipe model

Pipes play an important role in flowable energy transportation. In simulations, the pressure, temperature, and flow rate are vital quantitative indicators, while the length, diameter, and internal roughness of a pipeline are considered boundary conditions in simulation calculations. The physical parameters are consistent with the actual pipeline mode and equivalence relationships. Figure 1 illustrates how a stream is transferred from one end of a pipe to the other end under a combined action of gravity and pressure; S denotes the cross-sectional area, and \(\theta\) is the degree of inclination of the pipe. Moreover, the two flow equations in Fig. 1 represent the variation in fluid pressure and velocity over a length of dx.

Fig. 1
figure 1

Numerical model of pipeline. The fluid passing through the pipeline is discretized into several small segments, and each segment serves as the smallest numerical unit for solving the governing equations

Full size image

Figure 1 shows that the movement of the flow in a pipe is determined by the flow rate, pressure, temperature, and inclination degree of the pipe. All parameters of the pipe can be calculated by several mini-tubes. The detailed computational process and equations are represented in the following.

Unlike a large energy supply network, the pipeline diameter has been typically omitted in the analysis to reduce the complexity of solving fluid equations [13, 20]. The pipe component is formulated as an seven-tuple (\(pipe\_in\), \(pipe\_out\), length, \(\theta\), roughness, pressure, temperature), where \(pipe\_in\) and \(pipe\_out\) represent the inlet and outlet flow of a pipe; length and \(\theta\) are the spatial properties of the pipeline; roughness, pressure, and temperature are physical parameters of the pipe. When a real EPN is simulated, a one-dimensional NS equations can simulate the flow state in the pipe, and the spatial properties of an EPN might be ignored. In other words, the simulation of a pipe in an EPN can be treated as a one-dimensional model.

The flowing objects conform to the principles of conservation of mass, momentum, and energy. Equation 1 is the mass control equation, where \(\rho\) and u are the density and velocity of the fluid. The equation group in Eq. 2 is the Navier–Stokes equations, whose solution is quite time-consuming. \(\frac{Du}{Dt}\) is a material derivative, p is the pressure, f is the body force, and the \(\tau\) is the viscous stress. x, y, z represent the Cartesian spatial coordinates, the velocity components are indicated by u, v,  and w along x, y, and z directions, respectively.

$$\begin{aligned} \frac{\partial \rho }{\partial t}+ \frac{\partial (\rho u_x) }{\partial x} + \frac{\partial (\rho u_y) }{\partial y} + \frac{\partial (\rho u_z) }{\partial z} = 0 \end{aligned}$$
(1)
$$\begin{aligned} {\left\{ \begin{array}{ll} \rho \frac{Du}{Dt} = -\frac{\partial p_x }{\partial x} + \frac{\partial \tau _{xx} }{\partial x} + \frac{\partial \tau _{yx} }{\partial y} + \frac{\partial \tau _{zx} }{\partial z} + \rho f_x\\ \rho \frac{Dv}{Dt} = -\frac{\partial p_y }{\partial y} + \frac{\partial \tau _{xy} }{\partial x} + \frac{\partial \tau _{yy} }{\partial y} + \frac{\partial \tau _{zy} }{\partial z} + \rho f_y\\ \rho \frac{Dw}{Dt} = -\frac{\partial p_z }{\partial z} + \frac{\partial \tau _{xz} }{\partial x} + \frac{\partial \tau _{yz} }{\partial y} + \frac{\partial \tau _{zz} }{\partial z} + \rho f_z \end{array}\right. } \end{aligned}$$
(2)

Thus, the expressions given by Eqs. 1 and 2 might be simplified; the numerical control equations of the pipeline model are given by:

$$\begin{aligned}{} & {} \frac{\partial \rho }{\partial t} +\frac{\partial (\rho v)}{\partial x} =0 \end{aligned}$$
(3)
$$\begin{aligned}{} & {} \frac{\partial (\rho v)}{\partial t} +\frac{\partial (\rho v v)}{\partial x} + \frac{\partial p}{\partial x} +\frac{\lambda }{2} \frac{\rho v |v |}{d} +\rho g sin\theta = 0 \end{aligned}$$
(4)
$$\begin{aligned}{} & {} \frac{\partial }{\partial t} [(e+\frac{v^2}{2} +gE)\rho S]+\frac{\partial }{\partial t} [(h+\frac{v^2}{2} +gE)\rho v S]+q\rho v S = 0 \end{aligned}$$
(5)

The continuity, momentum, and energy equations, which are given in Eqs. 3–5, denote a group of basic governing equations of a pipeline flow derived based on a homogeneous, geometrical. Therefore, the flow status can be estimated anywhere on the pipeline using the aforementioned components. These equations have been widely applied to pipeline simulations [42, 43], and the detailed process of transformations can be found in [13]. At time t, according to the known pressure, flow rate, and fluid physical parameters through the above discrete equations, the motion state of the fluid in the pipeline can be obtained.

3.1.2 Non-pipe model

Non-pipe components of the EPN model represent junctions of different pipes, and they are used to construct the EPN topology. Generally, pipe fittings can be classified into two types according to their functions, connecting joints and pressure-regulating joints. The former mainly includes valves, two-way junctions, tees, and cross-joints, while the latter contains pressurized pumps and pressure-reducing valves. In Fig. 2, the graph vertices are classified into four categories that contain all junction components of the EPN. However, this study formulates them as a unified model, consisting of \(in\_flow\), \(out\_flow\), and parameter conversion functions. Each junction component in the UPSF is described with a triplet \((in\_flow,out\_flow,cnt\_flow)\), where \(in\_flow\) and \(out\_flow\) correspond to the pipes connected with a junction, and \(cnt\_flow\) is the number of pipes connected with the junction.

Fig. 2
figure 2

The general joint model. All types of components in (a) are abstracted into the corresponding models in (c); these individual joints are connected to each other and construct a connected graph such as subgraph (b). In the UPSF, a unified model corresponds to these components; the number of \(in\_flow\), \(out\_flow\) in (d) can be larger than one

Full size image

Because the hydraulic equations are solved in the pipe model, the joint objects in the UPSF are primarily used for calibrating the tube simulation results and verifying the convergence of the model’s accuracy in the proposed simulation approach. This means that the connecting junction instances of the UPSF transmit the input or output parameters between the nearby pipes, compare the errors, and then correct the values of the input or output characteristics of the pipelines using a specific physical model. This is the network parameter synchronizing mechanism, which can accelerate the convergence of simulation model errors. The \(in\_flow\) and \(out\_flow\) of the non-pipe model shown in Fig. 2 are related to the degree of the graph nodes. According to the laws of momentum and mass conservation, the mass of the material flowing through the non-pipe elements is conservative, and the mass constraint equations of the mathematical non-pipe model are given by

$$\begin{aligned}{} & {} \sum _{i=1}^{IN\_CNT} m_{in,i}^{k} - \sum _{j=1}^{OUT\_CNT} m_{out,j}^{k}=0 \end{aligned}$$
(6)
$$\begin{aligned}{} & {} Degree(k) = IN\_CNT + OUT\_CNT \end{aligned}$$
(7)

where \(IN\_CNT\) and \(OUT\_CNT\) are the in- and out-degree of the kth junction, respectively; their sum denotes the degree of a non-pipe node in the graph; m is the mass of the flow; k is the ID of the non-pipe node. As indicated by Eq. 6, there is no mass loss for any substance passing through the junctions.

This study uses the equations of the non-pipe model presented in [13] and modifies the pressure equation by adding a condition. Therefore, the modified pressure constraint equation of the mathematical non-pipe model is given by

$$\begin{aligned} p_{in,i}^{k}-\delta p_{out,j}^{k}=0,\ where\ i \in [1,IN\_CNT],j \in [1,OUT\_CNT],\ and\ \delta >0. \end{aligned}$$
(8)

Pressure is a variable parameter that is associated with the unique function of junctions, and it can be changed when the flow passes the non-pipe components. The expression in Eq. 8 relates to any component having factor \(\delta\) that is a non-negative number. If \(\delta \in (0,1)\), a junction is regarded as a pump; if \(\delta =1\), a junction is considered a common connection element; otherwise, it represents a pressure-reducing valve.

Temperature is another parameter in the mathematical model of non-pipes, and it is related to the materials transported through a pipeline and the function of junctions. For instance, the temperature is regarded as a constant in water supply networks; in a district heating system, the governing equation is much more complex. A detailed discussion of the energy transfer models for heating pipe networks has been provided by Wang et al. in [34, 35, 44]. Unlike water, natural gas is compressible, and its temperature can vary at connections. In [13, 17], the temperature mathematical equations for various junctions were derived. However, the temperature loss at the pipe links can be neglected, and the solving process of junctions can be simplified. Therefore, the temperature simulation results at a pipe’s junction can be calibrated by

$$\begin{aligned} T_{in,i}^{k}= T_{out,j}^{k},\ where\ i \in [1,IN\_CNT],j \in [1,OUT\_CNT]. \end{aligned}$$
(9)

Further, the EPN simulation system requires defining boundary conditions as the initial parameters in the hydrodynamic calculations. Normally, the pressure, flow rate, and temperature of the external components must be known before performing simulation calculations [43,44,45], and they correspond to the supply and consumption in actual EPNs. These conditions are transient, and three functions of time are used in this study to represent them by Eq. 10, where t is time.

$$\begin{aligned} \left\{ \begin{array}{cl} p=P(t)\\ m=M(t)\\ T=T(t) \end{array} \right. \end{aligned}$$
(10)

3.2 UPSF-based programmable model for EPN simulations

According to the component models discussed in Sect. 3.1, an EPN simulation system consists of several pipes and non-pipe elements, so there are extensive hydraulic numerical calculations to conduct. In practice, a discretization technique has been widely used to solve Eqs. 3–5, and numerous methods have been developed to resolve these equations [5, 13, 15, 27]. The solving process of these equations consumes a large number of computational resources. The programmable EPN component model is presented in Fig. 3. The EPN elements are defined using the OOM based on the component classification; thus, the simulation system becomes a set of instances of classes in the UPSF.

Fig. 3
figure 3

The programmable EPN model. This model works on the software architecture of computing nodes, which involves OS kernel and Libraries levels; it also provides a unified computing template for upper-layer application systems. UPSF encapsulates several class templates for the implementation of EPN simulations

Full size image

This subsection discusses how to use parallel techniques to accelerate the EPN simulation calculations. The hydraulic calculation for both pipe and non-pipe elements is encapsulated in programmable components, which splits the EPN simulation from the hydraulic calculation models for pipe networks. Different simulation scenarios, such as water supply networks, gas pipelines, and thermal networks, can be easily converted using different control equations. This computational paradigm can enhance the flexibility of UPSF. In [13, 24, 25], the feasibility of the OOM was validated by the EPN simulations. According to the initial parameters, all of the instances in the EPN simulation system could solve the hydraulic equations independently and synchronize the result of the individual components periodically until the parameter errors of the connected components converge. Therefore, the programmable component model shown in Fig. 3 has the potential to accelerate the simulation calculation on a parallelizable computing platform.

3.3 Parallel strategy for EPN simulation computation

In the traditional EPN simulation methods, a global fluid model is constructed by using the adjacency matrix of the pipe network topology and Kirchhoff’s law. Then, the variables of each network component are solved by substituting the boundary conditions into the set of equations corresponding to the hydraulic model. Finally, a simulation model equivalent to an actual EPN is established. This technique takes the entire pipeline network as a computational unit, which makes it challenging to decompose the numerical calculation process.

Fig. 4
figure 4

Parallel computing architecture for EPN simulations. Every component corresponds to an individual thread in the hydraulic calculation process, and there are many thread groups in the thread pool. The CUDA APIs and MPI (Message Passing Interface) complete the data exchange process through the shared memory between the host and the computational acceleration cards

Full size image

The UPSF provides a method of multi-component synchronous computation for EPN simulations, which makes allows an increase in the execution efficiency of a simulation system on a multi-processor computing platform. The EPN elements are abstracted into independent computing units in the simulation system. However, these individual computing units are interdependent in terms of boundary conditions. For instance, the output of a pipeline represents the input of its neighboring pipelines. The mutual transmission of boundary conditions between different pipes or non-pipe components restricts the parallel process.

This study proposes a parallel system architecture for EPN simulations based on shared memory. Each instance of an EPN simulation system periodically obtains the boundary parameters that are updated by other objects through the shared storage. This interaction mechanism enables the parallel calculation of hydraulic models that interrelate with each other. The parallel computing architecture for the EPN simulations is presented in Fig. 4. The memory is divided into two categories, namely host memory and device memory, and the calculation of each EPN component is executed in a single thread. The parameters of each thread are exchanged periodically through the shared memory of the host.

Generally, the number of available processors is significantly lower than the number of threads depicted in Fig. 4, and the threads responding to the EPN components constrain each other. A resource competitive relationship will arise during the computing process. To improve the computation execution performance, it is crucial to distribute the workload among processors evenly. Because the data dependencies between threads are consistent with the pipeline topology, it is needed to group the computational tasks in the proposed parallel strategy, which enables the interrelated calculation processes to be loaded on the same computational device. The graph partition technique can decrease the coupling degree of the computing devices, and the time spent on the shared memory accesses is consequently reduced. The graph partition method follows the “high cohesion, low coupling” division principle to group various simulation calculation objects in the pipe network and send them to separate processors or computing nodes for completion [23]. We have implemented a full discussion on this subject in [23].

4 EPN simulation system implementation

In this section, we firstly elaborate on the implementation of a serial EPN simulation process and then illustrate the parallel version of the serial program. Both processes of the simulation system are completed with the UPSF proposed above.

4.1 UPSF-based serial EPN simulation

Fig. 5
figure 5

Serial EPN simulation paradigm. The presented diagram is an example of a water supply network simulation, which includes numerous hydraulic equations that need to be solved during the simulation process

Full size image

To verify the feasibility of the UPSF for EPN simulations, we construct a serial simulation program, and its flow diagram is represented in Fig. 5. The first four steps of the flowchart denote the initialization of the program, and they complete the preparation task before the simulation calculation. For instance, the configure file of an EPN contains the connection relationship, roughness, pipeline length, and boundary conditions, and the error accuracy defines the number of iterations and the termination criteria for solving the fluid equations. Then, the equations of the pipeline model are solved, and the values of the parameters, such as pressure and flow rate, are obtained from the connected pipe instances. The errors of the pipe model parameters are compared in the non-pipe component instances. If the error is below the user-specified threshold, the program will be terminated, and the current model will be saved; otherwise, the parameters of the components are modified in Step 8, and the program starts the iterative process.

UPSF provides the templates of the EPN model’s components which is named the programmable component model in the previous section. The simulation system shown in Fig. 5 is decomposed into several interrelated objects whose numerical calculation is organized according to the network topology. Therefore, the UPSF can provide a distributed computing approach for associative computational objects in the EPN simulation calculation, which is the basis of the parallel solution.

4.2 Parallel algorithm of EPN simulation

The discretization and solving of the hydraulic equations, which are performed in Step 5 (refer to Fig. 5), denote the most time-consuming part of the serial program. Given that the calculation processes of each network component are encapsulated in each instance, they are executed separately during the operation of the simulation system. As such, a parallel version of the process presented in Fig. 5 can be obtained by implementing a parameter synchronization mechanism into the serial program.

Each component instance of the EPN is regarded as a minimum computational unit in the UPSF. In Fig. 2d, the solution to complex hydraulic equations contained in each component instance can be obtained independently. Once the boundary parameters, such as pressure, flow rate, and temperature, are obtained from the shared memory, as shown in Fig. 4, the calculation processes corresponding to these objects are performed in separate threads. That is why the simulation calculation processes of all the components in EPN model can be conducted simultaneously. The serial EPN simulation process is modified by importing an additional synchronization parameter, such as the number of DCU and the task grouping result, into the iterative computation, as shown in Fig. 5. Figure 6 is the parallelization flowchart of Fig. 5. The left of Fig. 6 is similar to Fig. 5. In a heterogeneous hardware architecture, these instructions run on the host side to control the flow of the entire program; the right part is the calculation parallel process after the program completes initialization, and the calculation process of all pipe network components is the calculation loaded into different DCUs or GPUs.

Fig. 6
figure 6

Parallel EPN simulation paradigm

Full size image

Every thread should update the calculation result to the shared memory after all hydraulic equations have been solved. The specific steps of the process are presented in Alg. 1.

Algorithm 1
figure a

Parallel algorithm based on UPSF

Full size image

Similar to the serial algorithm, the input of the parallel algorithm given in Alg.1 contains the pressure of sources, the flow rate of terminal users, and the EPN topology configuration, and a simulation model is the algorithm output. In Alg.1, variable \(conf\_file\) represents an EPN configuration file, including all parameters of each network component; \(boundary\_conditions\) relates to the pressure of sources and the flow rate of terminal users; \(err\_threshhold\) is the accuracy of simulation calculations; \(max\_times\) is the preset maximum number of iterations to solve the hydrodynamic equations. In the parallel algorithm, Lines 1–5 denote the initialization of the program, and they mainly include setting the network topology, boundary conditions, and construction of each component instance. The graph partition method is applied in Line 6, and it divides the instances into several groups that guide the distribution of the computational tasks in Line 11. The loop from Line 8 to Line 12 represents a parallel process for the simulation calculation of instances created in Line 2; Line 9 is the selection of a physical computing device, and Line 11 is the assignment of calculation tasks. Once the objects are loaded on a computing device, they are executed automatically. Line 13 shows the result of the synchronization of the hydraulic equations; this step involves accessing the shared memory. The errors are calculated in Line 14, and they are used to determine whether the while-loop continues or terminates. If the simulation error is convergent at the end of the while-loop, the model will be saved as an algorithm output in Line 17.

It should be noted that the component grouping operation is conducted in Line 6 in Algorithm 1. This is a critical stage for parallel policy scheduling in the algorithm. This study uses a graph partition method to cluster the components of the EPN model. The components in the same group can form a subnet that is considered to be an independently computable unit. Then, each unit can be assigned according to the different DCUs’ working statuses. Logically, the parameters of the pipelines between subnets become the boundary conditions of the parallel computable units; the algorithm needs to update them periodically through shared memory or MPI, as shown in Fig. 4.

5 Experiments and comparisons

In this section, the experiments of the EPN simulation system conducted in various scenarios on a high-performance computer are presented. A hybrid-architecture machine used in the experiments combined a CPU with four computational acceleration cards, which is widely deployed on the Songshan supercomputer. The host was equipped with a Hygon C86 7185 32-core processor, 128-GB RAM, and four vega20 GPUs, which denote AMD computational accelerators and are known as DCUs. All of the programs and algorithms were coded in the C programming language and compiled with gcc-7.3.1 and hipcc-4.0.1 on CenterOS-7.6.

The entire experimental program was divided into two parts. The experiments in Sect. 5.1 verify the applicability of the UPSF in the EPN simulation using three actual EPN cases. Three different versions of the simulation system were developed for a water supply network to test the performance of the parallel paradigm in the second part. The comparison and analysis of the experimental results are discussed in the following subsection.

5.1 Model validation

To verify the applicability and generality of the proposed approach for EPN simulations, three distinct practical scenarios in China, including an urban water supply network, a natural gas pipeline network, and a municipal thermal network, were selected. In the energy supply networks, some pressure gauges were intentionally installed to collect the pressure values on the pipes in real-time, and then, the obtained values were compared with the simulated values calculated by the simulation system.

5.1.1 Urban water supply network

This test considered a domestic water supply network of a county town in Jiangxi Province, China. This network contained 1875 pipelines and 1876 non-pipe components and included two water supply sources. There were 115 pressure gauges deployed throughout the network. The real-time pressure data were collected at 7:30 AM on July 10, 2021, 3:30 PM on July 12, 2021, and 10:30 AM on July 16, 2021. Then, the same boundary conditions were imported into the simulation system, and the pressure of the corresponding pipeline was compared with the values obtained from the pressure measuring equipment. The comparison results are presented in Fig. 7.

Fig. 7
figure 7

Comparison of pressure deviations under different conditions

Full size image
Fig. 8
figure 8

Comparison of pressure deviations in the water supply network under different conditions. The three dashed lines of different colors represent the deviation rate between the simulated pressure values and the measured values for each condition. The thicker blue curve represents the average deviation rate of pressure across the three conditions (color figure online)

Full size image

According to the deviation rate (DR) formula given by Eq. 11, it is determined how much the simulated values deviated from the real measured values.

$$\begin{aligned} DR(t) = \frac{|Sim\_V^{k}_{t} - Mea\_V^{k}_{t} |}{Mea\_V^{k}_{t}} \times 100\% \end{aligned}$$
(11)

In Eq. 11, DR is the function of t. Different moments relate to various boundary conditions. k is the ID of a gauge; \(Sim\_V\) and \(Mea\_V\) denote the simulated and measured values of the pipe pressure, respectively. \(Sim\_V^{k}_{t}\) and \(Mea\_V^{k}_{t}\) represent both pressure and temperature. Therefore, their DRs can be calculated using Eq. 11.

The pressure deviations under the three conditions are represented in Fig. 8, where it can be seen that the DRs fluctuate around 6%, and the pressure error deviations of most observation pipelines oscillated approximately 4% among different conditions. Since leakages and other failures are common in most practical urban water supply pipe networks [46], the pressure errors of the simulation result are unavoidable. If the failure rate of the network is kept at a low level, the DR value is likely to decrease.

5.1.2 Natural gas pipeline network

This case considered an urban gas distribution network in Zhejiang Province, China. The total length of the pipeline was 42.37 km; this contained a single source station, 131 valves, 69 plugs, 47 pressure regulators, 1513 two-way junctions, and 115 tee junctions. Similar to the water supply network, 31 pressure devices were used to obtain the gas pressure value, and simulated were compared with the values measured on Sep. 26, 2021, at 11:30 AM, Sep. 6, 2022, at 6:00 AM, and Sep. 20, 2022, at 8:00 PM. The results are presented in Fig. 9.

Fig. 9
figure 9

Comparison of gas pressure deviations under different conditions

Full size image

As shown in Fig. 9, most of the measurements at the observation pipelines were consistent with the values obtained by the simulations. These curves can demonstrate that the UPSF is applicable in the computational simulation of gas pipe network models.

Fig. 10
figure 10

Comparison of gas pressure deviations in the natural gas pipeline network under different conditions. The three dashed lines of different colors represent the deviation rate between the simulated gas pressure values and the measured values for each condition. The thicker blue curve represents the average deviation rate of gas pressure across the three conditions (color figure online)

Full size image

The DR curves are displayed in Fig. 10. Although the amplitudes of these curves are relatively large, their maximum errors still remain below 4.5%. As shown in Fig. 10, it can be observed that the DRs at 08:00 PM on Sep. 20, 2022, and 11:30 AM on Sep. 26, 2022, consistently remain below 1%, while the average deviations between them were below 2%.

5.1.3 Municipal thermal supply network

This case selects a municipal heating network of a district in northwestern China. The total pipeline length of this EPN is 190 km, and there is a heat source station, along with 136 thermal exchange stations, 205 valves, 44 plugs, 38 two-way junctions, 262 tee junctions, and 26 cross-junctions. Due to the bidirectional pipelines and heat loss, the structure of a thermal network is more complex than those of the other EPN types.

We deployed 75 measuring devices to gather the pressure and temperature of the pipes. The simulated results comparison of Feb. 23, 2022, namely 2:50 AM, 12:50 PM, and 9:50 PM, are presented in Figs. 11, 13, and 15, and the DRs comparison are displayed in Figs. 12, 13 and 14 respectively.

Fig. 11
figure 11

Comparison of the pressure in the thermal network under different conditions

Full size image
Fig. 12
figure 12

Comparison of pressure deviations in the thermal network under different conditions. The three dashed lines of different colors correspond to three different boundary conditions, and most DRs are near 5% (color figure online)

Full size image
Fig. 13
figure 13

Comparison of the supply water temperature in the thermal network under different conditions

Full size image
Fig. 14
figure 14

Comparison results of the supply water temperature in the thermal network under different conditions. Most DRs of the three conditions are less than 2%, and the mean DR is less than 1.5%

Full size image
Fig. 15
figure 15

Comparison of backwater temperature in thermal networks under different conditions

Full size image
Fig. 16
figure 16

Comparison results of the backwater temperature in the thermal network under different conditions. Most DRs of the three conditions are less than 2%, and the mean DR is less than 1.5%

Full size image

The curves in Figs. 11, 13, and 15 maintain a high level of fitting. Despite there are minor oscillations observed in certain data points of the curves in Figs. 13 and 15, the temperature relative deviation is negligible. This can be verified by DRs in Figs. 14 and 16. In the UPSF-based simulation of thermal network, the mean DR of pressure fluctuates around 6%, while the DRs of supply and back water temperatures remain below 2%, which can be seen as a positive consequence. These findings confirm the feasibility of the simulation framework we proposed in practical applications.

Overall, the comparison in Figs. 11, 13 and 15 show that the UPSF-based EPN simulation systems are accurate and available in multiple scenarios. The temperature DRs in Figs. 14 and 16 are below 2%, except for pressure DRs which oscillates around 5%. It should be noted that the DRs in Figs. 8 and 12 are higher than the DR of gas pressure in Fig. 10. That maybe resulted by some undetected leaks in the water pipelines network. However, the actual measured pressure values are closer to the values calculated in the simulation in Fig. 9 due to the high-level security in natural gas network.

5.2 Performance of parallel simulation for water supply network

The viability of the UPSF was tested in different scenarios in the previous subsections. In this section, the feasibility of the parallel approach for the EPN simulation is demonstrated using a model of the water supply network.

5.2.1 Test cases

Since the architecture of the computing node contained one CPU and four DCUs, the simulation program of the water supply network was conducted using different paradigms to verify the performance of the parallel simulation approach. The following three versions of the water EPN simulation system were used in these experiments.

  • Serial program (SP)

    This test case was consistent with the flowchart shown in Fig. 5. The hydraulic equations of each EPN component were solved in a certain order. In other words, these series of processes were carried out by a single thread;

  • Parallel program with multi-thread (MTPP)

    This case was based on multi-threading technology, where all calculation tasks were distributed to different cores of a computer in a multi-thread scheme, and the inter-thread communication technology could help individual network components synchronize data. The purpose of this experiment was to demonstrate the efficiency of the parallel component computing strategy in accelerating the EPN simulation computation speed;

  • Parallel program with GPU (GPUPP)

    The parallel program with a GPU used the robust processing capacity of GPUs and solved massive hydrological equations of parallel EPN components on a particular acceleration device. The viability of the proposed parallel approach was examined in this experiment.

The three programs were applied to networks of various sizes. Four water supply network models were used in the comparison experiments, and their parameters are listed in Table 1. Due to the differences in the types of network components, Table 1 shows the numbers of pipes and non-pipes. It should be noted that the computational workload of the simulation system would increase with the number of components.

Table 1 The water EPN models used in the tests
Full size table

Additionally, to unify the evaluation criteria of the experiments, the error threshold of all computational models in the experiments was set to \(10^{-6}\), and \(max\_times\) was set to \(3 \times 10^{5}\), both of them were used to determine whether the simulation model calculations had converged. Each network model was applied to a different version of the simulation program, and the executing times and acceleration ratios of all testing cases were compared, which is explained in detail in the following subsection.

5.2.2 Running times of different paradigms for water simulation system

To ensure the reliability of the program’s running time, all experiments were conducted in a way that exclusively occupied the entire computing machine. Table 2 compares the running times of different experimental programs. The results in Table 2 show that the calculation time of MTPP and GPUPP is significantly better than SP. Figure 17 can reflect their differences more vividly.

Table 2 Running times comparison among different paradigms
Full size table

The running time of the serial version is represented in Fig. 17a, where it can be seen that the duration of the program’s execution increased with the network size.

Fig. 17
figure 17

Running times of different simulation calculation paradigms of the four net models

Full size image

Due to the restriction of the number of cores in the computer hardware, the MTPPs using 4, 8, 16, 24, and 32 threads were launched simultaneously; the results are presented in Fig. 17b. The multi-thread version of the water supply network simulation performed significantly better in terms of running time than the SP. However, the result of the Net1 model was an unexpected outcome, which was caused by the fewer network components in the Net1 model. The synchronization of data between threads increased the running time of the program. An unanticipated result in Fig. 17b indicates that calculations performed with 32 threads took longer than those with 16 or 24 threads, which was due to the grouping of components. When the number of groups was too large, the data synchronization between threads increased the time consumption of the overall process. Thus, the case with 32 threads was not the fastest MTTP solution. Based on the results in Fig. 17b, the partition of the ENP simulation calculation tasks could optimize the calculation performance, which will be further discussed in future work.

The running time of the GPUPP is presented in Fig. 17c. All parallel testing cases with a GPU outperformed both the SP and the MTPP, which was mostly attributable to the GPU’s potent capability for parallel computation.

5.2.3 Comparison and analysis of acceleration ratio

The parallel approaches based on the UPSF performed better than the traditional serial version. This section discusses and compares the acceleration ratio (AR) results of different simulation calculation programs. The AR was calculated by

$$\begin{aligned} AR(pro_1,pro_2) =\frac{runtime\_pro_1}{runtime\_pro_2} \end{aligned}$$
(12)

Three ARs were calculated: SP versus MTPP, SP versus GPUPP, and MTPP versus GPUPP, as presented in Table 3. The running time in the third column of Table 2 represents the average result of each group in Fig. 17b. The running time in Table 2 confirms the UPSF’s practicability in a parallel computing simulation of EPNs. Although the ARs of SP versus MTPP (0.67) and SP versus GPUPP (0.85) were less than 1.0, that was an isolated case, which could be mainly attributed to the small amount of computation and the execution time wasted on data synchronization between threads.

Table 3 Comparison results of the AR for different paradigms
Full size table

However, compared to Net1, the rest ARs of the other networks performed relatively well, which was due to the larger network. The results in Table 3 show that the MTPP and GPUPP were significantly superior to the SP in terms of calculation time, and this advantage increased with the network size. Thus, the UPSF could provide an effective way to address the problem of accelerating the simulation calculation for large-scale EPNs, which was proven by the AR of SP versus GPUPP and AR of MTPP versus GPUPP for the Net4. When the size of EPN becomes enormous, the computing load on the simulation program increases, and the acceleration effect of GPUPP is more obvious. Compared to the running time of SP, the Net4's AR of SP versus GPUPP is 18.09, and the AR of MTPP versus GPUPP also reaches 4.86.

Because the GPUPP used four DCUs to complete the simulation of each network and considering the different numbers of physical cores in each test case, the runtimes of SP were compared with average runtimes on each core or DCU, as shown in Table 4. The average time was calculated by Eq. 13, where \(runtime\_pro\) is the program runtime, and \(core\_cnt\) is the core number of a physical computer.

$$\begin{aligned} Avg\_time = \frac{runtime\_pro}{core\_cnt} \end{aligned}$$
(13)

The running time of the SP was related to the fact that the program was executed on a single processor, as shown in Table 4. The results in Table 4 show that the acceleration effect of the proposed parallel scheme was relatively obvious and demonstrates the practicality of the UPSF for AMD DCUs.

Table 4 Comparison of average ARs among various paradigms
Full size table

6 Conclusion and future works

This study investigates the modeling and simulation methods for energy supply networks, modifies the non-pipe model for synchronizing pipeline status, and proposes a universal parallel simulation framework. It also carries out a series of simulation experiments of different EPN scenarios based on the UPSF with a high-performance DCU+CPU computing platform.

We compare the simulated data with the instrumentation value in the actual working scenarios, which can demonstrate the applicability of the proposed UPSF-based simulation system. Most deviation rates of each experiment are less than 5%. The experimental results of the three types of paradigms of the simulation system for a water supply network show that the parallel version of the program can achieve an 18-fold improvement in the execution time compared to the serial version. This demonstrates that the UPSF-based parallel approach exhibits a significant advantage in terms of computational speed for the simulation of large-scale EPNs.

Although the proposed parallel approach can achieve a better acceleration ratio compared with the serial paradigm, the HPC performance has not been completely exploited in this study. Therefore, simulations on larger-scale EPNs should be conducted in future work using the proposed method. In addition, it is necessary to study additional ways to distribute the tasks of solving hydraulic equations between different calculation accelerators. Finally, task partitioning between GPUs could be an essential task in future work.