## 1 Introduction  
Hydraulic transients are the fluctuations of pressure and flow in water (oil, gas) transmission or distribution systems caused by pump or valve operations, and their numerical simulations are essential for design, operation and management of engineering projects [1; 2; 3]. The existing simulation methods of hydraulic transients are adequate for many practical applications in terms of accuracy and efficiency [4]. But in some special cases, simulation efficiency must be improved to better meet demands. First, more and more water transmission systems longer than hundreds of kilometers (e.g., the Da Huofang project and the Liaoning northwest water supply project in China are 354 kilometers and 598.4 kilometers, respectively) and water distribution networks containing thousands of pipes, loops,and branches (e.g., the networks in Xiongan New Capital City and Wuhan New Yangtze River City, China) are newly built and planned, for which the hydraulic transient simulations are very time-consuming, normally needing several hours for calculating one working condition. Second, some of these long-distance systems are equipped with real-time monitoring and need online transient simulations for instant decision-making on operations, for which the existing methods are not competent [5]. Third, in design and optimization of many complex pipe systems, a large number of layouts and working conditions should be simulated in a limited time [6]. Therefore, it is necessary to find more efficient simulation methods of hydraulic transients. Two approaches on improving the computing efficiency of hydraulic transients have already been done. The first is to adopt more efficient algorithms. [PERSON] et al. [7; 8] proposed the wave characteristic method (WCM) to efficiently simulate the transient processes of water distribution systems. [PERSON] et al. [9] and [PERSON] et al. [10] further proved WCM is more efficient than the widely-used method of characteristics (MOC) because of its larger time-step interval and exemption from interior section computation. The second approach is to optimize or parallelize computing codes. [PERSON] and [PERSON] [11; 12] established a sub-procedure library for general hydraulic boundaries, which were invoked directly in simulations to avoid tiny segmentation and reduce resource occupancy. [PERSON] et al. [13] used a parallel genetic algorithm (PGA) in the transient simulations with cavitation and obtained about 10 times of speedup compared to the serial genetic algorithm. [PERSON] et al. [14] used an optimized genetic algorithm to redesign existing water distribution networks that do not operate properly. [PERSON] et al. [15] parallelized two computers to simulate the transients of several large hydraulic systems, and achieved an efficiency of 1.442 times to that on a single computer. However, these efficiency improvement are far from meeting the practical needs, therefore, new approaches should be explored. The graphic processing unit (GPU) has shown powerful capability in non-graphical computations and scientists have started using it to conduct large-scale scientific calculations at the beginning of the 21 century [16]. Many successful GPU applications in fluid dynamics calculations have shown quite promising performance gains. [PERSON] et al. [17] accelerated the simulations of fluid structure interaction (FSI) by implementing the immersed boundary-lattice Boltzmann coupling scheme on a single GPU chip, and the attained memory efficiencies for the kernels were up to 61%, 68%, and 66% for the testing case. [PERSON] et al. [18] accelerated the simulations of high-enthalpy flows by implementing the state-to-state approach on a single GPU chip, the comparison between GPU and CPU showed a speedup ratio up to 156, and such speedup ratio increased with the problem size and the problem complexity. [PERSON] et al.