RIVER RESEARCH AND APPLICATIONS
River Res. Applic. (2010) Published online in Wiley Online Library (wileyonlinelibrary.com) DOI: 10.1002/rra.1471
A DEM-BASED PARALLEL COMPUTING HYDRODYNAMIC AND TRANSPORT MODEL YAN JUN ZHANG,a,b* MANOJ JHA,b ROY GU,b LUO WENSHENG a and LEI ALIN c b
a College of Civil Engineering, Wuhan University, Wuhan, China Department of Civil, Construction, and Environmental Engineering, Iowa State University, Ames, Iowa, USA c Changjiang Water Resource Protection Institute, Wuhan, China
ABSTRACT
The sudden and accidental water pollution response system (SAWPRS) for Yangtze River in central China required to develop a hydrodynamic and transport model, which is readily available and capable of simulating a large river system within GIS environment. This study facilitates such effort by developing a parallel computing method based on digital elevation model (DEM) using overlapping domain decomposition approach (ODDA) and message passing interface (MPI) protocol. The hydrodynamic and transport model was redesigned using finite volume method for hydrodynamic and transport model dispersion, the SIMPLEC method for solving the flow field, and the pressure weighted interpolating method for the flow field modification. This modelling approach was verified in two experiments using different sets of computer clusters. The model output was evaluated against the measured data collected for the year 1998 for Wanzhou, an upstream river segment of Yangtze River. The relative error was found to be less than 10%. The performance of parallel computation was found excellent as evident from the cost efficiency values greater than 0.81 in both experiments and increased computation speed while increasing the number of computer clusters. Overall, the parallel computing modelling system developed here was found to meet all requirements of SAWPRS. Copyright # 2010 John Wiley & Sons, Ltd. key words: DEM; MPI; hydrodynamic and transport model; parallel computing; sudden and accidental pollution Received 17 January 2010; Revised 18 August 2010; Accepted 23 September 2010
INTRODUCTION Recently in 2007, the Yangtze water resources committee of China (YWRC) has planned to develop a sudden and accidental water pollution response system (SAWPRS) for Yangtze River that can efficiently deal with sudden contamination of chemicals due to accidents. This action was caused by several incidents that took place in recent years. The Jillin chemical plant was exploded on 13 November 2005 (http://en.wikipedia.org/wiki/Jilin_chemical_plant_explosions_2005) and almost a 1000 km of Songhua River was polluted with an estimated 100 tons of pollutants containing benzene and nitrobenzene. Just 1 month later, on 15 December 2005, Chromium polluted several hundred kilometres of Beijiang River in Guangdong Province (http://news.sina.com.cn/z/gdbjwr/index.shtml). Because of the nature of sudden pollution of rivers and its possible consequences, the YWRC required the modelling system to be capable of simulating a large river system (e.g. over 1000 km river length) with a GIS user friendly interface for efficient visualization and input/output presentation, and readily available. *Correspondence to: Yan Jun Zhang, College of Civil Engineering, Wuhan University, Wuhan 430072, China. E-mail:
[email protected]
Copyright # 2010 John Wiley & Sons, Ltd.
There are several studies done in the area of parallel computation for numerical modelling. Paglieri et al. (1997) developed a parallel shallow water flow model using finite element method, and concluded that the domain decomposition approach was a valid parallel computation method for fine computational meshes and short simulation time. Noel et al. (2000) integrated and parallelized the codes of the CE-QUAL-ICM family of the three-dimensional water quality models for accurate analyses of water quality. The results validated that the finer grid resolution lead to better performance in parallel computation. Similar validation studies were done in various research areas including open channel flow (Rao, 2004), flow field calculation in Wu River (Chunbo et al., 2005), transient shallow water in the presence of wetting/drying fronts (Latorre et al., 2007), groundwater flow and transport simulation (Huang et al., 2008) and sediment suspension methods (Charles et al., 2008). All of these studies simply adopted MPI and domain decomposition approach for parallel computation to decrease the execution time for a fixed computation load or domain size, but none of these focused on redesigning the model for extending the domain size while maintaining the same execution time. Moreover, these studies followed unstructured triangular or body-fitted grid, which makes
Y. J. ZHANG ET AL.
it difficult to use overlapping domain decomposition approach when the domain size is extended. Another requirement of SAWPRS to hydrodynamic and transport modelling is the use of readily available GIS data. Several hydrodynamic and water quality models have been using GIS tools in the pre-processing and post-processing of their modelling efforts, but none have used the GIS data (such as DEM) directly to modify the model configuration. For example, Pinho et al. (2004) created a hydroinformatic environment for solving environmental problems in coastal water; input data and results analyses were performed using the SMS software and ArcView; AUTOCAD was used to edit and digitize the graphic images; and the TRIANGLE software was used to generate the mesh. Also, Ng et al. (2009) presented an integration of GIS and three-dimensional hydrodynamic sediment and heavy metal transport model. Unlike hydrodynamic models, most distributed hydrological and non-point source pollution model uses DEM to design the model configuration. For example, the Soil and Water Assessment Tool model (Arnold et al., 1998) uses DEM to delineate subwatershed boundaries for modelling applications. This study attempts to address above explained two issues of simulation domain extension and direct use of GIS data in hydrodynamic and transport model. A parallel computing method based on DEM was developed using overlapping domain decomposition approach and MPI protocol. DEM consists of Cartesian grid, which is different than normally used body-fitted grid and unstructured triangular grid. It was used as a terrain data. The DEM grid resolution considered as a model grid of the simulation domain. Some preprocessing of DEM including water area recognition and generating water area boundary were added into the model. The overlapping domain decomposition approach was adopted to divide the whole simulation domain into several sub domains for parallel computing. MPI was applied to exchange the data between sub domains and build the parallel computation frame for the hydrodynamic and transport model. Furthermore, the model was redesigned (to incorporate grid approach) using finite volume method (FVM) (Patankar, 1980) for hydrodynamic and transport model dispersion, the SIMPLEC method (Van Doormal and Raithby, 1984) for solving the flow field, and the pressure weighted interpolating method (PWIM) (Miller and Schmidt, 1988) for the flow field modification. Figure 1 presents the flowchart of the modelling scheme. This modelling approach was verified in two experiments using different sets of computer clusters which are much cheaper than the super computer for Wanzhou, an upstream river segment of the Yangtze River. In brief, the objectives of this study were: (1) to use/integrate DEM data (Cartesian grid) into the parallel computing hydrodynamic and transport model; Copyright # 2010 John Wiley & Sons, Ltd.
Figure 1. Flow chart of model implements.
(2) to be able expand the simulation domain to a very large area in parallel computing environment and (3) to verify the developed DEM-based parallel computing algorithm in the Yangtz River system.
DEVELOPMENT OF THE ALGORITHM Parallel computing method Division. DEM has Cartesian grid format, which consists of many small rectangular grids with elevation information (Figure 2). This feature of DEM is advantageous in extending or reducing the size of the simulation domain
Figure 2. The 3D effect drawing of a DEM and its Cartesian grid. This figure is available in colour online at wileyonlinelibrary.com
River Res. Applic. (2010) DOI: 10.1002/rra
PARALLEL COMPUTING HYDRODYNAMIC AND TRANSPORT MODEL
as well as when combining with other DEMs as long as they have the same precision and the same reference point. This also makes it easy for overlapping domain decomposition method to divide the grids for parallel computation. The X and Y coordinates of DEM, and a time variable, T, form a three-dimensional domain. Ranges of each grid point coordinates are: X ¼ 1, 2, 3 . . . I; Y ¼ 1, 2, 3 . . . J and T ¼ 1, 2, 3 . . . N. These points are regarded as initial tasks with assigned variables: u (X-dimensional velocity), v (Ydimensional velocity), z (water stage) and c (concentration of water quality constituents such as DO, COD, etc.). Communication. Finite Volume Method (FVM) (Patankar, 1980) was used to solve hydrodynamic and transport model. This method uses five input channels and five output channels for each initial task. For example, ’iþ1;j;t , ’i;jþ1;t , ’i1;j;t , ’i;j1;t , ’i;j;t1 are needed to calculate ’i;j;t . Figure 3 shows the example communication mode used in the FVM method.
Congregating and mapping. Because of the large number of the initial tasks, it is impossible to calculate all of them at the same time, and so we chose to congregate initial tasks. Figure 4 depicts the congregation process. This was done in three steps: firstly, all the initial tasks which have same X, Y coordinates but different T coordinates merge; secondly, all the first step tasks which have same X coordinates but different Y coordinates merge and finally, a group of adjacent second step tasks merge to form a final task using the load balancing principle. Every final task is assigned a core (central processing unit) to calculate variables x, y, z and c. Every core calculates independently and communicates with its neighbour. This enables the use of same computer code to execute in all cores independently. Congregating procedure creates two new boundaries called exchange-in and exchange-out boundaries between cores, as shown in Figure 5. The u, v, z and c variables of one core, at the exchange-in boundary must be received from the exchange-out boundary of neighbouring core and vice versa. Performance of parallel computation. This study used two indices to measure the performance of the parallel computation algorithm: Speedup (SP) and cost efficiency (EP). These can be defined as: Sp ¼ TT1p
for fixed computation load; and
s Sp ¼ TTPs =D =DP
) (1)
for changed computation load
and, T1 EP ¼ pT p
for fixed computation load; and
sðTs =Ds Þ EP ¼ pðT p =Dp Þ
Figure 3. Communication mode used in the FVM method for initial tasks. This figure is available in colour online at wileyonlinelibrary.com
)
for changed computation load
(2)
where T and D represent the execution time and coverage length, T1 represents execution time of the sequential algorithm, TP represents execution time for P number of cores and subscript s represents the minimum number of cores involved in the computation.
Figure 4. Schematic diagram of three step congregation process of initial tasks. This figure is available in colour online at wileyonlinelibrary.com Copyright # 2010 John Wiley & Sons, Ltd.
River Res. Applic. (2010) DOI: 10.1002/rra
Y. J. ZHANG ET AL.
Two additional boundary types were generated after the simulation domain was divided into sub-domains: (1) exchange-out node (BS ¼ 5): nodes whose variables (u, v, z and c) are sent out to the neighbouring domain and (2) exchange-in node (BS ¼ 6): nodes whose variables (u, v, z and c) are received from the neighbouring domain. The precision of the DEM grid was set to less than 1/3 width of the most narrow water area for the computational accuracy of the result. This precision assured at least one normal node around the boundary node. The equations and solution
Figure 5. Schematic diagram showing the congregation and communication process between the cores in DEM grid (BS: Boundary Status). This figure is available in colour online at wileyonlinelibrary.com
Processing DEM DEM contains elevation and coordinates information but lacks other information such as water area and their boundary. In contrast, other grids such as unstructured grid and body-fitted grid, commonly used in hydrodynamic and transport model, water area and boundary are known but the grid size is unknown. These are often difficult to generate accurately. The submerge analysis method (Lihua and Yi, 2002) was used to define water area in the DEM grid using three main steps: (1) searching hollow, (2) delete pseudo hollow and (3) non-source flooding analysis. After recognizing water area, removing skin transformation (Moulting method) (Peng et al., 2006) was used to generate the water area boundary. Then, five types of nodes were considered for analysis (also shown in Figure 5): (1) invalid node (boundary status, BS ¼ 0): These nodes are located on the land and ignored in the calculation, (2) boundary node (BS ¼ 1): These are calculated by the boundary processing approach, (3) normal node (BS ¼ 2): These are covered by water, (4) inflow node (BS ¼ 3): Water flow into the area through these nodes. At this node, variables u, v and c are given and z is calculated. (5) outflow node (BS ¼ 4): Water flow out of the area through these nodes. At this node, variable z is given and u, v and c are calculated.
Copyright # 2010 John Wiley & Sons, Ltd.
Governing equations. We used two-dimensional (2D) shallow stream hydrodynamic and transport model, as follows (the shearing stress of wind and Coriolis force have been neglected): ~ ~ @G @U @F @G @F þ þ ¼ þ þb @t @x @y @x @y
(3)
2
3 2 3 2 3 h hu hv 6 hu 7 6 hu2 þ ghz 7 6 huv 7 7 6 7 6 7~ U¼6 4 hv 5 F ¼ 4 huv 5 G ¼ 4 hu2 þ ghz 5 F hci huci hvci 2 3 2 3 0 0 @u " h 6 "x h @u 7 6 x @x 7 @x 7 ~ 6 7b ¼ F ~ ¼6 G ¼ @v 4 "y h @y 5 4 "y h @v @y 5 i i Ey @hc Ex @hc @y @x 2 3 qpffiffiffiffiffiffiffiffiffi 2 þv2 6 gn2 u u1=3 7 phffiffiffiffiffiffiffiffiffi 7 ¼6 (4) 2 4 gn2 v u þv2 5 h1=3 hðSfi Ci þ S0i where X and Y represent vertical and horizontal length of water area; u and v represent vertical and horizontal velocity of the water; t is time; h is water depth; z is water level; hb is elevation of the surface below water; ci is concentration of the i water quality index (Carbonaceous Oxygen Demand (COD) or Dissolved Oxygen (DO)); "x and "y are vertical and horizontal eddy viscosity coefficients; g is the gravitational constant; Ex and Ey are the sum of vertical molecular diffusion coefficient, turbulent diffusion coefficient and dispersion coefficient (calculated by Jingxiu’s model (Jingxiu and Zhengli, 2000)); n is roughness coefficient; Q is interzone inflow quantity; Sfi is the first-order derivative of source and sink terms of water quality index i and Soi is constant for the source and sink terms of water quality index i.
River Res. Applic. (2010) DOI: 10.1002/rra
PARALLEL COMPUTING HYDRODYNAMIC AND TRANSPORT MODEL
where Qa is temperature coefficient (1.08) and v is average water velocity (m/s). The DO saturation (Cs) can be calculated as follow: lnCs ¼ 139:34 þ 1:5757 105 TK1 6:6423 107 TK2 þ 1:2438 1010 TK3 8:6219 1011 TK4 0:5535Sð0:03192919:428TK1 þ 3867:3TK2 Þ (8) where Tk is temperature (degree K) and S is salinity (mg/L). From Equations (5) and (6), following equations can be derived, to be used in Equations (3) and (4) above: Sf1 ¼ kd QT20 d
vs3 ð1fD5 Þ D
Sf2 ¼ k2 QT20 2
(10)
S01 ¼ 0
(11)
S02 ¼ k2 QT20 cS kd QT20 c2 2 d
Figure 6. DEM of Wanzhou section of the Yangtze River.
The COD and DO have been chosen to represent for the water quality indices. Their inter-relationship is explained by the Streeter–Phelps equations in a slightly modified form, also adopted in an in-stream water quality model WASP (Water quality Analysis Simulation Program) (http:// www.epa.gov/athens/wwqtsc/html/wasp.html): Sk1 ¼ kd QT20 c1 d
vs3 ð1fD5 Þc1 D
(5)
SODT D
Sk2 ¼
D
(6)
where Sk1 is source and sink terms of COD; Sk2 is source and sink terms of DO; kd is deoxygenation rate at 208C (0.210.16 day1); k2 is re-aeration rate; Cs is saturation of DO (mg O2/L); SOD is sediment oxygen demand (0.20.4 g/m2-day); D is average segment depth (m); Qd is temperature coefficient (1.047); Q2 is temperature coefficient (1.028); vs3 is organic matter settling velocity (m/day) and fD5 is fraction dissolved COD (0.5). Re-aeration rate (k2) is calculated using following equation: k2 ¼ 5:049v0:97 D1:67 Qa T20 Copyright # 2010 John Wiley & Sons, Ltd.
(7)
(12)
where Sf1is first-order derivate of source and sink terms of COD; Sf2 is first-order derivate of source and sink terms of DO; S01 is constant of source and sink terms of COD; S02 is constant of source and sink terms of DO and vs3 is organic matter settling velocity (m/day). Discretization and solution. The FVM method (Patankar, 1980) was used to discretize Equations (1) and (2). It is an implicit scheme which brings more calculation load, but the robustness and precision, which FVM brings, are necessary. The SIMPLEC method (Van Doormal and Raithby, 1984) was used to solve flow fields u, v and c. Because of the square DEM grid, most flow field modification method is useless for FVM. The PWIM method (Miller and Schmidt, 1988) was used to solve h in the collocated grid. Finally, the following equation can be formed: aP fP ¼ aE fE þ aW fW þ aN fN þ aS fS þ b
SODT k2 QT20 ðCS c2 Þkd QT20 c2 2 d
(9)
(13)
where a and b are coefficients and f represents u, v, z and c variables at a point (P) in east (E), west (W), north (N) and south (S) directions. These variables are derived in details and presented in Appendix A. Equation (13) was solved for each normal node using two methods: Tri-Diagonal Matrix Algorithm method (TDMA) and Alternating Direction Implicit method (ADI). The details on these methods can be found in Wenquan (2001). Initial and boundary conditions. For boundary nodes (status 1 nodes), vertical velocity is zero and the tangential velocity follow f ¼ af0, where f represents u, v and c 0 variables in boundary nodes, f represents the variables in its closest normal node (status 2 nodes), a is a coefficient (0–1, 0.8 in this paper) For normal nodes (status 2 nodes), River Res. Applic. (2010) DOI: 10.1002/rra
Y. J. ZHANG ET AL.
Figure 7. The domain division with (a) one core, (b) two cores, (c) three cores and (d) five cores.
Copyright # 2010 John Wiley & Sons, Ltd.
River Res. Applic. (2010) DOI: 10.1002/rra
PARALLEL COMPUTING HYDRODYNAMIC AND TRANSPORT MODEL
variables u and v are zero, c is initial concentration and z is obtained from interpolation. Moving boundary problem was solved using dry–wet method, as explained in Xia et al. (2006). In this method, if water depth (h) is less than zero in normal nodes, its status changed to invalid nodes and water drained into the neighbouring nodes. Similarly, if the water stage (z) is more than the elevation of node (hb) in invalid status, its status changed to normal nodes and the water drew-in from the neighbouring nodes. Status of all nodes is refreshed after applying dry–wet method in each iteration step.
COMPUTING SPEED OPTIMIZATION Sparse matrix compression It is important and helpful to compress the DEM grid (Sparse matrix compression) to reduce the calculation load as well as the storing space. A large number of invalid points (entries) whose boundary status is 0 will be generated after recognizing the water area in the domain grid (twodimensional sparse matrix, M). The compression was done by storing sparse metrics only with non-zero entries. This study used three one-dimension arrays to store a sparse m n matrix. The first array is A, which holds all non-zero entries (denoted by NNZ) of M in left-to-right top-to-bottom order. The second array is IA, which contains the row index of each element of the first array A. The third array is JA, which contains the column index of each element of the first array A. This sparse matrix compression decreased the calculation load significantly by removing null entries. Also, the storing space increased significantly as the total number of entries reduced. Iterations are normally set to stop if the number of iterations is more than a specified number or when the relative or the absolute error is less than the specified value. In a parallel computing environment, computation time is different for different parts due to the fact that the external disturbance has variable effect on each part, and so the absolute or relative error is considered to be the best choice for iteration termination.
Table I. The monitoring data of hydrology and water quality at Tuo Kou on the Yangtze River Date (1998)
Time
Water stage (m)
Discharge (m3/s)
COD (mg/L)
15 March
8 AM 2 PM 8 AM 2 PM
100.74 100.74 100.77 100.77
4150 — — 4330
2.14 — — 2.12
16 March
Source: Water quality monitoring centre of YWRC.
Administration of China (http://www.zhb.gov.cn/tech/hjbz/ bzwb/shjbh/shjzlbz/200206/t20020601_66497.htm). The model requires the terrain data (DEM), hydrology data and water quality data. These were collected from the water quality monitoring centre of YWRC. The solution of the DEM data was 20 m 20 m and the size of grid was 413 531. The DEM grid size was used directly as the grid of the simulation domain. The monitoring data of hydrology and COD at Tuo Kou was collected for 16 and 17 March 1998, and is presented in Table I. Table II lists all the effluent sources (also indicated in Figure 6) and the related hydrology and water quality data. The COD concentration data at Sai Wangba was collected to compare the model calculation (presented in Table III). Two experiments A and B were designed with two different sets of computer clusters. This was considered a viable option for developing countries like China where computer clusters are much cheaper than the super computer. In experiment A, 12.4 km of the river length (simulation domain) between Tou kou and Sai Wangba was divided into one (Figure 7a), two (Figure 7b) and three (Figure 7c) subdomain. Then, the simulated results of COD concentrations were compared with the measured data to check the accuracy and reliability of the model. Table II. Industrial effluent in the modelling domain on Yangtze River Pollution source MoDao Brook
Data and experimental program. The DEM-based parallel computing hydrodynamic and transport model, developed in this study, was applied to the Wanzhou section of the Yangtze River, which includes the area between Tou Kou and Sai Wangba as shown in Figure 6. The regions of Tou Kou and Sai Wangba are classified as CLASS 3 (COD level 4–6 mg/L) by the State Environmental Protection Copyright # 2010 John Wiley & Sons, Ltd.
Time
Discharge (m3/s)
COD (mg/L)
15 March
9:00 11:00 14:00 17:00 9:00 11:00 14:00 17:00 9:00 11:00 14:00 17:00
0.070
3108 3028 1295 1444 757 817 667 647 139 120 79.7 120
16 March
MODEL APPLICATION Model setup
Date (1998)
WanYuan Paper Mill
15 March 16 March
Zuxi River
15 March 16 March
0.290
0.085
Source: Water quality monitoring centre of YWRC. River Res. Applic. (2010) DOI: 10.1002/rra
Y. J. ZHANG ET AL.
Table III. The comparison of calculated and measured COD concentration at Sai Wangba Date and Time
Distance of the measured point to the right bank (m)
16 March 1998, 3
PM
17 March 1998, 9
AM
31 71 141 201 236 31 71 141 201 236
Additionally, the computation time among the subdomains were compared to check the performance of parallel computing. In experiment B, three computation areas were formed, one with three subdomains (Figure 7c), another with five subdomains (Figure 7d) and the last one with seven subdomains (Figure 7e). The first computation domain has three subdomains with 12.4 km of river length (same as experiment A), and the second computation domain has five subdomains with 20.2 km of river length (with added 7.8 km of the new domain area). Again, the computation times are compared to check the model’s ability to increase coverage area as the number of sub domains increases. Experiment A has a setup of three cores using two different PCs, one Intel P8400 processor (2 cores) and one Celeron M370 processor (1 core). These processors exchange message by a 100M Ethernet hub. Similarly, experiment B has a setup of five cores including three cores of experiment A and four cores of other two Intel P8400 processor. Also, the cores could be added according to the computation load needs. Parameters calibration. The major model parameters were calibrated using methods presented in Zhenli and Yuliang (2007). The final calibrated values for the roughness coefficient (n) was found to be 0.04, the vertical and horizontal eddy viscosity coefficient was found to be 22.6 Pa s, and the degradation coefficient (Kd) was found to be 0.23 (Table III). The maximum relative error was found to be less than 10%.
COD Concentration (mg/L)
Relative error (%)
Measured
Calculated
2.9 2.8 2.7 2.6 2.4 2.18 2.17 2.16 2.15 2.14
2.8 2.6 2.5 2.6 2.5 2.2 2.2 2.4 2.2 2.2
3.6 7.7 8.0 0.0 4.0 0.9 1.4 10.0 2.3 2.7
with enough confidence as the values of residuals converged to lower than 106. Table III shows the comparison between model calculated COD concentration with the measured data at Sai Wangba. The calculated values follow closely with the measured data. The relative error was found to be less than 10% in all cases with majority less than 4%. It showed that the DEM-based parallel computation hydrodynamic and transport model developed in this study was able to simulate COD concentration very well. The model results of flow velocity and COD concentration are plotted in Figures 9 and 10. The flow velocity found very high in the range of 1–2 m/s as shown in Figure 9. The velocity near the bank is slower than at the midstream. Figure 10 shows the COD concentration extended well beyond the CLASS 3 standard (4–6 mg/L) of the State Environmental Protection Administration of China, with length approaching nearly 8 km at MoDao Brook.
Results Model verification. The model was verified using the experimental setup A with three cores. Figure 8 shows the mode residuals plotted against the number of iterations. As the number of iterations increased, residuals decreased rapidly and then it changed to a gradual decrease after 300 iterations. The program was terminated after 1000 iterations Copyright # 2010 John Wiley & Sons, Ltd.
Figure 8. Model residuals against the number of iterations.
River Res. Applic. (2010) DOI: 10.1002/rra
PARALLEL COMPUTING HYDRODYNAMIC AND TRANSPORT MODEL
Figure 9. Flow velocity field at MoDao Brook and Paper Mill effluent sources. This figure is available in colour online at wileyonlinelibrary.com
Figure 10. Areal extent of COD concentration at MoDao Brook and Paper Mill effluent sources. This figure is available in colour online at wileyonlinelibrary.com
Parallel computation efficiency. Table IV shows the performance efficiencies of parallel computation algorithm used in the experiments A and B. For experiment A with fixed computation load (simulating 12.4 km of river length), Sp and EP was found to be 1.7 and 0.83 when two cores were used and 2.5 and 0.82 when three cores were used. The comparison of performance of above two sets of clusters validates that the parallel computation increase the speed of the computation while maintaining reasonable cost efficiency. In experiment B, computation load was increased from 12.4 to 20.2 km of river length. It can be seen from Table IV that the three cores simulated 12.4 km of river in 125 s, while five cores simulated 20.2 km of river in similar time, 127 s. The Ep and Sp were found to be 0.96 and 1.6, which explains that the computing speed with five cores increased tremendously while comparing with the three core
performances. This experiment validates that the parallel computing algorithm can simulate larger coverage area while maintaining the reasonable computation time or reduces the computation time significantly for the same coverage area. CONCLUSIONS AND DISCUSSION Recently in 2007, the Yangtze water resources committee of China (YWRC) has planned to develop a sudden and accidental water pollution response system (SAWPRS) for Yangtze River that can efficiently deal with sudden contamination of chemicals due to accidents. The modelling system should be capable of simulating a large river system (e.g. over 1000 km river length) with a GIS user-friendly interface for efficient visualization and input/output presentation, and readily available. A parallel computing method
Table IV. The parallel computation efficiency of computer cluster for experiment A with fixed coverage area and for experiment B with increased coverage area Number of cores (P) Experiment A Experiment B
One Two Three Three Five Seven
Copyright # 2010 John Wiley & Sons, Ltd.
Scheme
Time for 1000 iteration (s)
River length (km)
Speedup (SP)
Cost efficiency (EP)
Figure 7a Figure 7b Figure 7c Figure 7c Figure 7d Figure 7e
306 183 125 125 127 129
12.4 12.4 12.4 12.4 20.2 27.8
1.0 1.7 2.5 1.0 1.6 2.2
1.00 0.83 0.82 1.00 0.96 0.94
River Res. Applic. (2010) DOI: 10.1002/rra
Y. J. ZHANG ET AL.
Figure 11. Schematic diagram of TGWPERS running. This figure is available in colour online at wileyonlinelibrary.com
based on DEM was developed using overlapping domain decomposition approach and MPI protocol. The implementation of this method into the SAWPRS system is exemplified in Figure 11. The model was redesigned (to incorporate grid approach) using finite volume method for hydrodynamic and transport model dispersion, the SIMPLEC method for solving the flow field, and the pressure weighted interpolating method (PWIM) for the flow field modification. The developed modelling approach was verified in two experiments using two different sets of computer clusters for Wanzhou, an upstream river segment of the Yangtze River. The comparison of the model results with the measurement for the year 1998 yielded less than 10% relative error. The performance of parallel computation was found excellent as evident from the cost efficiency values greater than 0.81 in both experiments and increased computation speed (or simulation domain) while increasing the number of computer clusters. This study also emphasized the use of DEM grid over body-fitted or triangular grids for extending coverage area because it is easy to overlap different DEMs together as long as they have the same precision and same reference point. The developed algorithm automated the pre-processing of the DEM grid such as water area recognition and generating water area boundary. This facilitates users to quickly use the algorithm in the new area if only DEM, hydrology data and water quality data are available. Overall, it can be concluded Copyright # 2010 John Wiley & Sons, Ltd.
that the parallel computational hydrodynamic and transport model developed here can be used to simulate a large domain within a GIS environment.
REFERENCES Arnold JG, Srinivasan R, Muttiah RS, Williams JR. 1998. Large-area hydrologic modeling and assessment: Part I model development. Journal of American Water Resource Association 34(1): 73–89. Charles WM, van den Berg E, Lin HX, Heemink AW, Verlaan1 M. 2008. Parallel and distributed simulation of sediment dynamics in shallow water using particle decomposition approach. Journal of Parallel and Distributed Computing 68(6): 717–728. Chunbo J, Kai, L., Ning, L., Qinghai, Z., 2005. Implicit parallel FEM analysis of shallow water equations. Tsinghua Science and Technology 10(3): 364–371. Huang J, Christ JA, Goltz MN. 2008. An assembly model for simulation of large-scale ground water flow and transport. Ground Water 46(6): 882– 892. Jingxiu L, Zhengli H. 2000. Study on longitudinal dispersion coefficient of the three gorges. Journal of Hydraulic Engineering (8): 84–87. (Chinese). Latorre B, Murillo J, Garcı´a-Navarro P. 2007. Parallel computation of unsteady inundation shallow-water flow. International Workshop on Numerical Modelling of Hydrodynamics for Water Resources. Zaragoza: Spain; June 18-21. Lihua G, Yi L. 2002. Analysis of flood submerging based on DEM. Bulletin of Surveying and Mapping 11: 25–30. (Chinese). Miller TF, Schmidt FW. 1988. Use of pressure weighted interpolation method for the solution of incompressible Navier-Stokes equations on a nonstaggered grid system. Numerical Heat Transfer 14(2): 213–233. River Res. Applic. (2010) DOI: 10.1002/rra
PARALLEL COMPUTING HYDRODYNAMIC AND TRANSPORT MODEL
Ng SMY, Wai OWH, Li. Y, Li. L, Jiang Z, Y., 2009. Integration of a GIS and a complex three-dimensional hydrodynamic, sediment and heavy metal transport numerical model. Advances in Engineering Software 40(6): 391–401. Noel MR, Gerald TK, Cerco CF. 2000. Parallel performance and benchmarking of the CE-QUAL-ICM family of three-dimensional water quality models. Applications of High-Performance Computing in Engineering VI: 301–309. Paglieri L, Ambrosi D, Formaggia L, Quarteroni A, Scheinine AL. 1997. Parallel computation for shallow water flow: a domain decomposition approach. Parallel Computing 23(9): 1261–1277. Patankar SV. 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corp: Washington DC. Peng H, Lian Y, Chuanyong Y. 2006. Map Algebra. Wuhan University Press: Wuhan, China; 137–140. (Chinese).
Pinho JLS, Pereira JM, Vieiraa JS, do Carmob A. 2004. Hydroinformatic environment for coastal waters hydrodynamics and water quality modeling. Advances in Engineering Software 35(3–4): 205–222. Rao P. 2004. A parallel hydrodynamic model for shallow water equations. Applied Mathematics and Computation 150(1): 291–302. Van Doormal JP, Raithby GD. 1984. Enhancement of the SIMPLE method for predicting incompressible fluid flows. Numerical Heat Transfer 7(1): 147–163. Wenquan T. 2001. Numerical Heat Transfer (2nd edn). Xi’an Communication University Press: Xi’an, China; (Chinese). Xia S, Yao Q, Wang P. 2006. Study and application of moving boundary in POM. Marine Science Bulletin 25(1): 1–6. (Chinese). Zhenli H, Yuliang L. 2007. Water Quality Prediction and Water Environmental Carrying Capacity Calculation for Three Georges Reservoir. China Water Power Press: Beijing, China (Chinese).
DISCRETIZATION EQUATIONS Following equations are derived from Equation (13): 0 1 u Bv C C f¼B @ h0 A ci 0 1 ðaE þ aW þ aN þ aS þ aop Sp DXDYÞ=a B ðaE þ aW þ aN þ aS þ aop Sp DXDYÞ=a C C ap ¼ B @ A aE þ aW þ aN þ aS o ðaE þ aW þ aN þ aS þ ap Sp DXDYÞ=a 0 1 De AðjPe jÞ þ ½½Fe ; 0 B D AðjP jÞ þ ½½F ; 0 C e e B e C aE ¼ B C @ A d e he Dy De AðjPe jÞ þ ½½Fe ; 0 0 1 DW AðjPW jÞ þ ½½FW ; 0 B DW AðjPW jÞ þ ½½FW ; 0 C B C aW ¼ B C @ A d w hw Dy DW AðjPW jÞ þ ½½FW ; 0 1 Dn AðjPn jÞ þ ½½Fn ; 0 B D AðjP jÞ þ ½½F ; 0 C n n B n C aN ¼ B C @ A dn hn Dx 0
Dn AðjPn jÞ þ ½½Fn ; 0 0
1 DS AðjPS jÞ þ ½½FS ; 0 B DS AðjPS jÞ þ ½½FS ; 0 C B C as ¼ B C @ A ds hs Dx DS AðjPS jÞ þ ½½FS ; 0 0
Sc DXDY þ fop aop þ ð1aÞ=aap up
1
C B Sc DXDY þ fop aop þ ð1aÞ=aap up C B C b¼B B u he Dy þ u hw Dyv hn Dx þ v hs Dx þ ðh0 h ÞDxDy=Dt þ qDxDy C w n s p p A @ e Sc DXDY þ fop aop þ ð1aÞ=aap up Copyright # 2010 John Wiley & Sons, Ltd.
River Res. Applic. (2010) DOI: 10.1002/rra
Y. J. ZHANG ET AL.
0
1 hop DXDY=Dt B ho DXDY=Dt C B C aop ¼ B p C @ A hop DXDY=Dt 0 1 0 1 0 1 0 1 ðuhÞw Dy ðvhÞs Dx ðuhÞe Dy ðvhÞn Dx B ðuhÞ Dy C B ðvhÞ Dx C B ðuhÞ Dy C B ðvhÞ Dx C B C B C C B C w s e n Fe ¼ B @ A Fw ¼ @ A Fn ¼ @ A Fs ¼ @ A ðuhÞw Dy ðvhÞn Dx 0 1 ð"x hÞe Dy=Dx ð"x hÞw Dy=Dx B ð"x hÞ Dy=Dx C B ð"x hÞ Dy=Dx C C Dw ¼ B C e w De ¼ B @ A @ A 0
ðuhÞe Dy
1
ðvhÞs Dx
ð"x hÞw Dy=Dx 1 ð"y hÞs Dx=Dy ð"y hÞn Dx=Dy B ð"y hÞ Dx=Dy C B C C Ds ¼ B ð"y hÞs Dx=Dy C n Dn ¼ B @ @ A A 0
ð"x hÞe Dy=Dx
1
0
ð"y hÞn Dx=Dy ð"y hÞs Dx=Dy 1 Fe =De B F =D C e eC Pe ¼ B @ A Fe =De 0 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gn2 u2 þ v2 =h1=3 q B pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1=3 C 2 2 2 B C Sp ¼ B gn u þ v =h q C @ A Sfi H 0 1 ghðZe Zw Þ=Dx B ghðZ Z Þ=Dx C e w C Sc ¼ B @ A S0i H 0 1 B C C ue ¼ B @ A P ðð anb unb þ bÞ=aP Þe þ ðghDy=aP Þe ðZP ZE Þ 0 1 B C C de ¼ B @ A ðDy=aP ÞP ðdxÞeþ =ðdxÞe þ ðDy=aP ÞE ðdxÞe =ðdxÞe 0
where a is the relaxation factor of SIMPLEC method, usually, a < 1.
Copyright # 2010 John Wiley & Sons, Ltd.
River Res. Applic. (2010) DOI: 10.1002/rra
