## Archive for the ‘fluid simulation’ tag

## Deep Opacity Maps

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) break; // 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) break; 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:

## 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; EmitVertex(); gl_Position = gl_in[1].gl_Position; EmitVertex(); gl_Position = gl_in[2].gl_Position; EmitVertex(); EndPrimitive(); }

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;

Simple!

## Downloads

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!