Simulating Eulerian Fluids Using a Random Basis
The Simulation Method
(These are fundamentally old ideas, but I haven’t seen a simple explanation of this before – Just writing this down for posterity.)
Recall that in fluid simulation, the object of interest is the timeevolution of the velocity vector field \( \mathbf u(\mathbf x, t) \). This object tells us at any point in space and time how fast and in what direction matter is being transported.
A Special Basis
The basic idea of this method is to represent fluid motion over time as a set of a linear superposition of vector fields whose weights are time varying, but the actual vector fields are not. This will give us a decoupling that’s useful later. That is, assume the solution at any point in space and time can be expressed as
We will require two things of our vector field basis functions:
 They must all be divergencefree.
 They must form an orthonormal set. (Really just orthogonal, but normality is convenient.)
These weights are just coordinates in the basis provided by the \(\mathbf \phi\)s.
Integrating the Basis Into the Dynamics
Setting aside for a moment how to actually compute a good basis with the abovementioned properties, the basic gist of the method works like this…
Take the equations for inviscid, incompressible fluid motion
We’ll now substitute using the basis described above, drop the implied spatial and time function arguments, and use Einstein summation to make things easier to follow:
In the first term, because our vector field basis is a constant w.r.t. time, the time derivative of the velocity field turns into a time derivative on the weights. That is, \( \partial_t (\mathbf \phi_i w_i) = \mathbf \phi_i \dot w_i \).
The advection term looks a little more complicated, but it really isn’t. We can commute the weights to the left and we’re left with
In order to get an explicit equation for the dynamics of the reducedbasis weights, we’ll form an inner product with all of the other basis elements:
The first term simplifies due to our orthogonality condition:
The second might have some cool structure (still looking into that…), but it generally looks like the rank3 tensor, \( \mathbf \Pi_{ijp}\):
The fullyassembled, reducedbasis equations are now
Properties of This Construction
Altogether, we don’t have a PDE anymore. The dynamics are now described by a new ODE (still firstorder, quadratic).

Potentially much fewer number of unknowns. Our new dynamics have nothing to do with the number of unknowns in the original vector space. We can choose the number of basis vector fields to tradeoff between fast simulation and our ability to express more varied solutions.

All of our solutions are divergencefree, which means all our weightedsum will be divergencefree, which means… We no longer need to spend time enforcing the divergencefree condition during simulation. We can compute the \(\mathbf\Pi\) tensor before simulation starts and then evaluate the timederivative of the weights in \(\mathcal O(r^2)\), with \(r\) being the number of basis vectors.

Unfortunately… Like I just said, evaluation of \(\dot w\) is still quadratic in the number of basis vectors. This entire technique only works so long as the number of basis is ~\(\sqrt{N}\), with \(N\) being the number of unknowns in the fulldimensional problem.

The same logic above can be used to project forces and viscous effects.
Current Problems

For some reason, a simple Euler Method integration of this seems to be stable, but looks nothing like I’d expect… Hmm..

Is there a special integrator for this type of ODE?