Real Time Fluids
This Jul (Swedish for Christmas) I’m not only eating Julbord, Julfika, Julgröt and drinking Julöl (basically put ‘Jul’ before any word and there you go), I’m also doing some Julcoding.
Some time ago I ported a c implementation of a Naiver-Stokes fluid simulation to AS3 and I thought I should share it. It’s all based on this paper ‘Real Time Fluid Dynamics for Games‘ from by Joe Stam. If you want to understand the details about this I suggest you read it, or look around on the web. But the concept is to set up a velocity and density field to which you add external velocities and densities. For each update you solve the velocities by iterating over the grid and then move the densities accordingly.
On top of that I added some balls that you can play around with by moving a ball with the mouse, it’s visually not very impressive but it should give an idea of what you can do with it.
It seems like the flash bottleneck is the memory access (all the _v[index] etc), I think using Alchemy for memory access would speed it up a lot, and you could probably do really cool simulations. If you decide to try it, let me know how it works!
Use the method ‘getDensities’ for rendering (just iterate it and plot it). Here is the AS3 code for the fluid solver:
package com.playchilla.algorithm.fluid
{
public class FluidSolver
{
public function FluidSolver(
width:int,
height:int,
solveIterations:int,
diffusion:Number,
viscosity:Number)
{
_width = width;
_height = height;
_solveIterations = solveIterations;
_diff = diffusion;
_visc = viscosity;
_n2 = _width * _height;
_size = (_width + 2) * (_height + 2);
_d = _createArr();
_d0 = _createArr();
_u = _createArr();
_v = _createArr();
_u0 = _createArr();
_v0 = _createArr();
_rowAdder = _width + 2;
}
public function update(dt:Number):void
{
// Add sources from previous to current (last += new)
_addSource(_u, _u0, dt);
_addSource(_v, _v0, dt);
// Diffuse new with old
_diffuse(1, _u0, _u, _visc, dt);
_diffuse(2, _v0, _v, _visc, dt);
// Project velocities onto new velocities
_project(_u0, _v0, _u, _v);
// enforce new velocities to uv
_advect(1, _u, _u0, _u0, _v0, dt);
_advect(2, _v, _v0, _u0, _v0, dt);
_project(_u, _v, _u0, _v0);
// Density step
_addSource(_d, _d0, dt);
_diffuse(0, _d0, _d, _diff, dt);
_advect(0, _d, _d0, _u, _v, dt);
for (var i:int = 0; i < _size ; i++)
_d0[i] = _u0[i] = _v0[i] = 0;
}
public function getDensities():Vector.<Number> { return _d; }
public function addDensity(x:int, y:int, amount:Number):void { _d0[_index(x,y)] += amount; }
public function addVelocity(x:int, y:int, velocityX:Number, velocityY:Number):void
{
const index:int = _index(x, y);
_u0[index] += velocityX;
_v0[index] += velocityY;
}
public function setVelocity(x:int, y:int, velocityX:Number, velocityY:Number):void
{
const index:int = _index(x, y);
_u0[index] = velocityX;
_v0[index] = velocityY;
}
public function getWidth():Number { return _width; }
public function getHeight():Number { return _height; }
private function _addSource(x:Vector.<Number>, s:Vector.<Number>, dt:Number):void
{
for (var i:int = 0; i < _size; ++i)
x[i] += dt * s[i];
}
private function _diffuse(b:int, x:Vector.<Number>, x0:Vector.<Number>, diff:Number, dt:Number):void
{
const a:Number = dt * diff * _n2;
_linearlySolve(x, x0, a, b, 1 + 4 * a);
}
private function _linearlySolve(x:Vector.<Number>, x0:Vector.<Number>, a:Number, b:Number, c:Number):void
{
for (var k:int = 0; k < _solveIterations; ++k)
{
for (var j:int = 1; j <= _height; ++j)
{
var current:int = j * _rowAdder;
for (var i:int = 1 ; i <= _width; ++i)
{
++current;
x[current] =
(x0[current] +
a *
(x[current-1] +
x[current+1] +
x[current-_rowAdder] +
x[current+_rowAdder])) / c;
}
}
_setBnd(b, x);
}
}
private function _advect(b:Number, d:Vector.<Number>, d0:Vector.<Number>, u:Vector.<Number>, v:Vector.<Number>, dt:Number):void
{
const dt0:Number = dt * _width; // FIX
for (var j:int = 1; j <= _height; ++j)
{
const current:int = j * _rowAdder;
for (var i:int = 1 ; i <= _width; ++i)
{
++current;
const x:Number = Math.min(Math.max(i - dt0 * u[current], 0.5), _width + 0.5);
const y:Number = Math.min(Math.max(j - dt0 * v[current], 0.5), _height + 0.5);
const i0:int = int(x);
const i1:int = i0 + 1;
const j0:int = int(y);
const j1:int = j0+1;
const s1:Number = x - i0;
const s0:Number = 1 - s1;
const t1:Number = y - j0;
const t0:Number = 1 - t1;
d[current] = s0 *
(t0 * d0[_index(i0, j0)] +
t1 * d0[_index(i0, j1)]) + s1 *
(t0 * d0[_index(i1, j0)] +
t1 * d0[_index(i1, j1)]);
}
}
_setBnd(b, d);
}
private function _project(u:Vector.<Number>, v:Vector.<Number>, p:Vector.<Number>, div:Vector.<Number>):void
{
for (var j:int = 1; j <= _height; ++j)
{
var current:int = j * _rowAdder;
for (var i:int = 1; i <= _width; ++i)
{
++current;
div[current] = -0.5 *
(u[current+1] -
u[current-1] +
v[current+_rowAdder] -
v[current-_rowAdder]) / _width; // FIX /_N
p[current] = 0;
}
}
_setBnd(0, div);
_setBnd(0, p);
_linearlySolve(p, div, 1, 0, 4);
for (j = 1; j <= _height; ++j)
{
current = j * _rowAdder;
for (i = 1 ; i <= _width; ++i)
{
++current;
u[current] -= 0.5 * _width * (p[current+1] - p[current-1]);
v[current] -= 0.5 * _height * (p[current+_rowAdder] - p[current-_rowAdder]);
}
}
_setBnd(1, u);
_setBnd(2, v);
}
private function _setBnd(b:int, x:Vector.<Number>):void
{
for (var i:int = 1; i <= _height; ++i)
{
x[_index(0, i)] = b==1 ? -x[_index(1, i)] : x[_index(1, i)];
x[_index(_width+1, i)] = b==1 ? -x[_index(_width, i)] : x[_index(_width, i)];
}
for (i = 0; i <= _width; ++i)
{
x[_index(i, 0)] = b == 2 ? -x[_index(i, 1)] : x[_index(i, 1)];
x[_index(i,_height + 1)] = b == 2 ? -x[_index(i, _height)] : x[_index(i, _height)];
}
x[_index(0, 0)] = 0.5 * (x[_index(1,0 )]+x[_index(0 ,1)]);
x[_index(0, _height + 1)] = 0.5 * (x[_index(1, _height + 1)] + x[_index(0 , _height)]);
x[_index(_width + 1, 0)] = 0.5 * (x[_index(_width, 0)] + x[_index(_width + 1, 1)]);
x[_index(_width + 1, _height + 1)] = 0.5 * (x[_index(_width, _height + 1)] + x[_index(_width+1, _height)]);
}
private function _index(x:int, y:int):int { return x + y * _rowAdder; }
private function _createArr():Vector.<Number>
{
const arr:Vector.<Number> = new Vector.<Number>(_size, true);
for (var i:int = 0; i < _size; ++i)
arr[i] = 0;
return arr;
}
private var _solveIterations:int;
private var _diff:Number;
private var _visc:Number;
private var _width:int;
private var _height:int;
private var _n2:int;
private var _size:int;
private var _rowAdder:int;
private var _d:Vector.<Number>; // dens
private var _d0:Vector.<Number>; // prev dens
private var _u:Vector.<Number>; // vel x
private var _u0:Vector.<Number>;
private var _v:Vector.<Number>; // vel y
private var _v0:Vector.<Number>;
}
}
Merry Christmas!