Particle Advection by Reaction Diffusion in SideFX Houdini

After watching this excellent tutorial that discusses advecting particles by magnetic fields to create an animation of the sun, I was inspired to use the same technique to advect particles by fields that are a function of reaction diffusion systems.

The source for my reaction diffusion vector field is a geometry node. This field this geometry node outputs defines the concentrations of the component chemical-species that will later be used to advect the particles:

The VDB node in this network defines the name of the field, I've used rdField, and the type: a Vector Float. The Volume Wrangle defines the initial state: I want a central cube where the field's z is high surrounded by an area where the x is high:

if (@P.x > -0.1 && @P.x < 0.1 &&
@P.y > -0.1 && @P.y < 0.1 &&
@P.z > -0.1 && @P.z < 0.1) {

v@foo.x *= 0.1;
v@foo.z += rand(@P * 11) * 5;
}
else {
v@foo.z *= 0;
v@foo.x += rand(@P * 11) * 5;
}

Inside the solver, I iterate over the volume, applying the VEX to create the reaction diffusion simulation. I've done this many times in two dimensions and adding a third dimension is easy for VDBs in Houdini:

float delta = volumevoxeldiameter(0, "rdField") / sqrt(3);

vector laplacian = set(0, 0, 0);
laplacian += volumesamplev(0, "rdField", v@P + set( delta, 0, 0));
laplacian += volumesamplev(0, "rdField", v@P + set( -delta, 0, 0));
laplacian += volumesamplev(0, "rdField", v@P + set( 0, delta, 0));
laplacian += volumesamplev(0, "rdField", v@P + set( 0, -delta, 0));
laplacian += volumesamplev(0, "rdField", v@P + set( 0, 0, delta));
laplacian += volumesamplev(0, "rdField", v@P + set( 0, 0, -delta));
laplacian -= (@rdField * 6.0);

float dU = 0.189;
float dV = 0.03;
float k = 0.082;
float f = 0.0425;
float timeStep = 0.75;

float u = v@foo.x;
float v = v@foo.z;

float reactionRate = u * v * v;
float deltaU = (dU * laplacian.x) - reactionRate + f * (1 - u);
float deltaV = (dV * laplacian.z) + reactionRate - (f + k) * v;

float newU = clamp(u + deltaU * timeStep, 0, 1);
float newV = clamp(v + deltaV * timeStep, 0, 1);

v@rdField = set(newU, newU, newV);

The Laplacian is generated by sampling the neighboring voxels at each side of the current voxel. The VEX functions volumevoxeldiameter and volumesamplev make this pretty simple. Dividing the former by the square root of 3 gives the voxel side length.

To assist during development, I added a volume slice. Sure enough, after a handful of iterations, my little box in the centre of the volume has now "evolved" based on the Gray-Scott model:

The next step is to use the shelf tools to create a particle system and a DOP Network:

Houdini does have a POP Advect by Volumes node, but that expects volumes where the vector fields define velocity. My vector field defines values to which I want particles advected towards. So, the advecting happens inside a POP Wrangle node with the Input 2 set to the output of the geometry node discussed above. The VEX in this node looks at the surrounding voxels and sets the force based on the difference of opposite neighboring voxels:

float mult = 10;
float delta = volumevoxeldiameter(1, "rdField") / sqrt(3);

v@force.x += volumesamplev(1, "rdField", v@P + set( delta, 0, 0)).b * mult;
v@force.x -= volumesamplev(1, "rdField", v@P + set( -delta, 0, 0)).b * mult;

v@force.y += volumesamplev(1, "rdField", v@P + set( 0, delta, 0)).b * mult;
v@force.y -= volumesamplev(1, "rdField", v@P + set( 0, -delta, 0)).b * mult;

v@force.z += volumesamplev(1, "rdField", v@P + set( 0, 0, delta)).b * mult;
v@force.z -= volumesamplev(1, "rdField", v@P + set( 0, 0, -delta)).b * mult;

A few additional nodes keep the particles within the bounding box and add a speed limit and some drag.

The second clip takes things a little further. First of all, the Gray-Scott parameters are slightly tweaked to yield a solitons rather than a coral effect:

float dU = 0.2097;
float dV = 0.080;
float k = 0.062;
float f = 0.025;
float timeStep = 0.5;

I then extended the POP Wrangle in the DOP Network to add an additional force: not only are the particles advected by the relative nearby chemical strengths in the volume, they are also advected towards other nearby particles. The Attributes from Volume node copies the rdField into the particle's Cd and the VEX uses that:

float mult = 10;
float delta = volumevoxeldiameter(1, "rdField") / sqrt(3);
v@force.x += volumesamplev(1, "rdField", v@P + set( delta, 0, 0)).b * mult;
v@force.x -= volumesamplev(1, "rdField", v@P + set( -delta, 0, 0)).b * mult;
v@force.y += volumesamplev(1, "rdField", v@P + set( 0, delta, 0)).b * mult;
v@force.y -= volumesamplev(1, "rdField", v@P + set( 0, -delta, 0)).b * mult;
v@force.z += volumesamplev(1, "rdField", v@P + set( 0, 0, delta)).b * mult;
v@force.z -= volumesamplev(1, "rdField", v@P + set( 0, 0, -delta)).b * mult;

int point;
int count = 48;
int neighbours[] = nearpoints(0, v@P, 2, count);

foreach(point; neighbours) {
vector otherColor = point(0, "Cd", point);
vector otherPosition = point(0, "P", point);
float dist = distance(@P, otherPosition);
vector direction = normalize(otherPosition - @P);

if (dist < 0.1) {
dist = 0.5;
direction = -direction;
}

v@force += 0.04 * direction * (1 / pow(dist,2)) * otherColor.b;
}

The particles are rendered as little strands by giving them trails that feed into a Wireframe node:

The final clip in the video above is based on the Belousov-Zhabotinsky Reaction, so the solve VEX is rather different:

float alpha = 1.3;
float beta = 1.55; //1.55;
float gamma = 1.65;
float delta = volumevoxeldiameter(0, "rdField") / sqrt(3);
vector accumulator = set(0, 0, 0);

// sides...
accumulator += volumesamplev(0, "rdField", v@P + set( delta, 0, 0));
accumulator += volumesamplev(0, "rdField", v@P + set( -delta, 0, 0));
accumulator += volumesamplev(0, "rdField", v@P + set( 0, delta, 0));
accumulator += volumesamplev(0, "rdField", v@P + set( 0, -delta, 0));
accumulator += volumesamplev(0, "rdField", v@P + set( 0, 0, delta));
accumulator += volumesamplev(0, "rdField", v@P + set( 0, 0, -delta));
// top corners...
accumulator += volumesamplev(0, "rdField", v@P + set( delta, -delta, delta));
accumulator += volumesamplev(0, "rdField", v@P + set( -delta, -delta, delta));
accumulator += volumesamplev(0, "rdField", v@P + set( delta, -delta, -delta));
accumulator += volumesamplev(0, "rdField", v@P + set( -delta, -delta, -delta));
// bottom corners...
accumulator += volumesamplev(0, "rdField", v@P + set( delta, delta, delta));
accumulator += volumesamplev(0, "rdField", v@P + set( -delta, delta, delta));
accumulator += volumesamplev(0, "rdField", v@P + set( delta, delta, -delta));
accumulator += volumesamplev(0, "rdField", v@P + set( -delta, delta, -delta));
accumulator += @rdField;
vector average = ((@rdField * 2) + (accumulator / 15)) / 3;

float a = average.x + average.x * (alpha * gamma * average.y) - average.z;
float b = average.y + average.y * ((beta * average.z) - (alpha * average.x));
float c = average.z + average.z * ((gamma * average.x) - (beta * average.y));

a = clamp(a, 0, 1);
b = clamp(b, 0, 1);
c = clamp(c, 0, 1);

v@rdField = set(a, b, c);
v@rdField.r += (rand(@P * @Time) * 0.05) - 0.025;
v@rdField.g += (rand(13 * @P * @Time) * 0.05) - 0.025;
v@rdField.b += (rand(17 * @P * @Time) * 0.05) - 0.025;

In this example, the positions of the trails are smoothed out to get rid of any sharp angles:

FlexMonkey

Creative Coding in Swift

View comments