The Fluid Simulation Algorithm used in Pixelflow is based on Jos Stam’s paper Real-Time Fluid Dynamics for Games.

A very detailed article for a GPU-Version can be found here: http://http.developer.nvidia.com/GPUGems/gpugems_ch38.html

My Implementation is basically a port of the Cg-Shaders to GLSL-Shaders. So this posting just covers a quick summary of the techniques used and what differs to the original article.

Fluid Solver
update step
{
  advection: velocity, density, temperature, particles, etc.
  forces: buoyancy, vorticity, user input, etc.
  projection
    divergence
    pressure (jacobi solver)
    gradient-subtraction
}
Advection

Advection “is the process by which a fluid’s velocity transports itself and other quantities in the fluid”.

At each grid-cell, the current velocity-vector is used for stepping back and looking up the quantity to advect. Velocity itself is advected this way and replaces the old velocity at the current grid-cell.

The GPU-Gems article executes a separate Diffusion-step, which affects the fluids viscosity (velocity diffusion). Same as the advection step, diffusion can be applied for each quantity. However, in Pixelflow, diffusion is disabled to save performance, and instead a simple damping parameter is used during the advection step, called dissipation. It gives similar results, but avoids a lot of additional jacobi-iterations.

Forces

This is where additional quantities are added to the fluid-textures. External velocity, density, temperature, vorticity, buoyancy, etc…

Projection

divergence “is the rate at which “density” exits a given region of space. In the Navier-Stokes equations, it is applied to the velocity of the flow, and it measures the net change in velocity across a surface surrounding a small piece of the fluid.” The divergence of a vector field (fluid velocity) is a scalar field. It is used as input for the Jacobi solver to compute the Pressure. The jacobi solver requires multiple iterations to give decent results. Finally the current Velocity gets updated by subtracting the Pressure-Gradients from the intermediate Velocity.

 

OpenGL Textures (PING Pong)

 

Particles (custom advection)
StreamLines

The Streamlines-example is similar to the GPU particle simulation, but instead of points (or sprites) curves, based on the fluid-velocity, get updated and displayed.

At certain locations (regular grid points, e.g. every 5-10 pixels) a curve has a start point. The subsequent curve vertices are computed by looking up the velocity from the fluid-solver at the current location. In the end this is quite an expensive task in multiple ways.
To keep the VRAM usage low, the vertex locations are computed on the fly, so only one RGBA32F texel is required per curve. This also means, that the drawing and vertex updating happens in an iterative + alternating fashion.

 

 

GLSL-Shader – Advection
/**
 * 
 * PixelFlow | Copyright (C) 2016 Thomas Diewald - http://thomasdiewald.com
 * 
 * A Processing/Java library for high performance GPU-Computing (GLSL).
 * MIT License: https://opensource.org/licenses/MIT
 * 
 */

#version 150 

precision mediump float;
precision mediump int;

out vec4 glFragColor;

uniform sampler2D tex_velocity;
uniform sampler2D tex_source;
uniform sampler2D tex_obstacleC;

uniform vec2  wh_inv;
uniform float timestep;
uniform float rdx;
uniform float dissipation;

void main(){
  vec2 posn = gl_FragCoord.xy * wh_inv;
  
  float oC = texture(tex_obstacleC, posn).x;
  if (oC == 1.0) {
    glFragColor = vec4(0);
  } else {
    vec2 velocity = texture(tex_velocity, posn).xy;
    vec2 posn_back = posn - timestep * rdx * velocity * wh_inv;
    glFragColor = dissipation * texture(tex_source, posn_back);
  }
}
GLSL-Shader – Divergence
/**
 * 
 * PixelFlow | Copyright (C) 2016 Thomas Diewald - http://thomasdiewald.com
 * 
 * A Processing/Java library for high performance GPU-Computing (GLSL).
 * MIT License: https://opensource.org/licenses/MIT
 * 
 */


#version 150

precision mediump float;
precision mediump int;

out float glFragColor;

uniform sampler2D tex_velocity;
uniform sampler2D tex_obstacleC;
uniform sampler2D tex_obstacleN;

uniform vec2  wh_inv;
uniform float halfrdx;

void main(){
  
  vec2 posn = wh_inv * gl_FragCoord.xy;
  
  float oC = texture(tex_obstacleC, posn).x;
  if (oC == 1.0) { 
    glFragColor = 0.0; 
    return;
  }
  
  // velocity
  vec2 vT = textureOffset(tex_velocity, posn, + ivec2(0,1)).xy;
  vec2 vB = textureOffset(tex_velocity, posn, - ivec2(0,1)).xy;
  vec2 vR = textureOffset(tex_velocity, posn, + ivec2(1,0)).xy;
  vec2 vL = textureOffset(tex_velocity, posn, - ivec2(1,0)).xy;
  vec2 vC = texture      (tex_velocity, posn).xy;
  
  // no-slip (zero) velocity boundary conditions
  // use negative center velocity if neighbor is an obstacle
  vec4 oN = texture(tex_obstacleN, posn);
  vT = mix(vT, -vC, oN.x);  // if(oT > 0.0) vT = -vC;
  vB = mix(vB, -vC, oN.y);  // if(oB > 0.0) vB = -vC;
  vR = mix(vR, -vC, oN.z);  // if(oR > 0.0) vR = -vC;
  vL = mix(vL, -vC, oN.w);  // if(oL > 0.0) vL = -vC;
  
  glFragColor = halfrdx  * ((vR.x - vL.x) + (vT.y - vB.y));
}
GLSL-Shader – Jacobi (Pressure)
/**
 * 
 * PixelFlow | Copyright (C) 2016 Thomas Diewald - http://thomasdiewald.com
 * 
 * A Processing/Java library for high performance GPU-Computing (GLSL).
 * MIT License: https://opensource.org/licenses/MIT
 * 
 */


#version 150

precision mediump float;
precision mediump int;

out vec4 glFragColor;

uniform sampler2D tex_x;
uniform sampler2D tex_b;
uniform sampler2D tex_obstacleC;
uniform sampler2D tex_obstacleN;

uniform vec2  wh_inv;
uniform float alpha;
uniform float rBeta;

void main(){

  vec2 posn = wh_inv * gl_FragCoord.xy;
  
  float oC = texture(tex_obstacleC, posn).x;
  if (oC == 1.0) { 
    glFragColor = vec4(0); 
    return;
  }
  
  // tex b
  vec4 bC = texture(tex_b, posn);
  
  // tex x
  vec4 xT = textureOffset(tex_x, posn, + ivec2(0,1));
  vec4 xB = textureOffset(tex_x, posn, - ivec2(0,1));
  vec4 xR = textureOffset(tex_x, posn, + ivec2(1,0));
  vec4 xL = textureOffset(tex_x, posn, - ivec2(1,0));
  vec4 xC = texture      (tex_x, posn);

  // pure Neumann pressure boundary
  // use center x (pressure) if neighbor is an obstacle
  vec4 oN = texture(tex_obstacleN, posn);
  xT = mix(xT, xC, oN.x);  // if (oT > 0.0) xT = xC;
  xB = mix(xB, xC, oN.y);  // if (oB > 0.0) xB = xC;
  xR = mix(xR, xC, oN.z);  // if (oR > 0.0) xR = xC;
  xL = mix(xL, xC, oN.w);  // if (oL > 0.0) xL = xC;
  
  glFragColor = (xL + xR + xB + xT + alpha * bC) * rBeta;
}
GLSL-Shader – Gradient Subtraction
/**
 * 
 * PixelFlow | Copyright (C) 2016 Thomas Diewald - http://thomasdiewald.com
 * 
 * A Processing/Java library for high performance GPU-Computing (GLSL).
 * MIT License: https://opensource.org/licenses/MIT
 * 
 */


#version 150

precision mediump float;
precision mediump int;

out vec2 glFragColor;

uniform sampler2D tex_velocity;
uniform sampler2D tex_pressure;
uniform sampler2D tex_obstacleC;
uniform sampler2D tex_obstacleN;

uniform vec2  wh_inv;
uniform float halfrdx;

void main(){
    
  vec2 posn = wh_inv * gl_FragCoord.xy;
  
  float oC = texture(tex_obstacleC, posn).x;
  if (oC == 1.0) { 
    glFragColor = vec2(0); 
    return;
  }
  
  // pressure
  float pT = textureOffset(tex_pressure, posn, + ivec2(0,1)).x;
  float pB = textureOffset(tex_pressure, posn, - ivec2(0,1)).x;
  float pR = textureOffset(tex_pressure, posn, + ivec2(1,0)).x;
  float pL = textureOffset(tex_pressure, posn, - ivec2(1,0)).x;
  float pC = texture      (tex_pressure, posn).x;
  
  // pure Neumann pressure boundary
  // use center pressure if neighbor is an obstacle
  vec4 oN = texture(tex_obstacleN, posn);
  pT = mix(pT, pC, oN.x);  // if (oT > 0.0) xT = xC;
  pB = mix(pB, pC, oN.y);  // if (oB > 0.0) xB = xC;
  pR = mix(pR, pC, oN.z);  // if (oR > 0.0) xR = xC;
  pL = mix(pL, pC, oN.w);  // if (oL > 0.0) xL = xC;
  
  // gradient subtract
  vec2 grad = halfrdx * vec2(pR - pL, pT - pB);
  vec2 vOld = texture(tex_velocity, posn).xy;

  glFragColor = vOld - grad;  
}