Research on Seismic Connectivity Reliability Analysis of Water Distribution System Based on CUDA

: To improve the seismic connectivity reliability (SCR) analysis efﬁciency of water distribution systems (WDS) based on Monte Carlo (MC) simulation, the quasi-Monte Carlo (QMC) method sampled by a low-discrepancy sequence is applied. Furthermore, a parallel algorithm combined with the breadth-ﬁrst search algorithm for SCR analysis of WDS based on the QMC method and Compute Uniﬁed Device Architecture (CUDA) platform was proposed. A city WDS was taken as a computational example, the accuracy and efﬁciency of the traditional MC algorithm and parallel algorithm were compared, and the inﬂuence of the Sobol sequence and pseudo-random number sequence was analysed. The analysis results show that when 1,000,000 simulations are performed, the maximum error of the calculation results of the two methods is 0.2%, and the parallel method can obtain a six-fold speedup ratio compared with the serial method, indicating that the proposed parallel method is correct, meets the accuracy requirements, and helps to improve the SCR analysis efﬁciency. When the number of simulations is the same, the simulation results based on the Sobol sequence are more accurate than those based on the pseudo-random number sequence. The proposed parallel method also achieves a good acceleration effect in the SCR analysis of large-scale WDS.


Introduction
Urban water supply systems, gas supply systems, power supply systems and other infrastructure systems are the lifeblood to maintain the economic and social functions of the city, and people call them lifeline projects. With the rapid development of the national economy, the scale of the city increases, the degree of modernisation is increasingly high, and the dependence on the lifeline project is increasingly strong. However, in increasingly frequent natural disasters, especially under the action of a strong earthquake [1][2][3], lifeline projects suffered serious damage, and water supply, power supply, communication and other services were interrupted, which caused secondary disasters, and even led to paralysis of the function of the whole city, resulting in a large number of life and property losses. The water distribution system (WDS) is an indispensable infrastructure in urban lifeline projects, with the characteristics of wide spatial distribution and strong damage correlation, which has been severely damaged in previous destructive earthquakes, resulting in partial or complete failure of the urban water supply function and seriously affecting people's lives [4][5][6][7][8][9]. Therefore, it is necessary to carry out seismic connectivity reliability (SCR) analysis to understand the ability of a WDS to perform its function in the first instance and provide scientific support for the pre-earthquake reinforcement and post-earthquake reconstruction of a WDS, and as a preliminary analysis of the water network resilience. The SCR analysis of a WDS involves researching its seismic reliability at the system level and analysing the connection probability between user nodes and water sources [10][11][12].
At present, research on SCR analysis of WDS is quite mature. However, with the continuous progress of urbanisation, the scale of the urban WDS is increasing. How to improve the efficiency of SCR analysis of WDS is the focus of current research. In the present paper, combining graph theory and the QMC simulation technique, an efficient method for the SCR analysis of WDS is proposed. Furthermore, a CUDA-based parallel algorithm is developed to further improve the SCR analysis efficiency of WDS.

Graph Theory
Graph theory is a branch of mathematics that deals with graphs [41]. A graph is composed of a set of vertices and a set of binary relations between vertices. It is often used to describe a specific relationship between things, where vertices represent things and edges represent the relationship between corresponding things. Actually, it is an important data structure that can appropriately describe or model problems in many areas of the natural and social sciences. To use a computer program to analyse and study a graph, the graph must be stored in some form in the computer memory. There are many ways to store graphs on a computer, including the adjacency matrix [42], incidence matrix [43], and adjacency list [44]. The WDS is a complex network system composed of water supply sources (water plants, water wells, highland pools, etc.), nodes (concentrated water points, user demand terminals, etc.), and pipelines. In essence, it can be simplified as a network topology composed of points and edges, and its mathematical model can be represented by a graph. It should be noted that in a WDS the flow in the pipe could be unidirectional or bi-directional, thus the WDS can be abstracted as a directed graph.

QMC Method Principle
The basic idea of the QMC method [10,45] is to establish a probabilistic model or stochastic process related to the solution of the problem, whose parameters are equal to the solution of the problem; then, the statistical characteristics of the desired parameters are calculated through the observation or sampling test of the probability model or random process; finally, the approximate solution of the problem is given according to the requirements of solution accuracy or test times. Compared with the traditional MC method, the key of the QMC method is to use LDS to sample and simulate the problem. For example, when the QMC method is used to solve the integral on the unit hypercube I s = [0, 1] s , the calculation formula can be written as follows: where N is the number of sampling times, x k is the random of the kth sampling and x = {x 1 , x 2 , . . . , x N } is LDS. According to the Koksma-Hlawka inequality, the order of the error of the QMC integral can be obtained as follows: where D N * is the discrepancy degree of the LDS, which can be defined as follows: where V is any hypercube contained in I s , U (x 1 , x 2 , . . . , x n ) is the number of points in LDS{x 1 , x 2 , . . . , x N } contained in region V, and λ(V) is the Lebesgue measure of region V. Equation (2) gives the upper limit of the integral error of the QMC method, indicating that the smaller the discrepancy degree of the sequence used is, the more accurate the result of integration. For an LDS containing N points, the discrepancy degree of the sequence is O(N −1 (logN) s−1 ), and the order of the error convergence rate using the QMC method is O(N −1 ). However, the order of the error convergence rate using the MC method is O(N −1/2 ), which shows that the calculation result of the QMC method has higher accuracy than that of the MC method when the same sampling times are sampled [46].

Sobol Sequence
The construction of the LDS is the key content of the QMC method. The advantage of the LDS is that the samples generated tend to be more evenly distributed. At present, the commonly used LDS include the Halton sequence [47], Sobol sequence [48] and Latin hypercube sequence [49]. In this study, the Sobol sequence is used to sample the failure probability of WDS pipelines. The construction method of the Sobol sequence is briefly introduced below.
Suppose m i is a positive odd number less than 2 i , then the direction number d i is calculated as follows: where w is the number of direction numbers; d i is generated by a simple polynomial, which can be expressed as follows: where a q equals 0 or 1. When i > q, d i can be expressed by the following recursive formula: where ⊕ represents the bitwise XOR of the binary system. For example, 0 ⊕ 1 = 1 ⊕ 0 = 1, For m i , the equivalent recurrence formula is expressed as follows: Then the kth number in the Sobol sequence is generated as follows: where (· · · c 3 c 2 c 1 ) is the binary representation of k.
To improve the generation speed of Sobol sequences, Antonov and Saleev proposed an improved calculation method [50], which uses the following recursive: where · · · g 2 (k)g 1 (k)g 0 (k) is a binary representation of gray code with a variable k − 1. j is equal to the smallest n that satisfies g n (k) = 0. Figure 1a,b show the two-dimensional distribution of the pseudo-random sequence and Sobol sequence sampled 2000 times in the range of 0 to 1, respectively. It can be seen that the sampling points of the pseudo-random sequence are prone to local clustering, while the Sobol sequence has better uniformity than the pseudo-random sequence.

SCR Analysis of WDS Based on the QMC Method
The basic idea of using the QMC method to analyse the SCR of a WDS is to gene a large number of quasi-random numbers and compare them with the failure probab of pipelines to determine the failure state of each pipeline in the WDS. In each simula the connectivity between the water supply source and demand nodes is judged. Then frequency of each node in the connected state after multiple simulations is calculated, the probability evaluation of its connectivity reliability is obtained.
When analysing the SCR of a WDS, it is ordinarily assumed that the water su source and demand nodes will not be damaged under earthquake action. Without con ering the change in water pressure at the node after the earthquake, the water in the always flows from the high hydraulic pressure side to the low hydraulic pressure sid this case, there is only one unique path between two nodes. Consequently, the topo of the WDS can be abstracted as a directed graph with edge weights. In this article adjacency table is used as the storage method of the pipe network, and the breadth search (BFS) algorithm [51] is used to judge the connectivity of the pipe network. Ass ing that there are m nodes in the WDS, the steps of SCR analysis based on the Q method are as follows: 1. According to the seismic vulnerability analysis of pipelines, the damage probab Pij of each pipeline in the WDS under seismic action is determined, see Equation 2. Generate a low-discrepancy random number Rij between 0 and 1 and compare it the damage probability Pij of the corresponding pipeline. If Rij ≥ Pij, the pipelin considered to work normally; otherwise, the pipeline is considered to be ineffec The pipeline status aij is shown in Equation (11), the numbers 0 and 1 represen pipeline in "failures state" and "reliable state", respectively. Based on this met the connection state of all pipelines can be obtained, and then the adjacency tab the pipe network can be generated.
where i and j are upstream and downstream nodes and m is the total numbe nodes. 3. The BFS algorithm is used to traverse the adjacency table of the pipe network, the water source point is set as the starting point. First, put the water source p into queue Q, and search for all the nodes directly adjacent to the water source p Then, put the adjacent nodes of the water source point into queue Q, and delete first node in the current queue Q, that is, the water source point is out of the qu Next, continue searching with the first vertex in queue Q as the starting point,

SCR Analysis of WDS Based on the QMC Method
The basic idea of using the QMC method to analyse the SCR of a WDS is to generate a large number of quasi-random numbers and compare them with the failure probability of pipelines to determine the failure state of each pipeline in the WDS. In each simulation, the connectivity between the water supply source and demand nodes is judged. Then, the frequency of each node in the connected state after multiple simulations is calculated, and the probability evaluation of its connectivity reliability is obtained.
When analysing the SCR of a WDS, it is ordinarily assumed that the water supply source and demand nodes will not be damaged under earthquake action. Without considering the change in water pressure at the node after the earthquake, the water in the pipe always flows from the high hydraulic pressure side to the low hydraulic pressure side. In this case, there is only one unique path between two nodes. Consequently, the topology of the WDS can be abstracted as a directed graph with edge weights. In this article, an adjacency table is used as the storage method of the pipe network, and the breadth-first search (BFS) algorithm [51] is used to judge the connectivity of the pipe network. Assuming that there are m nodes in the WDS, the steps of SCR analysis based on the QMC method are as follows: 1.
According to the seismic vulnerability analysis of pipelines, the damage probability P ij of each pipeline in the WDS under seismic action is determined, see Equation (14).

2.
Generate a low-discrepancy random number R ij between 0 and 1 and compare it with the damage probability P ij of the corresponding pipeline. If R ij ≥ P ij , the pipeline is considered to work normally; otherwise, the pipeline is considered to be ineffective. The pipeline status a ij is shown in Equation (11), the numbers 0 and 1 represent the pipeline in "failures state" and "reliable state", respectively. Based on this method, the connection state of all pipelines can be obtained, and then the adjacency table of the pipe network can be generated.
where i and j are upstream and downstream nodes and m is the total number of nodes.

3.
The BFS algorithm is used to traverse the adjacency table of the pipe network, and the water source point is set as the starting point. First, put the water source point into queue Q, and search for all the nodes directly adjacent to the water source point. Then, put the adjacent nodes of the water source point into queue Q, and delete the first node in the current queue Q, that is, the water source point is out of the queue. Next, continue searching with the first vertex in queue Q as the starting point, and repeat the above steps until queue Q is empty. In this process, it should be stated that the vector Z = [z i ] (i = 1, 2, . . . , m) is used to record the access status of the nodes. For node i, if it was ever put into queue Q, the value of z i is set to 1, indicating that node i is connected with the water source point, otherwise it is set to 0. For a multi-source WDS, a virtual source point can be added, and assuming that the failure probability of the pipeline connecting the virtual source point to each water source point is 0, the multi-source system can be transformed into a single-source WDS. 4.
The vector Z obtained from each simulation is added up and the result is saved in the vector T initialised to 0:

5.
Repeat steps (2)-(4) K times to obtain the connectivity probability vector P c = {P c1 , P c2 , . . . , P cm } between the water source point and each node, P ci can be calculated as follows:

Parallel Algorithm for SCR Analysis of WDS Based on CUDA
With the increase in simulation times and pipe network scale, the calculation time of SCR analysis of WDS increases significantly. Therefore, combined with GPU general computing technology, a parallel algorithm based on CUDA for SCR analysis of WDS was designed in the current article.

Parallel Scheme Design
It can be seen from Section 3.3 that the SCR analysis of the WDS requires K QMC simulations. Since each QMC simulation is carried out independently and does not interfere with each other, it has natural parallelism. In addition, in each QMC simulation, when the BFS algorithm is used to calculate the connectivity of the pipe network graph, the adjacent nodes of each layer of vertices can be searched in parallel by multiple threads to shorten the time of traversing the graph. Therefore, this article uses multithreading to process K QMC simulations and BFS traversals of graphs in parallel. The algorithm is based on CUDA using GPU-CPU heterogeneous computing architecture, divided into a CPU processing part and a GPU processing part. The CPU processing part includes the steps of importing basic data of the WDS, generating the adjacency table of the pipe network, transferring data between the host and the device, and allocating the memory of the host and the device. The GPU processing part processes K times for parallel QMC simulations, BFS parallel traversals, and LDS generation to achieve an efficient solution of the SCR analysis of WDS. The logic flow chart of the algorithm of CUDA-based SCR analysis of WDS is shown in Figure 2. The main flow of the algorithm is as follows, for detailed algorithm implementation, see Section 4.2.3.
(1) Input the basic data such as the total number of nodes m, the total number of pipelines n, the starting point and ending point number of the pipeline, the failure probability of the pipeline under earthquake, the QMC simulation times K in the CPU, and store the pipe network diagram in the format of compressed row storage (CRS) [52]. (2) Use the "cudaMalloc()" function provided by CUDA to allocate the space in the global graphics memory, and copy the data (i.e., pipe network diagram, number of nodes and pipelines in the pipeline network) from the host memory to the graphics memory through the "cudaMemcpy()" function. (3) Generate a uniformly distributed Sobol sequence between 0 and 1 on the GPU through the "curandCreateGenerator()" and "curandGenerateUniform()" functions, containing n*K random numbers. Then, the Sobol sequence is compared with the failure probability of pipelines, and the CRS of K graphs is generated. (4) Start the kernel function "bfsVisit()" and implement the BFS traversal of k pipeline network graphs in parallel on the GPU. It should be noted that the kernel function "bfsVisit()" needs to be restarted repeatedly until all the valid nodes are traversed, each time it only searches the valid adjacency nodes of all valid nodes in a certain layer of the pipeline network graph in parallel. Assuming that the adjacency edge of a valid node at the ith layer of the graph is in a reliable state and the corresponding node of the adjacency edge has not been traversed, this node is called a valid adjacency node. The valid adjacency node is defined as the valid node at the (i + 1)th layer of the graph. (5) Start the kernel function "avgVisitedCuda()" to calculate the final connection probability of each node in parallel on the GPU. (6) Copy the connectivity probability of each node to the CPU and output it.
Water 2023, 15, x FOR PEER REVIEW 7 of 18 (4) Start the kernel function "bfsVisit()" and implement the BFS traversal of k pipeline network graphs in parallel on the GPU. It should be noted that the kernel function "bfsVisit()" needs to be restarted repeatedly until all the valid nodes are traversed, each time it only searches the valid adjacency nodes of all valid nodes in a certain layer of the pipeline network graph in parallel. Assuming that the adjacency edge of a valid node at the ith layer of the graph is in a reliable state and the corresponding node of the adjacency edge has not been traversed, this node is called a valid adjacency node. The valid adjacency node is defined as the valid node at the (i + 1)th layer of the graph. (5) Start the kernel function "avgVisitedCuda()" to calculate the final connection probability of each node in parallel on the GPU. (6) Copy the connectivity probability of each node to the CPU and output it.

Storage of Graph on the GPU
Storing graph data in the GPU requires three considerations: first, the relatively limited capacity of the GPU's video memory compared to the host's memory; second, CUDA uses a CPU-GPU heterogeneous computing model in which data are transferred back and forth between the host side and the device side; third, in the CUDA storage system, device memory and host memory have different address spaces, which makes it inconvenient to operate pointer data on the CUDA platform. In addition, since the device memory is Storing graph data in the GPU requires three considerations: first, the relatively limited capacity of the GPU's video memory compared to the host's memory; second, CUDA uses a CPU-GPU heterogeneous computing model in which data are transferred back and forth between the host side and the device side; third, in the CUDA storage system, device memory and host memory have different address spaces, which makes it inconvenient to operate pointer data on the CUDA platform. In addition, since the device memory is abstracted as linear memory, the data structure in the form of an array has the highest access efficiency. Therefore, the CRS format is used to store the pipeline network diagram in the parallel algorithm in the current article. Figure 3c shows the CRS format of the directed edge-weighted graph in Figure 3a. The CRS format uses three arrays to store directed edgeweighted graphs, val array stores the values of nonzero elements in a row-wise fashion, col_ind array stores the corresponding column indices of the elements in the val array, and row_ptr array stores the locations in the val array and col_ind array that start a row. access efficiency. Therefore, the CRS format is used to store the pipeline network diagram in the parallel algorithm in the current article. Figure 3c shows the CRS format of the directed edge-weighted graph in Figure 3a. The CRS format uses three arrays to store directed edge-weighted graphs, val array stores the values of nonzero elements in a rowwise fashion, col_ind array stores the corresponding column indices of the elements in the val array, and row_ptr array stores the locations in the val array and col_ind array that start a row.

Sobol Sequence Generation
In the simulation, a large number of quasi-random numbers are used. The cuRAND function library in the CUDA platform provides a pseudo-random number generator and a quasi-random number generator, which can be used to generate random numbers on the host and device. Sobol sequence is the only quasi-random number sequence supported by the host and device API [53]. Therefore, the Sobol sequence is used as the LDS for QMC simulation in this study.

Implementation of the CUDA-Based Reliability Algorithm for WDS
The essence of the parallel algorithm in this study is to use multiple threads to perform multiple BFS traversals on the pipeline network graph at the same time. The pseudocode of the algorithm executed on the host is shown in Algorithm 1. Three arrays, Ra, Ca and Va are created to store the pipeline network graph G(V, E, W). A float array Sa with the size of eNum * Sn is used to store low-discrepancy random numbers, and an int array Pa and another float array Na with the size of vNum * Sn are used to store the seismic reliability of each node and the number of layers where the node is located in the pipe network, respectively. First, the number of layers of the source node in each graph is set to 1 through the kernel function "initVisit()" (see Algorithm 2), and then the source nodes are the starting point when traversing pipe network graphs, and the number of layers of other nodes is set to 0. Next, initialise the Boolean variable F to false and the integer variable Ct to 1. In each iteration, it is necessary to determine whether the value of the flag variable F

Sobol Sequence Generation
In the simulation, a large number of quasi-random numbers are used. The cuRAND function library in the CUDA platform provides a pseudo-random number generator and a quasi-random number generator, which can be used to generate random numbers on the host and device. Sobol sequence is the only quasi-random number sequence supported by the host and device API [53]. Therefore, the Sobol sequence is used as the LDS for QMC simulation in this study.

Implementation of the CUDA-Based Reliability Algorithm for WDS
The essence of the parallel algorithm in this study is to use multiple threads to perform multiple BFS traversals on the pipeline network graph at the same time. The pseudocode of the algorithm executed on the host is shown in Algorithm 1. Three arrays, R a , C a and V a are created to store the pipeline network graph G(V, E, W). A float array S a with the size of eNum * S n is used to store low-discrepancy random numbers, and an int array P a and another float array N a with the size of vNum * S n are used to store the seismic reliability of each node and the number of layers where the node is located in the pipe network, respectively. First, the number of layers of the source node in each graph is set to 1 through the kernel function "initVisit()" (see Algorithm 2), and then the source nodes are the starting point when traversing pipe network graphs, and the number of layers of other nodes is set to 0. Next, initialise the Boolean variable F to false and the integer variable C t to 1. In each iteration, it is necessary to determine whether the value of the flag variable F is false, and if so, invoke the kernel function "bfsVisit()" (see Algorithm 3). Then, the flag variable F is set to true, and the nodes whose layer numbers are equal to C t are treated as the new starting point to traverse the next layer of the pipe network graph. If the connectivity probability of the adjacent edges of the nodes of the current layer is greater than the corresponding low-discrepancy random number, and the number of layers of the adjacent nodes is 0, then the layer of the corresponding adjacent nodes is set to C t +1, and flag variable F is set to false. It should be noted that the above process is repeated several times until F is true. Finally, the kernel function "avgVisitedCuda()" is activated to calculate the seismic reliability of each node (see Algorithm 4). 1: Create location array R a , column indices array C a and weight array V a from G(V, E, W) 2: Create Sobol random numbers array S a of size eNum * S n , node traversal array N a of size vNum * S n , connectivity probability array P a of size V 3: Invoke CUDA initVisit (N a , vNum, S 0 ) on the grid 4: F ← Flag ← false, C t ← current_Traversal ← 1 5: while F is false 6: Invoke CUDA bfsVisit (R a , C a , V a , S a , C t , N a , F, vNum, eNum) on the grid 7: C t ← C t + 1 8: end while 9: Invoke CUDA avgVisitedCuda (P a , N a , vNum, S n ) on the grid

Numerical Example
To verify the applicability of the parallel method in this paper, a WDS is taken as an example to analyse its SCR, and the layout of the WDS is shown in Figure 4, which consists of 40 nodes and 64 pipelines. The water supply source is node 40, whose water level is constant at 60 m. The demand flow rate and ground elevation of other nodes can also be seen in Figure 4. Ductile cast iron pipes are adopted in this WDS, and information on pipeline length and diameter is shown in Figure 5. To obtain the probability of different failure states of the pipeline, three kinds of analysis methods are generally adopted, in-cluding empirical methods, expert judgment methods and analytic methods. In this study, an empirical method is adopted, and its core idea is to obtain the pipe repair rate (RR) considering pipe type, pipe diameter, site conditions, soil characteristics and other factors by statistical regression of historical earthquake damage data. Furthermore, assuming that the pipeline seismic damage occurs randomly and independently along the pipeline and follows the Poisson distribution, the failure probability P d of the pipeline can be expressed as follows: where P d is the failure probability of the pipeline, including the probability of moderate damage and severe damage; RR is the pipe repair rate, expressed as repairs/km; L is the length (km) of the pipeline.   As shown in Equation (14), the key to using an empirical method to solve pipeline failure probability lies in the determination of the RR. In this study, the RR model proposed by Isoyama et al. is adopted [54], which relates the RR to the peak ground acceleration (PGA) by multiplying some correction factors for the pipeline characteristics, surrounding soil characteristics and topology (see Table 1), as shown in the following equations: R 0 = 2.88 × 10 −6 × (PGA − 100) 1.97 (16) where C d , C p , C g , and C l represent the correction factors for the pipe diameter, pipe material, topography, and liquefaction, respectively. R 0 is the standard pipeline RR.  Based on Equations (14)- (16), the failure probabilities of the pipelines under the action of a 9-degree earthquake (0.4 g) are shown in Figure 6. Then the parallel algorithm proposed in Section 4 is used to analyse the SCR of the WDS and compare the results with the traditional serial simulation method based on the Monte Carlo method. It should be noted that both the traditional serial algorithm and the parallel algorithm in this article are implemented in the Visual Studio platform using the C++ programming language. The parallel algorithm is implemented based on the Nvidia CUDA5.0 development kit, and the operating environment is Windows 7. A 2.20-GHz Inter(R)Xeon(R)E5-2407 processor with 4 GB memory and an NVIDIA GeForce GT 610 with 2 GB graphics memory are employed. To verify the accuracy of the parallel method, the serial program and parallel program developed by the author were used to perform 1,000,000 simulations to calculate the seismic connection probability between each node of the pipeline network and the water source under the action of a 9-degree earthquake. As shown in Figure 7, the calculation results of the two methods are basically consistent, and the maximum error is 0.2%, which shows that the parallel method in this paper is correct and meets the accuracy requirements.   Based on Equations (14)- (16), the failure probabilities of the pipelines under the action of a 9-degree earthquake (0.4 g) are shown in Figure 6. Then the parallel algorithm proposed in Section 4 is used to analyse the SCR of the WDS and compare the results with the traditional serial simulation method based on the Monte Carlo method. It should be noted that both the traditional serial algorithm and the parallel algorithm in this article are implemented in the Visual Studio platform using the C++ programming language. The parallel algorithm is implemented based on the Nvidia CUDA5.0 development kit, and the operating environment is Windows 7. A 2.20-GHz Inter(R)Xeon(R)E5-2407 processor with 4 GB memory and an NVIDIA GeForce GT 610 with 2 GB graphics memory are employed. To verify the accuracy of the parallel method, the serial program and parallel program developed by the author were used to perform 1,000,000 simulations to calculate the seismic connection probability between each node of the pipeline network and the water source under the action of a 9-degree earthquake. As shown in Figure 7, the calculation results of the two methods are basically consistent, and the maximum error is 0.2%, which shows that the parallel method in this paper is correct and meets the accuracy requirements.
ployed. To verify the accuracy of the parallel method, the serial program and parallel program developed by the author were used to perform 1,000,000 simulations to calculate the seismic connection probability between each node of the pipeline network and the water source under the action of a 9-degree earthquake. As shown in Figure 7, the calculation results of the two methods are basically consistent, and the maximum error is 0.2%, which shows that the parallel method in this paper is correct and meets the accuracy requirements. To test the computational efficiency of the parallel method in this study, the serial program and the parallel program are used to simulate 3000 times, 5000 times, 10,000 times, 100,000 times, and 1,000,000 times. The simulation times and speedup ratios of the two methods are shown in Table 2 and Figure 8, respectively. It can be seen that with the increase in the number of simulations, the program calculation time and speedup ratio of the two methods increase accordingly. When simulating 1,000,000 times, the parallel method can obtain a 6-fold speedup over the serial method, and shorten the calculation time by 11 s, which shows that the parallel method is helpful to improve the efficiency of the SCR analysis of the pipeline network.  To test the computational efficiency of the parallel method in this study, the serial program and the parallel program are used to simulate 3000 times, 5000 times, 10,000 times, 100,000 times, and 1,000,000 times. The simulation times and speedup ratios of the two methods are shown in Table 2 and Figure 8, respectively. It can be seen that with the increase in the number of simulations, the program calculation time and speedup ratio of the two methods increase accordingly. When simulating 1,000,000 times, the parallel method can obtain a 6-fold speedup over the serial method, and shorten the calculation time by 11 s, which shows that the parallel method is helpful to improve the efficiency of the SCR analysis of the pipeline network.  To verify the computational efficiency of the MC and QMC method, the pseudo-random number generator and the Sobol quasi-random number generator are used to generate random number sequences when using parallel method to analyse the SCR of WDS. The MC and QMC methods were used for 10,000, 20,000, and 30,000 simulations, and the SCR analysis result is shown in Figure 9. It needs to be stated that the connection reliability result when the serial program simulates 1,000,000 times is taken as the analytical value, and the influence of the two random number sequences on the error of the calculation result is analysed. As the number of simulations increases, the simulation accuracy of the two methods is improved. When the number of simulations is the same, the simulation result based on the Sobol sequence is more accurate than that based on the pseudo-random number sequence. When simulating 30,000 times, the errors of the parallel method and the parallel method are 0.82% and 1.4%, respectively. When simulating 10,000 times based on the Sobol sequence, the simulation error is 0.95%, which is smaller than the error of simulating 300,000 times based on the pseudo-random number sequence. In general, both methods can meet the simulation accuracy requirements when simulating 10,000 times, while the Sobol sequence can obtain higher precision results with fewer simulations. To verify the computational efficiency of the MC and QMC method, the pseudorandom number generator and the Sobol quasi-random number generator are used to generate random number sequences when using parallel method to analyse the SCR of WDS. The MC and QMC methods were used for 10,000, 20,000, and 30,000 simulations, and the SCR analysis result is shown in Figure 9. It needs to be stated that the connection reliability result when the serial program simulates 1,000,000 times is taken as the analytical value, and the influence of the two random number sequences on the error of the calculation result is analysed. As the number of simulations increases, the simulation accuracy of the two methods is improved. When the number of simulations is the same, the simulation result based on the Sobol sequence is more accurate than that based on the pseudo-random number sequence. When simulating 30,000 times, the errors of the parallel method and the parallel method are 0.82% and 1.4%, respectively. When simulating 10,000 times based on the Sobol sequence, the simulation error is 0.95%, which is smaller than the error of simulating 300,000 times based on the pseudo-random number sequence. In general, both methods can meet the simulation accuracy requirements when simulating 10,000 times, while the Sobol sequence can obtain higher precision results with fewer simulations.
To verify the applicability of the parallel method in the SCR analysis of a large-scale pipe network, the pipe network in Figure 4 is spliced, and the scale of the pipe networks after splicing is shown in Table 3. Serial and parallel programs are used to analyse the SCR of the spliced pipe network. The simulation time of the two methods when simulating 10,000 times are shown in Table 3. It can be seen from Table 3 that when analysing different scale pipe networks, the parallel algorithm can achieve effective acceleration compared with the serial method. When the number of nodes is 1001, a 4-fold speedup ratio can be obtained.  1  40  64  58  180  5  201  325  250  590  10  401  650  500  1089  15  601  975  745  1568  20  801  1300  1032  2567  25  1001  1625  1165  4524 Water 2023, 15, x FOR PEER REVIEW 15 of 18 To verify the applicability of the parallel method in the SCR analysis of a large-scale pipe network, the pipe network in Figure 4 is spliced, and the scale of the pipe networks after splicing is shown in Table 3. Serial and parallel programs are used to analyse the SCR of the spliced pipe network. The simulation time of the two methods when simulating 10,000 times are shown in Table 3. It can be seen from Table 3 that when analysing different scale pipe networks, the parallel algorithm can achieve effective acceleration compared with the serial method. When the number of nodes is 1001, a 4-fold speedup ratio can be obtained. Table 3. Comparison of the time consumption of the two methods to calculate the SCR of different scale pipe networks.  1  40  64  58  180  5  201  325  250  590  10  401  650  500  1089  15  601  975  745  1568  20  801  1300  1032  2567  25 1001 1625 1165 4524

Conclusions
With the continuous expansion of the WDS scale, it is of great significance to accurately and quickly predict and evaluate the SCR of pipe networks under earthquake action. In this paper, we leverage the immense computational power of GPU to devise an efficient algorithm for SCR analysis of WDS, enabling rapid assessment of a WDS's functional capacity at first glance.
Based on graph theory, the WDS is abstracted into a directed edge-weight graph, and combined with the BFS algorithm, the QMC method of replacing the traditional pseudorandom number sequence with the Sobol sequence is applied to the SCR analysis of the WDS. Through example analysis, it is shown that the QMC method has higher computational efficiency and accuracy than the traditional MC method.
The parallelism of the SCR analysis of the WDS is analysed, and a parallel algorithm for SCR analysis based on the QMC method is proposed and programmed by combining

Conclusions
With the continuous expansion of the WDS scale, it is of great significance to accurately and quickly predict and evaluate the SCR of pipe networks under earthquake action. In this paper, we leverage the immense computational power of GPU to devise an efficient algorithm for SCR analysis of WDS, enabling rapid assessment of a WDS's functional capacity at first glance.
Based on graph theory, the WDS is abstracted into a directed edge-weight graph, and combined with the BFS algorithm, the QMC method of replacing the traditional pseudo-random number sequence with the Sobol sequence is applied to the SCR analysis of the WDS. Through example analysis, it is shown that the QMC method has higher computational efficiency and accuracy than the traditional MC method.
The parallelism of the SCR analysis of the WDS is analysed, and a parallel algorithm for SCR analysis based on the QMC method is proposed and programmed by combining the CUDA platform. Example analysis shows that the results calculated by the parallel algorithm proposed in this paper and the traditional serial method are almost consistent, which shows that the parallel method proposed is correct and meets the accuracy requirements and verifies the rationality of the parallel algorithm. When the scale of the pipe network remains unchanged, with the increase in the number of simulations, the acceleration effect of the parallel algorithm is more obvious than that of the serial method. When the number of simulations remains unchanged, the parallel algorithm has more advantages as the scale of the pipe network increases.
The parallel method proposed in this article is versatile and can also be applied to the seismic reliability analysis of other lifeline systems (gas supply system, road and bridge system, power system, etc.). Data Availability Statement: Data will be made available on request.

Conflicts of Interest:
The authors declare no conflict of interest.