Brandon Pelfrey Mathematics and Computational Curiosities

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 time-evolution 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:

  1. They must all be divergence-free.
  1. 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 above-mentioned 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 reduced-basis 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 rank-3 tensor, \( \mathbf \Pi_{ijp}\):

The fully-assembled, reduced-basis 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 first-order, 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 trade-off between fast simulation and our ability to express more varied solutions.

  • All of our solutions are divergence-free, which means all our weighted-sum will be divergence-free, which means… We no longer need to spend time enforcing the divergence-free condition during simulation. We can compute the \(\mathbf\Pi\) tensor before simulation starts and then evaluate the time-derivative 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 full-dimensional 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?