This video contains five clips using SideFX Houdini's Grain Solver with an attached POP Wrangle that uses VEX to generate custom forces. Here's a quick rundown of the VEX I used for each clip (please forgive the use of a variable named oomph).

Clip One "Twin Peaks"

Here, I compare each grain's current angle to the scene's origin to the current time. If the grain's angle is within a certain distance of the current time, I calculate a vertical force based on the Gaussian function of the distance of the particle from the origin The sigma of the Gaussian function is based on the sine of the time to give a rotating and pulsating force.

float e = 2.7182818284590452353602874713527;  

float tau = $PI * 2; float width = 0.15; 

vector force = v@force;  

float particleAngle = atan2(v@P.x, v@P.z); 

float timeAngle = ($PI + @Time); 

float angleDist = abs(particleAngle - timeAngle) % tau;  

float oomph = 30;        

float angleDist2 = abs(particleAngle - timeAngle - $PI) % tau; 

if (angleDist < width || angleDist2 < width) {

    float dist =  clamp(length(@P) - 1.0, -1.0, 1.0);     

    float variance = pow(0.125 + abs(sin(@Time * 2.25)) * 0.35, 2.0);     

    float gaus = (1.0 / sqrt(tau * variance)) * pow(e, -pow(dist, 2) / (2 * variance));        

    oomph *= gaus;         

    @force = set(force.x, force.y + oomph, force.z); 

}

Clip Two: "Three Peaks"

This VEX is not dissimilar to the previous one, but rather than calculating the difference between each grain's angle and the current time, the force is based on the sine of the angle added to the current time. 

The final render is a top down view with a little depth-of-field added.

float e = 2.7182818284590452353602874713527;  

float tau = $PI * 2;  

vector force = v@force;  

float particleAngle = (atan2(v@P.x, v@P.z) * 3) + (@Time * 10); 

float oomph = 20 * max(0, sin(particleAngle));  

float dist =  clamp(length(@P) - 1.0, -1.0, 1.0);  

float variance = pow(0.125 + abs(sin(@Time * 1.5)) * 0.35, 2.0);  

float gaus = (1.0 / sqrt(tau * variance)) * pow(e, -pow(dist, 2) / (2 * variance));  

oomph *= gaus;  

@force = set(force.x, force.y + oomph, force.z);

Clip Three: "Rotating Ridges"

The upwards force is a function of each grain's angle to the scene's origin multiplied by the sine of its distance to the origin:

float tau = $PI * 2; 

float width = 0.3;  

vector force = v@force;     

float particleAngle = atan2(v@P.x, v@P.z); 

float timeAngle = ($PI + @Time);   

float dist = abs(particleAngle - timeAngle) % tau;   

if (dist < width) { 

    float oomph = 25;     

    oomph *= abs(sin((@Time * -2) + length(@P * 6)));

    @force = set(force.x, force.y + oomph, force.z); 

}

Clip Four: "Ripples 1"

The upwards force is based on the difference between each grain's distance from the scene's origin and the current time multiplied by some curlnoise. The force is also tempered by each grain's vertical position so that it decreases with height.

vector force = v@force;    

float oomph = 0;     

int dist = int(length(@P) * 20); 

int time = int(@Frame % 40);      

if (dist == time) {    

    oomph = 90; 

}     

oomph *= (1 - clamp(v@P.y, 0.0, 1.0));     

vector4 noiseSource = set(v@P.x * 0.5, v@P.y * 0.5, 

v@P.z * 0.5, @Time * 5);  

vector noise = curlnoise(noiseSource);  

oomph *= 1.0 + (0.25 * noise.x);      

@force = set(force.x, force.y + oomph, force.z);

Clip Five: "Ripples 1"

The upwards force is based on the sine of the distance of each grain from the scene's origin added to the current time.

float phase = length(@P)  * 5.0; 

phase = max(0.0, sin(phase + (@Time * -3.0)));    

float amplitude = 15 * (1 - clamp(v@P.y, 0.0, 1.0)) * (2.5 - clamp(length(@P), 0.0, 2.5));  

vector force = v@force;       

@force = set(force.x, force.y + (amplitude * phase), force.z);