3D Eulerian Gas Simulation

Howdy! I'm Joshua Zimmerman and welcome to my website! This page will help break down what I did for my final project in CSCE 450.

Presentation Video

Goal of This Project

I set out to:

• Make a 3D Eulerian simulation in C++ for gasses to mimic examples I had seen on youtube like this
• Be able to export simulation data to be rendered by an animation software like Blender
• Define Simulation constants in one place to change parameters on the fly depending on a user's preferences
• Have Blender volume shader compatibility to be able to change the gas color, apply lighting, and fully render an animation or stillframe using Cycles or EEVEE

Libraries and External Software Used

Libraries I used in my program and why:

• Open VDB: Used for openvdb::FloatGrid to export each cell's density as a value in a widely used volumetric data format
• Eigen: Used for Vectors, general Vector operations, Sparse Matrices, and ConjugateGradient solvers
• Header Only Perlin Noise: Used to randomize injected density to help the simulation not be so static
• Blender: Used to visualize and render exported openvdb grid objects as frame sequences
• I also used Cmake for compatibility, Visual Studio, and Cycles within Blender for renders

A Program Overview

A breakdown of core program functionality:

• This is a timestep based program where each frame is a small step forward in the simulation.
• The class structure is as follows:
  • Cell - Holds Temperature (scalar), Velocity (vector), Density (scalar), and Pressure (scalar)
  • Grid - Holds GRID_HEIGHT * GRID_WIDTH * GRID_DEPTH number of cells. This has methods for get/set/trilinearly sampling at a certain x y z. Also has a grid to vdb conversion function.
  • Simulator - Holds a grid and houses private methods for manipulating the grid
• A simulation step consists of the following processes in order:
• Injection Step - Uses Perlin noise to inject semi-random density at the bottom of the Grid. A constant temperature is also given to all non-zero density cells.
• External Force Step (Optional) - I simulate wind using Perlin noise as a variable strength force in a predetermined and slightly wavering direction.
• Buoyancy - Apply an upward velocity to cells based on the buoyancy density scale and buoyancy temperature scale. The temperature component is decreased depending on how close it is to the ambient temperature of the environment.
• Vorticity Confinement - Adds small scale swirls back into the simulation that Semi-Lagrangian advection typically smooths out. The algorithm i use computes the curl of the gas then adds a force perpendicular to both the curl and the gradient of the curl.
• Diffuse Velocity - This step spreads the velocity out over time. I use the linear system
  (I - α∇²)uⁿ⁺¹ = uⁿ
  to smooth the spread of velocity from grids to its neighbors. Velocity diffusion mimics viscosity.
• Project Pressure - Enforce incompressibility after changing velocity throughout the simulation. Without incompressibility, the simulation has divergence meaning gas will disappear over time.
• Advect Velocity - This step is what moves velocity with the gas. Think of a gas carrying its own motion vector which is what makes it move through the grid!
• Project Pressure - Enforce incompressibility after changing velocity again.
• Diffuse Density (Optional) - Spread out the density of the gas into neighboring cells. Causes dissipation making gas sometimes appear to vanish.
• Diffuse Temperature - Spread out the temperature of the gas into neighboring cells.
• Advect Density - Now I move the density field along the velocity, moving the gas visually in the simulation.
• Advect Temperature - Finally I move the temperature field along the velocity, shifting the temperature of the gas with its movement.
• Enforce Boundary Conditions (Optional) - Either contain the gas or let it escape.
• The grid after each simulation step is converted to a vdb file and written to a specified directory. This is done NUM_FRAMES number of times.

Some Development Notes

Some tips for someone who wants to implement a similar system:

• The main difficulty of this program is discretizing the Laplace operator. Check out this source for an introduction to the idea. Notice that in 2D, we can use a typical 2D sparse matrix to represent a matrix of laplace operators but in 3D we need a flattened coordinate system and then use triplet construction on top of that.
• A 7 point stencil for 3D using Eigen this looks a little something like:

         Eigen::SparseMatrix<double> Amatrix(gridSize, gridSize);
          std::vector<Eigen::Triplet<double>> coefficients;
          for x y z in grid {
              uint32_t i = flattenIndex(x, y, z);
              coefficients.emplace_back(i, i, 6);
              coefficients.emplace_back(i, flattenIndex(x + 1, y, z), -1);
              coefficients.emplace_back(i, flattenIndex(x - 1, y, z), -1);
              coefficients.emplace_back(i, flattenIndex(x, y + 1, z), -1);
              coefficients.emplace_back(i, flattenIndex(x, y - 1, z), -1);
              coefficients.emplace_back(i, flattenIndex(x, y, z + 1), -1);
              coefficients.emplace_back(i, flattenIndex(x, y, z - 1), -1);
          }
          Amatrix.setFromTriplets(coefficients.begin(), coefficients.end());

• Pressure projection enforces incompressibility by solving a Poisson equation derived from the velocity divergence:
  ∇²p = ∇ · u*
  where:
  • ∇ · u* is the divergence of the intermediate velocity field (after advection and diffusion)
  • To estimate velocity divergence using a scalar at a grid cell using finite differences:

    
              Eigen::Vector3d vxp = grid.getVelocity(x + 1, y, z);
              Eigen::Vector3d vxm = grid.getVelocity(x - 1, y, z);
              Eigen::Vector3d vyp = grid.getVelocity(x, y + 1, z);
              Eigen::Vector3d vym = grid.getVelocity(x, y - 1, z);
              Eigen::Vector3d vzp = grid.getVelocity(x, y, z + 1);
              Eigen::Vector3d vzm = grid.getVelocity(x, y, z - 1);
              double divergence = 0.5 * (
                  vxp.x() - vxm.x() +  // ∂u/∂x
                  vyp.y() - vym.y() +  // ∂v/∂y
                  vzp.z() - vzm.z()    // ∂w/∂z
              );
              b[flattenIndex(x, y, z)] = -divergence; // b in Ax = b system
                

  • ∇²p is the Laplacian of the pressure field, which we solve for
  • Subtracting ∇p from u* gives the updated, divergence-free velocity: uⁿ⁺¹ = u* - ∇p
• The diffusion step solves a linear system based on the discretized diffusion equation:
  (I - α∇²)uⁿ⁺¹ = uⁿ
  where:
  • I is the identity matrix
  • α = Δt · ν (Δt is the timestep, ν is the viscosity)
  • uⁿ⁺¹ is the velocity after diffusion
• Once your Ax = b system is set up for diffusion or pressure projection, you can use Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> as a solver for your linear system to get a Eigen::VectorXd that you can access with solutionVector[flattenIndex(x, y, z)]
• For Vorticity confinement, imagine a helper function:

     Eigen::Vector3d Simulator::computeVorticity(int x, int y, int z) {
      Eigen::Vector3d vxp = grid.getVelocity(x + 1, y, z);
      Eigen::Vector3d vxm = grid.getVelocity(x - 1, y, z);
      Eigen::Vector3d vyp = grid.getVelocity(x, y + 1, z);
      Eigen::Vector3d vym = grid.getVelocity(x, y - 1, z);
      Eigen::Vector3d vzp = grid.getVelocity(x, y, z + 1);
      Eigen::Vector3d vzm = grid.getVelocity(x, y, z - 1);
  
      double dw_dy = (vzp.y() - vzm.y()) * 0.5;
      double dv_dz = (vyp.z() - vym.z()) * 0.5;
      double du_dz = (vxp.z() - vxm.z()) * 0.5;
      double dw_dx = (vzp.x() - vzm.x()) * 0.5;
      double dv_dx = (vyp.x() - vym.x()) * 0.5;
      double du_dy = (vxp.y() - vxm.y()) * 0.5;
  
      return Eigen::Vector3d(dv_dz - dw_dy, dw_dx - du_dz, du_dy - dv_dx);
    }

  And vorticity force equation: f = ε ( ∇|ω| × ω ) where:
  • ε is a vorticity confinement scale (strength of vorticity)
  • ω is the vorticity vector
  • ∇|ω| is the normalized vorticity magnitude gradient
  Then we can compute the cofinement force as follows:

     Eigen::Vector3d vort = computeVorticity(x, y, z);
      double nx = computeVorticity(x + 1, y, z).norm() - computeVorticity(x - 1, y, z).norm();
      double ny = computeVorticity(x, y + 1, z).norm() - computeVorticity(x, y - 1, z).norm();
      double nz = computeVorticity(x, y, z + 1).norm() - computeVorticity(x, y, z - 1).norm();
      Eigen::Vector3d gradMag(nx, ny, nz);
      gradMag = 0.5;
      Eigen::Vector3d confinement = gradMag.cross(vort) VORTICITY_CONFINEMENT_SCALE * dt;

Future Work

If I had another semester with this project I would aim to:

• Implement specific wind or temperature based scenes. Doing this could create something cool like modeling a tornado or hurricane
• Implement simple gas-object collisions (think gas moving around a sphere) and fix boundary condition behavior
• Export velocity in each cell as its own volumetric grid and make a render where blender colors gas based on its speed
• Automate rendering and uploading frames to blender with scripts
• Refine buoyancy force to allow gas to properly sink after cooling past a certain point