The Little Grasshopper

Archive for the ‘fluid simulation’ tag

Deep Opacity Maps

without comments

In the video, the light source is slowly spinning around the smoke stack, showing that a deep opacity map is generated at every frame. For now it’s being computed rather inefficiently, by performing a full raycast over each voxel in the light map. Here’s my fragment shader for generating the opacity map:

in float gLayer;
out float FragColor;

uniform sampler3D Density;
uniform vec3 LightPosition = vec3(1.0, 1.0, 2.0);
uniform float LightIntensity = 10.0;
uniform float Absorption = 10.0;
uniform float LightStep;
uniform int LightSamples;
uniform vec3 InverseSize;

void main()
    vec3 pos = InverseSize * vec3(gl_FragCoord.xy, gLayer);
    vec3 lightDir = normalize(LightPosition-pos) * LightStep;
    float Tl = 1.0;
    vec3 lpos = pos + lightDir;
    for (int s = 0; s < LightSamples; ++s) {
        float ld = texture(Density, pos).x;
        Tl *= 1.0 - Absorption * LightStep * ld;
        if (Tl <= 0.01)

        // Would be faster if this conditional is replaced with a tighter loop
        if (lpos.x < 0 || lpos.y < 0 || lpos.z < 0 ||
            lpos.x > 1 || lpos.y > 1 || lpos.z > 1)

        lpos += lightDir;

    float Li = LightIntensity*Tl;
    FragColor = Li;

It would be much more efficient to align the light volume with the light direction, allowing you to remove the loop and accumulate results from the previous layer. Normally you can’t sample from the same texture that you’re rendering to, but in this case it would be a different layer of the FBO, so it would be legal.

Another potential optimization is using MRT with 8 render targets, and packing 4 depths into each layer; this would generate 32 layers per instance rather than just one!

Here’s another video just for fun, taken with a higher resolution grid. It doesn’t run at interactive rates (I forced the video to 30 fps) but the crepuscular shadow rays sure look cool:

Written by Philip Rideout

February 16th, 2011 at 2:55 pm

Posted in OpenGL

Tagged with

3D Eulerian Grid

without comments

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