## Time evolution

Simulations evolve forward in time in discrete steps. Time
evolution schemes include families of explicit and implicit techniques that
have various accuracy and stability properties. The most straightforward scheme
entails explicitly and directly integrating values forward in time, given
information about previous times. The same schemes apply to ordinary
differential equations used to simulate particle and rigid body motion.

Explicit methods solve for the future using information
about the past. Varieties include Euler (Figure 8a), Runge-Kutta, and Midpoint
(Figure 8b). The simplest method -- the forward Euler method -- is defined as

where is the quantity being evolved, is the time step, is the first time derivative, and is time. The forward Euler method is analogous
to using a forward difference scheme (as in Figure 7a). The example PDE above,
discretized in both time (using forward difference) and in space (using
centered difference), becomes

**Figure 8: Explicit
timestepping. (a) Explicit Euler. (b) Midpoint "leapfrog".
**

Runge-Kutta methods combine multiple smaller steps to
create a bigger step whose error is smaller than the combined errors of the
smaller steps. The Midpoint method is analogous to using the centered-difference
scheme (as in Figure 7c). It has higher-order accuracy than Euler with similar
computational cost. Position and velocity "leapfrog" each other, so you would
effectively keep track of two staggered simulations.

Applying this to each node *i* makes a system of linear algebraic equations. You could express
that system of equations in matrix-vector form. The resulting matrix would be
sparse (that is, have a lot of elements that are zero). Specialized linear algebra
algorithms exist that exploit sparseness: You could use them to solve that
system of equations.

Explicit methods are simple, but they can be tricky or slow
to make stable. The Courant-Friedrichs-Lewy (CFL) condition establishes a
relationship between spatial resolution and time step in some PDEs. Informally,
this means that the time step must be smaller than the duration it takes for a
bit of fluid to move from one grid cell to another (as Figure 9 shows). That
implies that the faster the fluid motion, the smaller the time step, which
implies more time steps, which implies more computational effort and a slower
simulation.

If a numerical solution is unstable, errors grow rapidly
without bound. One way to mitigate this is to damp out oscillations (for
example, using viscosity). This prevents the solution from "exploding" but also
restricts fluid motion to have high viscosity -- for example, syrup instead of
water.

You can mitigate the results of overdamping by injecting
fine-grain detail lost as a result of excessive viscosity that was introduced
to stabilize the simulation. Such techniques include large-eddy simulation
(LES), Reynolds-averaged Navier-Stokes (RANS), detached eddy, and vorticity
confinement.

Semi-Lagrangian techniques (described above) avoid the
instability by using a particle-based scheme to advect fluid (as Figure 9c
depicts). The other terms of the fluid dynamics equations are computed using an
Eulerian view. Backtracking dodges the instability even when fluid parcels move
farther than a grid cell per time step.

**Figure 9: CFL
condition. (a) CFL violated, (b) CFL satisfied with smaller timesteps, and (c)
instability avoided using backtracking.
**

Implicit methods derive from writing the equations of
motion so that the future solution depends on both future and past solutions. For
example, the backward Euler method is defined as

This might seem absurd: You would apparently need to know
the future to solve the problem. But in practice, it is possible to find a solution
implicitly by using a numerical linear algebraic solver, such as Jacobi,
Gauss-Seidel, successive over-relaxation (SOR), or conjugate gradient (CG). When
using a finite difference discretization and its corresponding linear system of
equations (described above), using an implicit scheme would only require
changing the values in the linear algebraic system's matrix; the solver code
would remain identical.

Implicit methods are often stable when explicit methods are
not, even if the simulation apparently violates the CFL condition. But this
stability comes with a cost: Effectively, the simulation becomes damped,
behaving somewhat as though it cranked up viscosity. This phenomenon is called numerical viscosity.
The final results might end up looking like those from an explicit solver with
very high viscosity -- and the implicit solver can take more computational effort.

**Enforcing Solenoidality
**

Usually for visual effects, fluid simulations use the
"incompressible" approximation, meaning mass cannot converge or diverge. In
this case, there is no separate equation for pressure: It is coupled directly
to velocity via incompressibility. The advection step can violate
incompressibility, so it needs to be restored. One solution is to use a "scalar
Poisson solver" to project the divergent field to obtain its solenoidal
component. This is what Stam uses in "stable fluids." In contrast, Mick West
constructs the advection step to conform to the incompressibility constraint
and so avoids this separate step.

An analogous problem arises in vortex-based simulations: The
vorticity field can accumulate divergence, even though mathematically vorticity
(or any curl) cannot have divergence. This divergence does not influence the
velocity field, but it can alter vortex stretching. The simulation must
formulate vortex stretching such that divergent vorticity does not cause undue
problems. This requires some straightforward algebraic manipulation, but for
the sake of brevity, this article omits a detailed description.