The Little Grasshopper

3D Eulerian Grid

I finally extended my OpenGL-based Eulerian fluid simulation to 3D. It’s not as pretty as Miles Macklin’s stuff (here), but at least it runs at interactive rates on my GeForce GTS 450, which isn’t exactly top-of-the-line. My original blog entry, “Simple Fluid Simulation”, goes into much more detail than this post.

I thought it would be fun to invert buoyancy and gravity to make it look more like liquid nitrogen than smoke. The source code can be downloaded at the end of this post.

Extending to three dimensions was fairly straightforward. I’m still using a half-float format for my velocity texture, but I’ve obviously made it 3D and changed its format to GL_RGB16F (previously it was GL_RG16F).

The volumetric image processing is achieved by ping-ponging layered FBOs, using instanced rendering to update all slices of a volume with a single draw call and a tiny 4-vert VBO, like this:

glDrawArraysInstanced(GL_TRIANGLE_STRIP, 0, 4, numLayers);

Here are the vertex and geometry shaders that I use to make this work:

-- Vertex Shader

in vec4 Position;
out int vInstance;

void main()
    gl_Position = Position;
    vInstance = gl_InstanceID;

-- Geometry Shader

layout(triangles) in;
layout(triangle_strip, max_vertices = 3) out;
in int vInstance[3];
out float gLayer;
uniform float InverseSize;
void main()
    gl_Layer = vInstance[0];
    gLayer = float(gl_Layer) + 0.5;
    gl_Position = gl_in[0].gl_Position;
    gl_Position = gl_in[1].gl_Position;
    gl_Position = gl_in[2].gl_Position;

Note that I offset the layer by 0.5. I’m too lazy to pass the entire texture coordinate to the fragment shader, so my fragment shader uses a gl_FragCoord trick to get to the right texture coordinate:

-- Fragment Shader

uniform vec3 InverseSize;
uniform sampler3D VelocityTexture;

void main()
    vec3 tc = InverseSize * vec3(gl_FragCoord.xy, gLayer);
    vec3 velocity = texture(VelocityTexture, tc).xyz;
    // ...

Since the fragment shader gets evaluated at pixel centers, I needed the 0.5 offset on the layer coordinate to prevent the fluid from “drifting” along the Z axis.

The fragment shaders are otherwise pretty damn similar to their 2D counterparts. As an example, I’ll show you how I extended the SubtractGradient shader, which computes the gradient vector for pressure. Previously it did something like this:

float N = texelFetchOffset(Pressure, T, 0, ivec2(0, 1)).r;
float S = texelFetchOffset(Pressure, T, 0, ivec2(0, -1)).r;
float E = texelFetchOffset(Pressure, T, 0, ivec2(1, 0)).r;
float W = texelFetchOffset(Pressure, T, 0, ivec2(-1, 0)).r;
vec2 oldVelocity = texelFetch(Velocity, T, 0).xy;
vec2 gradient = vec2(E - W, N - S) * GradientScale;
vec2 newVelocity = oldVelocity - gradient;

The new code does something like this:

float N = texelFetchOffset(Pressure, T, 0, ivec3(0, 1, 0)).r;
float S = texelFetchOffset(Pressure, T, 0, ivec3(0, -1, 0)).r;
float E = texelFetchOffset(Pressure, T, 0, ivec3(1, 0, 0)).r;
float W = texelFetchOffset(Pressure, T, 0, ivec3(-1, 0, 0)).r;
float U = texelFetchOffset(Pressure, T, 0, ivec3(0, 0, 1)).r;
float D = texelFetchOffset(Pressure, T, 0, ivec3(0, 0, -1)).r;
vec3 oldVelocity = texelFetch(Velocity, T, 0).xyz;
vec3 gradient = vec3(E - W, N - S, U - D) * GradientScale;
vec3 newVelocity = oldVelocity - gradient;



Liquid Nitrogen

I’ve tested the code with Visual Studio 2010. It uses CMake for the build system. If you have trouble building it or running it, let me warn you that you need a decent graphics card and decent OpenGL drivers.

I consider this code to be on the public domain. Enjoy!

Written by Philip Rideout

January 23rd, 2011 at 6:06 pm