The Equations Incompressible Navier-Stokes: t = ( r u ) u 1 u - PowerPoint PPT Presentation

F LUID S IMULATION Kristofer Schlachter

The Equations § Incompressible Navier-Stokes: ∂ t = � ( r · u ) u � 1 ∂ u ρ r p + v r 2 u + F § “Incompressibility condition” r · u = 0

Breakdown ∂ u The derivative of velocity with respect to time . ∂ t Calculated at each grid point each time step. � ( r · u ) u The convection term . This is the self advection term. In the code we will use the backward parti- cle trace for this term. ≈ v r 2 u or The viscosity term . We are actually going to ig- nore this term. When you do that you are actually v r · r u using the Euler Equations. External force . Any external forces including F External force . Any external forces including gravity.

Breakdown cle trace for this term. � 1 The pressure term . ρ is the density of the fluid ρ r p a and p is the pressure. p is whatever it takes to make the velocity field divergence free. The sim- ulator will solve for a pressure that makes our fluid incompressible at each time step. a density of water ρ ≈ 1000 kg m 2

3 Incompressible Euler Equations If you drop the viscosity term from the incompressible Navier Stokes equations we get: ∂ t + ( r · u ) u + 1 ∂ u ρ r p = F (4) r · u = 0 (5) Such an ideal fluid with no viscosity is called inviscid . These are the equations we are going to use.

How do we discretize the equations?

A Simple Grid § We could put all our fluid variables at the nodes of a regular grid § But this causes some major problems ∂ u § In 1D: incompressibility means ∂ x = 0 u i + 1 − u i − 1 § Approximate at a grid point: = 0 2 Δ x § Note the velocity at the grid point isn’t involved! [Bridson 07]

A Simple Grid Disaster ∂ u § The only solutions to are u=constant ∂ x = 0 § But our numerical version has other solutions: u x [Bridson 07]

Staggered Grids § Problem is solved if we don’t skip over grid points § To make it unbiased, we stagger the grid: put velocities halfway between grid points § In 1D, we estimate divergence at a grid point as: u i + 12 − u i − 12 ∂ u ∂ x ( x i ) ≈ Δ x § Problem solved! [Bridson 07]

The MAC Grid § From the Marker-and-Cell (MAC) method [Harlow&Welch’65] § A particular staggering of variables in 2D/3D that works well for incompressible fluids: § Grid cell (i,j,k) has pressure p i,j,k at its center § x-part of velocity u i+1/2,jk in middle of x-face between grid cells (i,j,k) and (i+1,j,k) § y-part of velocity v i,j+1/2,k in middle of y-face § z-part of velocity w i,j,k+1/2 in middle of z-face [Bridson 07]

MAC Grid in 2D p i , j + 1 v i , j + 12 p i − 1, j p i , j p i + 1, j u i − 12 , j u i + 12 , j v i , j − 12 p i , j − 1 [Bridson 07]

Math on a MAC Grid The partial derivative ( ∂ ). tral difference: ∂ f ( x, y, z ) = f ( x, y + 1 , z ) � f ( x, y, z ) ∂ y The gradient operator ( r ). 0 1 f ( x + 1 , y, z ) � f ( x, y, z ) , f ( x, y + 1 , z ) � f ( x, y, z ) , r f ( x, y, z ) = B C @ A f ( x, y, z + 1) � f ( x, y, z )

Math on a MAC Grid The divergence of a vector field ( r · ) produces r · u ( x, y, z ) =( u x ( x + 1 , y, z ) � u x ( x, y, z ))+ ( u y ( x, y + 1 , z ) � u y ( x, y, z ))+ ( u z ( x, y, z + 1) � u z ( x, y, z )) The Laplacian operator ( r 2 ) or ( r · r ). gradient operators. r 2 f ( x, y, z ) = f ( x + 1 , y, z ) + f ( x � 1 , y, z )+ f ( x, y + 1 , z ) + f ( x, y � 1 , z )+ f ( x, y, z + 1) + f ( x, y, z � 1) � 6 f ( x, y, z ) *There might be a one over h squared term missing

Simulation – Set u A = advect ( u n , ∆ t, u n ) – Add u B = u A + ∆ tF – Set u n +1 = project ( ∆ t, u B ) *I am ignoring choosing a time step in this presentation

Advection Figure 3: Basic idea behind the advection step. Instead of moving the cell centers forward in time (b) through the velocity field shown in (a), we look for the particles which end up exactly at the cell centers by tracing backwards in time from the cell centers (c). [Stam 2003]

Advection x mid = x G � 1 2 ∆ tu ( x G ) x p = x G � ∆ tu ( x mid ) q n +1 = interpolate ( q n , x p ) G confused by the Cline[David Cline ] that

Projection The project ( ∆ t, u ) routine does the following: • Calculate the negative divergence b (the right-hand side) • Set the entries of A • If using CG - Construct the MIC(0) preconditioner. • Solve Ap = b with a linear solver. If using CG then solve with MICCG(0), i.e., the PCG algorithm with MIC(0) as pre- conditioner. • Compute the new velocities n n +1 according to the pressure- gradient update to u .

Projection Setting up the divergence vector b (the right hand side) r · u ( x, y, z ) =( u x ( x + 1 , y, z ) � u x ( x, y, z ))+ ( u y ( x, y + 1 , z ) � u y ( x, y, z ))+ ( u z ( x, y, z + 1) � u z ( x, y, z ))

Projection Setting up the matrix � Ω 1 β 1 , 2 . . . β 1 ,n 2 3 p 1 � D 1 2 3 2 3 . . p 2 � D 2 6 7 β 2 , 1 � Ω 2 . 6 7 6 7 6 7 5 = . . 6 7 6 7 6 7 . . . � ... . . . 6 7 4 4 5 . β n − 1 ,n − 1 4 5 p n � D n � Ω n β n, 1 β n,n − 1 · · · Where D i corresponds to the divergences through cell i . Ω i is the number of non- solid neighbors of cell i , and β i,j takes values based upon the equation: ( 1 if cell i is a neighbor of cell j β i,j = 0 otherwise This matrix is well known. It is called a 7 Point Laplacian Matrix

Projection Calculate the pressure gradient and subtract it from the velocity field to ensure it is divergence free: 2 ,j,k � ∆ t 1 p i +1 ,j,k � p i,j,k u n +1 2 ,j,k = u n i + 1 i + 1 ρ ∆ x 2 ,k � ∆ t 1 p i,j +1 ,k � p i,j,k v n +1 2 ,k = v n i,j + 1 i,j + 1 ρ ∆ x 2 � ∆ t 1 p i,j,k + � p i,j,k w n +1 2 = w n i,j,k + 1 i,j,k + 1 ∆ x ρ

Demo

Implementation Mapping to Mechanisms: Stencils Common example (“5-point stencil”): = 1 u n +1 h 2 ( − 4 u n i , j + u n i − 1 , j + u n i +1 , j i , j + u n i , j − 1 + u n i , j +1 ) • Sequential • OpenMP? • MPI? • GPU — 2D? • GPU — 3D? Remember this slide from the last lecture? Visualization Parallel Patterns

Implementation Thoughts on optimization for parallel execution The most time consuming part of the sim is the pressure solve so the optimization should start with solving it quickly. Optimizations for serial & CPU: Blocking: Cache coherent. Can send blocks to different cores Optimizations for GPU: Do blocking again using local memory to store each block. For 3D use slicing, 3 Slices at a time putting middle slice (block) in local memory.