Fluid Simulations on GPU with Complex Boundary Conditions Youquan Liu1,3, Xuehui Liu1, Enhua Wu1, 2 Laboratory of Computer Science, Institute of Software, Chinese Academy of Sciences, Beijing, China 2 Department of Computer and Information Science, Faculty of Science and Technology, University of Macau, Macao, China 3 Graduate School of the Chinese Academy of Sciences, Beijing, China [email protected], [email protected], [email protected] http://lcs.ios.ac.cn/~lyq/demos/gpgpu/gpgpu.htm 1

1 Introduction Fluid simulation is always a hot topic in computer graphics community. There is much work focused on this [1,2,3], trying to simulate water flowing, cloud animation, fire, smoke etc. We have implemented a system to simulate fluid flowing in 2D and 3D domain, which can process complex boundary conditions easily and efficiently on GPU.

2 Fluid Simulations on GPU To depict the flowing effects more exactly the Navier-Stokes Equations (NSEs) which mainly include two parts, the continuous equation (1) to ensure mass conservation, and the momentum equation (2) to ensure the momentum conservation are used. ∇ ⋅u = 0 (1) (2) ∂ u ∂ t = − (u ⋅ ∇ )u + v∇ 2u -∇ p + f Here u is the velocity vector, v is the coefficient for control of the diffusion, f is the external force and p is the pressure. Because of its good stability, we adopt semi-Lagrangian method [4] to solve the advection of NSEs without pressure. We compute the Possion equation of pressure, and then correct the velocity field with the pressure. To illustrate the flowing effects, the two scalar variables, density ρ and temperature T are introduced which are passively convected by the velocity field u, similar to the momentum equation. (3) ∂ ρ ∂ t = − (u ⋅ ∇ ) ρ + k ρ∇ 2 ρ + S ρ (4) ∂ T ∂ t = − ( u ⋅ ∇ )T + k T ∇ 2 T + S T To take full advantage of parallelism of GPU, we pack the density and temperature variables into RGBA channels with velocity variable together for 2D problems in our system. For 3D problems, just the density and temperature are packed together and flat 3D texture [2] is used to reduce the rendering pass number.

according to the obstacle information. Node’s Type(i,j) = Obstacle(i,j,k)*64 + Obstacle(i+1,j,k)*1 + Obstacle(i-1,j,k)*2+ Obstacle(i,j+1,k)*8 + Obstacle(i,j-1,k)*4+ Obstacle(i,j,k+1)*16 + Obstacle(i,j,k-1)*32

In this way, the boundary condition can be rewritten as φboundary = dφoffset + e , where d and e are determined from the

equation (1), which are packed to one texture. This texture is used to modify the velocity and pressure on boundaries.

4 Results After computing the whole field, we render the flowing effects with density distribution. So the whole simulation totally happens on GPU, and the CPU is just used to send commands and several quads, or to control the display of result. The test results prove the efficiency of our method. Figure 1 gives an example to illustrate the Von Karman vortex street effect in a classic example in computational fluid dynamics (CFD) field. Compared with the same algorithm on CPU, this method on GPU will speed up to 14 times on our machine. Figure 2 provides another example flowing in a maze, which can be used to simulate the smoke flowing in buildings. More demos can be visited at http://lcs.ios.ac.cn/~lyq/demos/gpgpu/gpgpu.htm. And the 2D application is also downloadable on the site.

Figure 1: Flowing around an obstacle in 2D domain

3 Boundary Processing It’s important to process the boundary conditions for those real fluid problems. The boundary conditions include Dirichlet (prescribed value), Neumann (prescribed derivative), and mixed (coupled value and derivative), as expressed by the equation (1): (1) aφ + b ∂φ ∂n = c Where φ stands for velocity, density, temperature or pressure etc., a, b and c are the coefficients. In order to process the complex boundary conditions on GPU, we take the obstacle information into an image, with 1 standing for fluid, 0 for obstacle. By doing so, we can take advantage of higher parallelism of GPU. It’s easy to implement in 2D domain, but in 3D domain, we clip the bounding box with stencil buffer to form the correct obstacle information that is different from [2]. We discrete equation (1) in first order, then the value of each voxel on boundaries can be determined by one nearest voxel. So we can generate one texture that consists of the voxel offsets that determine the value of the current voxel according to the node type with one fragment program. But before computing the offsets, we firstly generate one texture standing for the type of each node by the following coding scheme (2D problems have the similar form) with one fragment program

Figure 2: Flowing in maze in 3D domain

Acknowledgements This work is financially supported by the NSFC (60223005, 60033010), the National Grand Fundamental Research Grant of Science and Technology (973 Project: 2002CB312102), and the Research Grant of University of Macau.

References 1. 2.

3.

Ronald Fedkiw, Jos Stam and Henrik Wann Jensen. Visual Simulation of Smoke. In Proceedings of SIGGRAPH, pages15-22, August 2001. Wei Li, Zhe Fan, Xiaoming Wei, and Arie Kaufman, GPU-Based Flow Simulation with Complex Boundaries. Technical Report 031105, Computer Science Department, SUNY at Stony Brook, Nov 2003. Jos Stam. Stable Fluids. In Proceedings of SIGGRAPH, pages121-128, July 1999.