# Heterogeneous Parallel Implementation of Large-Scale Numerical Simulation of Saint-Venant Equations
## 2 Governing Equations and Numerical Method  
The governing equations used in our research are as follows:  
\[\frac{\partial z}{\partial t}+\frac{\partial(zu)}{\partial x}+\frac{\partial( zv)}{\partial y}=0 \tag{1}\]  
\[\frac{\partial(zu)}{\partial t}+\frac{\partial\Big{(}zu^{2}+\frac{gz^{2}}{2} \Big{)}}{\partial x}+\frac{\partial(zuv)}{\partial y}=gz\Big{(}S_{0x}-S_{fx} \Big{)} \tag{2}\]  
\[\frac{\partial(zv)}{\partial t}+\frac{\partial(zuv)}{\partial x}+\frac{ \partial\Big{(}zv^{2}+\frac{gz^{2}}{2}\Big{)}}{\partial y}=gz\Big{(}S_{0x}-S_{ fx}\Big{)} \tag{3}\]  
Equation (1) is derived from the conservation of mass, and Equations (2) and (3) are derived from the conservation of momentum in the \(x\) and \(y\) directions, respectively. Among them, \(z\) refers to the elevation (depth or height) of the water flow in the open channel, \(u\) and \(v\) are the water velocity in the \(x\) and \(y\) directions respectively, \(t\) is the time, g is the acceleration of gravity. \(S_{0}\) and \(S_{f}\) refer to the water surface gradient and frictional resistance. In order to use this system of equations more conveniently in the computer, we need to use the product rule of differentiation to further simplify the system of equations. We then discretize the equations using an explicit finite difference scheme, where the time and space derivatives are respectively discretized by the following expressions:  
\[\frac{\partial u}{\partial t}\approx\frac{u_{i,j}^{k+1}-u_{i,j}^{k-1}}{2 \Delta t},\frac{\partial u}{\partial x}\approx\frac{u_{i+1,j}^{k}-u_{i-1,j}^{ k}}{2\Delta x},\frac{\partial u}{\partial y}\approx\frac{u_{i,j+1}^{k}-u_{i,j-1}^{ k}}{2\Delta y}\]  
This format is the central difference format, in which the central difference in time is also called the leapfrog format, and its advantage is that it can enhance the stability of the calculation. In this way, the discretized equations become the following form:  
\[\frac{z_{i,j}^{k+1}-z_{i,j}^{k-1}}{2\Delta t}+u_{i,j}^{k}\frac{z_{i+1,j}^{k}- z_{i-1,j}^{k}}{2\Delta x}+z_{i,j}^{k}\frac{u_{i+1,j}^{k}-u_{i-1,j}^{k}}{2 \Delta x}+v_{i,j}^{k}\frac{z_{i,j+1}^{k}-z_{i,j-1}^{k}}{2\Delta y}+z_{i,j}^{k} \frac{v_{i,j+1}^{k}-v_{i,j-1}^{k}}{2\Delta y}=0 \tag{4}\]  
\[ \begin{split} u_{i,j}^{k}\frac{z_{i,j}^{k+1}-z_{i,j}^{k-1}}{2 \Delta t}+z_{i,j}^{k}\frac{u_{i,j}^{k+1}-u_{i,j}^{k-1}}{2\Delta t}+\Big{(}u_{ i,j}^{k}\Big{)}\frac{z_{i+1,j}^{k+1}-z_{i,j}^{k-1,j}}{2\Delta x}+2z_{i,j}^{k}u_{i,j}^{k }\frac{u_{i+1,j}^{k+1}-u_{i,j}^{k-1}}{2\Delta x}+gz_{i,j}^{k}\frac{z_{i+1,j}^ {k}-z_{i,j}^{k}-z_{i,j}^{k}}{2\Delta x}+\\ u_{i,j}^{k}\frac{z_{i,j+1}^{k}-z_{i,j-1}^{k}}{2\Delta x}+z_{i,j}^{k} v_{i,j}^{k}\frac{u_{i+1,j}^{k+1}-u_{i,j}^{k-1}}{2\Delta y}+z_{i,j}^{k}u_{i,j}^{k }\frac{v_{i,j+1}^{k+1}-v_{i,j-1}^{k}}{2\Delta y}=gz_{i,j}^{k}\Big{(}S_{0x}-S_ {fx}\Big{)}\end{split}  \tag{5}\]\[ \begin{split} v_{i,j}^{k}\frac{z_{i,j}^{k+1}-z_{i,j}^{k-1}}{2\Delta t }+z_{i,j}^{k+1}\frac{v_{i,j}^{k+1}-v_{i,j}^{k-1}}{2\Delta t}+\left(v_{i,j}^{k} \right)^{2}\frac{z_{i,j+1}^{k}-z_{i,j-1}^{k}}{2\Delta y}+2z_{i,j}^{k}v_{i,j}^{k} \frac{v_{i,j+1}^{k}-v_{i,j-1}^{k}}{2\Delta y}+gz_{i,j}^{k}\frac{z_{i,j+1}^{k}- z_{i,j-1}^{k}}{2\Delta y}+\\ u_{i,j}^{k}v_{i,j}^{k}\frac{z_{i,j}^{k}-z_{i,j-1}^{k}}{2\Delta x}+z _{i,j}^{k}v_{i,j}^{k}\frac{z_{i,j+1}^{k}-z_{i,j-1}^{k}}{2\Delta x}+z_{i,j}^{k} u_{i,j}^{k}\frac{v_{i,j+1}^{k}-v_{i,j-1}^{k}}{2\Delta x}=gz_{i,j}^{k}\left(S_{0y}-S_{ fy}\right)\end{split}  \tag{6}\]  
Then, after simply shifting the term and removing the denominator, the three variables \(z\), \(u\), and \(v\) can be solved iteratively. As shown in Figure 1, each time an iteration value is calculated, the results of the previous two iterations are used, which is also a feature of the leapfrog format. Because of this, the stability of the calculation process is strengthened. It is also because of the solution value range shown in Figure 1 that the boundary of the grid will be exceeded when calculating the edge value. At this time, we need to introduce boundary conditions and add a circle of ghost cells around the standard grid.