## 2 Fluid Models
### SPH Model  
The SPH framework followed in this work was based on a parallel processing, just as in [15]. Its pseudocode for a single timestep cycle is illustrated in Figure 2. OpenGL library and its Utility Toolkit (GLUT) [16; 17] were used to handle graphics, and OpenCL library [18] to handle GPU scripting that was processed by a graphics board (NVIDIA GFORCE GTX660 ti with 1344 GPU cores) in this work. The left part in Figure 2 is the sequential main loop processed by the CPU, and in the right part, the parallel functions (kernels) processed by the GPU. Particle data in the CPU is stored in three \(N\)-sized vector arrays: position, velocity and force. \(N\) is the number of particles running in the fluid simulation. The CPU functions form a number \(N\) of GPU kernels, passing them into the vector arrays to compute data in parallel. Each CPU function block has to wait for output data from the former block so parallel processing is implemented only for each CPU function block. After vector array initialization, a nearest neighbor search algorithm (NNS) [15] is needed to provide the surrounding particles that will contribute to the force acting on particle \(i\). Initially, in NNS, the 2D domain is divided into cells using a uniform grid so a cell index can store each particle where it belongs. Then, a spatial hash value is created using the previous cell index, and particles are sorted according to this hash value. A bi-tonic sorting algorithm is used to take advantage of the GPU processing [19]. Finally, the vector arrays of position,  
Figure 2: One time-step CPU-GPU implementation of smoothed particle hydrodynamics (SPH) fluid motion. To the left are the sequential tasks processed by the CPU. To the right are the parallel tasks (kernels) processed by the GPU. velocity and force are permuted according to their sorted hash values. The NSS at the end will provide a space map on which particles are located so that finding neighbors can be iterated as quickly and straightforwardly as possible. The less time-consuming the NNS algorithm is, the more helpful to the next parts of the SPH implementation it will be. This is why GPU processing turns out to be fundamental to computing the hashing, the sorting and the permuting. OpenCL threads can be synchronized to optimize global and local GPU memory transfers. The implementations of this work are available from NVIDIA [19] and [20]. The next block in Figure 2 is the calculation of liquid density from the individual SPH particles. Equation (1) defines density as the sum of masses weighted by a \(W\) function. \(W\) is the interpolation kernel that computes the position of particles \(i\) and \(j\) within the distance \(h\), also defined as the smoothing distance by SPH formalisms. A cubic spline function is used for \(W\); this provides compact support in its derivatives, just as in [7], and efficiently smooths away the acting neighboring forces with a relatively minor computational cost. Density is the first to compute because it is needed by the next parallel computed forces; see Figure 2.  
\[\rho_{i}=\sum\
olimits_{j}m_{i}W\big{(}\tau_{ij},h\big{)}. \tag{1}\]  
After processing density for every fluid particle using Equation (1), a force vector for each particle is computed for pressure, viscosity and external (gravity and collisions) force contributions. Pressure \(p\) and viscosity \(v\) come from the SPH interpolation formulation, just as in [21]; see the Equations (3) and (4) used. Using the same \(W\) function as in Equation (1) but with first and second derivative order, the gravity force is a \(y\)-axis constant vector [0 9.81 0] m/s\({}^{2}\).