Simulating Fluids in CG

1
Simulating Fluids in CG

• Advection Algorithms

2
Last Time

• The equations of fluids
• Pressure
• Hydrostatic equilibrium
• Differential operators Gradient, divergence
• Incompressible Euler Equations
• Gauss Theorems

3
This Lecture

• HOW do we solve the Incompressible Euler
  equations
• Numerical methods
• Time integration
• Splitting
• Finite differences
• The MAC grid
• Advection algorithms
• Adding forces

4
Recap
5
Time Integration I

• We are dealing with equations of the form
• Given an initial condition we want to
  solve this by computing at discrete times
  with
• The time step is the length of time between these
  discrete times

6
Time Integration II

• The simplest thing to do is to assume that is
  constant during the entire time step
• This is called forward Euler integration
• Only first-order accurate

7
Time Step Size

• Stability
• Many time integration schemes are only
  conditionally stable
• If ?t is chosen too large, the computed solution
  blows up exponentially even though the exact
  solution stays bounded
• Convergence
• A large ?t can give us the wrong solution (this
  is less critical in CG, as long as it looks good)

8
Splitting I

• Consider the simple ordinary differential
  equation
• The solution is but lets try splitting it
  using forward Euler

9
Splitting II

• Lets generalize
• And split
• F and G are specialized integrators

10
Splitting the Euler Equations
(advection / transport)
(body forces)
(pressure / incompressibility)
11
Body forces

• Forces accelerate the fluid. In the model we
  consider it is the only way to initiate and
  control the movement of the fluid.

We are solving the equation
Like this
12
Basic fluid algorithm

• Start with initial, divergence-free velocity
  field v0
• For time n 0, 1, 2,
• Determine time step ?t from time tn to time tn1
• vA advect(vn, ?t, vn) (advection)
• vB vA ?t f (body forces)
• vn1 project(?t, vB) (incompressibility)

13
Advection

• We need to develop an algorithm advect to solve
  the equation

14
Discretizing the Differential Operators
Question We have just considered the
differential operators in R3, but how do we apply
these operators in the discrete simulation space?
Answer Finite difference approximations based on
Taylor series.
15
Taylor Series

• 1D
• 2D
• Generalizes to ND.

16
Finite Differences I
The Taylor series
Can be re-arranged into
If we assume that ?x is the grid spacing (always
positive) we can form Finite differences like
this
17
Finite Differences II

• Assume we have a function of two variables
• We can find finite difference approximations to
  mixed partial derivatives

18
Finite Differences III

• The collection of discrete points that we use to
  compute the finite differences is called the
  stencil.

19
The MAC Grid I

• In a standard grid we sample pressure and
  velocities in the cell centers

20
The MAC Grid II

• The staggered/MAC moves the components of the
  velocities to the cell faces

21
The MAC Grid III

• Convention (now in 3D)

22
Back to Advection

• Lets try applying finite differences

23
Problems!!

• Unconditionally unstable!
• It will always blow up
• High-frequency components do not get advected as
  fast as low-frequency components
• Artifacts such as oscillations and wiggles

24
Semi-Lagrangian Advection I

• Recall the Lagrangian framework
• Advection is trivial in this view
• Simply moving the particles solves the equation
• And to get the new value of a quantity q at
  position x, one could just find the particle that
  ends up at x and look at its value of q.

25
Semi-Lagrangian Advection II
To find the value of a given voxel at time t?t,
we trace the velocity field backwards in time to
time t, and take the value from there.
The Velocity Field
26
Semi-Lagrangian Advection III
So we will trace the velocity backwards.
We consider two ways of doing this
p
1 A first order accurate backwards Euler
time-step
2 A second order accurate backwards Runge Kutta
time-step
To find the value of a given voxel at time t?t,
we trace the velocity field backwards in time to
time t, and take the value from there.
27
Trilinear Interpolation
Note that the back-traced position may not lie
exactly at the center of a cell or on a face, so
we need to Interpolate surrounding values. If the
back-traced position ends up outside the grid, we
just assume the value is zero.
We look up the values at the eight closest cell
centers/faces and interpolate linearly in each
dimension hence the name tri-linear
interpolation in 3D.
The interpolation simplifies to the
following formula where (x,y,z) is the offset
(measured in cells, ie x,y and z is between zero
and one) from q000
q
q
q
q
q
q
q
q
28
Numerical Dissipation

• With each advection step, tri-linear
  interpolation performs an averaging (i.e.
  blurring) operation
• Blurring every time step gets bad!
• Analyzing the effect, shows us that we are in
  fact introducing an additional viscosity-like
  term to the equation were solving!

29
Higher-order Interpolation

• Cubic Hermite spline (in 1D)
• Derivatives are approximated by central finite
  differences

30
Stability concern
31
Advection Loop

• Procedure advect ()
• for (each voxel in grid at time t?t)
• Back-trace voxel-position to time t using current
  velocity.
• Lookup values in grid at time t at neighboring
  voxels of back-traced voxel-position.
• Do some interpolation between these neighboring
  values.
• Assign the result of the interpolation to the
  current voxel at time t?t.
• end
• end

32
Smoke densities

• Advection of velocity field components v and
  smoke densities d is done in essentially the same
  way

33
Loop For Moving Smoke Densities

• Keep two density grids in memory. Needed for
  storing intermediate results. Remember to swap
  the density grids.
• For (number of simulation iterations)
• addDensitySource()
• advectDensity()
• End

34
Additional Resources

• Book Mathematical Methods for Physicists
• Web www.mathworld.com and http//en.wikipedia.org
  

35
Next Time

• Solving for incompressibility
• Exercises in advection