64K Particles

This is basically snowplow, but cranked up to 256x256 (64K) grid of particles, and the particle diameter and physics toned down for stability. You can't see the particles clearly at the 512x512 F texture resolution, but this number of particles exceeds anything an O(n²) simulation could handle interactively.

Er, JavaScript or WebGL doesn't seem to be running, so basically all you're going to see is the bare code...

// JavaScript setup: runs once at startup

// Create 2D texture to store the particles
sim.nP=256*256; // number of particles (size of texture)
sim.P=new gVector1D(sim.nP); // current particle positions
sim.Pnext=new gVector1D(sim.nP); // positions in next frame
sim.V=new gVector1D(sim.nP); // current particle velocities
sim.Vnext=new gVector1D(sim.nP); // velocities in next frame

// For debugging, show P and V onscreen
sim.P.debugTo(scene,10.0,-6.0,10.0,6.0); // P left
sim.V.debugTo(scene,10.0,+6.0,10.0,6.0); // V right

/* Return a renderable object used to show
   the srcP texture's positions onscreen. 
   
 scale affects the center locations.
 size affects the diameter of each splat.
 targetcoords is "xyz" for X-Y plane; "xzy" for X-Z plane splats.
*/
function makeSplatter(srcP,scale,size,targetcoords)
{
// Extract simulated positions in a vertex shader.
//   Basically we look up each vertex's 3D position from P.
var Pgeom=new THREE.Geometry();

// Loop over the 2D pixel centers
var dx=1.0/srcP.tex.width;
var dy=1.0/srcP.tex.height;
// Cen is the location of the particle
for (var cx=0.5*dx;cx<1.0;cx+=dx)
for (var cy=0.5*dy;cy<1.0;cy+=dy) {
  var cen=new vec3(cx,cy,0.0);
  var idx=[];
  
  // Del is deviation around particle center
  for (var l=0.0;l<=1.0;l++)
  for (var t=0.0;t<=1.0;t++) {
    idx.push(Pgeom.vertices.length);
    var del=new vec3(l,t,0.0);
    // HACK: pass del as position
    Pgeom.vertices.push(del); 
  }
  Pgeom.faces.push(new THREE.Face4(
         idx[0],idx[1],idx[3],idx[2],
         cen // HACK: pass cen as normal
  ));
}

// Turn off view and backface culling
Pgeom.boundingSphere = new THREE.Sphere( 
         new THREE.Vector3(), 100.0 /* huge! */ );
Pgeom.uniforms= {
        readP: {type: "t", value: srcP.tex }
};
var PgeomShader=new THREE.ShaderMaterial( {
     side: THREE.DoubleSide,
     blending: THREE.AdditiveBlending,
     transparent:true, // makes blending work (?!)
     depthWrite:false, // don't mess with Z buffer
     uniforms:Pgeom.uniforms,

/* The vertex shader fetches our particle's position
   from sim.P.tex */
     vertexShader:
"varying vec2 cen,del;\n"+
"uniform sampler2D readP;\n"+
"void main(void) {\n"+
"   cen=vec2(normal); \n"+
"   del=vec2(position); \n"+
"   vec3 P=vec3(texture2D(readP,cen));\n"+
"   P.xy+="+size+"*(vec2(del)-vec2(0.5));\n"+
"   P."+targetcoords+"=P.xyz;\n"+
"   P*="+scale+"; // scale up!  \n"+
"   gl_Position=projectionMatrix * modelViewMatrix * vec4( P, 1.0 );\n"+
"}\n"
,
/* The fragment shader draws a soft smooth curve */
    fragmentShader:
"varying vec2 cen,del;\n"+
"void main(void) {\n"+
"  //gl_FragColor=vec4(1,0,0,1); return; //debug \n"+
"  float radius=length(del-vec2(0.5,0.5));\n"+
"  if (radius>0.5) discard;\n"+
"  float alpha=0.8*(1.0-2.0*radius);\n"+
//"  float alpha=(cos(2.0*3.14159*radius)+1.0)*0.5;\n"+
"  // alpha=fract(4.0*alpha);\n"+
"  gl_FragColor=vec4(alpha); \n"+
"}\n"
} );

var P3D=new THREE.Mesh(Pgeom, PgeomShader);
P3D.castShadow=true; // why dis no work?
P3D.uniforms=Pgeom.uniforms;
return P3D;
}

// Draw the particle locations out in 3D (not needed)
// scene.add(makeSplatter(sim.P,10.0001,0.01,"xzy"));

// Draw the particle locations offscreen.
//   The shader defines the boundary conditions.
sim.F=new gTexture2D(512,512,
    "varying vec2 texcoords;\n"+
    "uniform float time;\n"+
    "void main() { \n"+
    "  vec2 tc=texcoords;\n"+
    "  float d=min(min(tc.x,1.0-tc.x),\n"+
    "              min(tc.y,1.0-tc.y));\n"+
    "  vec2 sphere=vec2(sin(0.1*time),0.5);\n"+
    "  d=min(d,length(tc-sphere));\n"+
    "  d=0.8*clamp(1.0-20.0*d,0.0,1.0);\n"+
    "  gl_FragColor=vec4(d,d,d,d);\n"+
    "}\n");
sim.F.debugTo(scene,10.0,0.0,8.0,6.0);
var diam=0.01;
sim.F.splat=makeSplatter(sim.P,1.000001,diam,"xyz")
sim.F.scene.add(sim.F.splat);

// Ground
var groundTex=THREE.ImageUtils.loadTexture( "textures/checkerboard_noisy.jpg" );
groundTex.wrapS=groundTex.wrapT=THREE.RepeatWrapping;
groundTex.repeat=new vec2(25.0,25.0);

sim.ground = new THREE.Mesh(
  new THREE.CubeGeometry(50,50,0.0001),
  new THREE.MeshLambertMaterial(
    {color: 0xffccaa, opacity: 1, map:groundTex})
);
sim.ground.position.z=-0.001;
sim.ground.receiveShadow=true;
scene.add(sim.ground);

// Sun-like spotlight
var l=new THREE.SpotLight();
sim.light = l;
l.position.set( -50, -100, 100 );
l.castShadow=true;
l.shadowCameraNear = 50; 
l.shadowCameraFar = 500; 
l.shadowCameraFov = 10; // degrees field of view
scene.add(l);

// First time, compile shaders and initialize buffers.
if (!sim.fnV) {
  sim.fnV=new gFunction1D(PixAnvil.loadTab('V'));
  sim.fnP=new gFunction1D(PixAnvil.loadTab('P'));
  
  // Initialize particle positions to a nice grid
  new gFunction1D(
    "vec4 run(float i) {\n"+
    "  vec4 P=fract( vec4(i/256.0,i/(256.0*256.0),1.0,1.0) );\n"+
    "  return P; // -0.5;\n"+
    "}\n"
  ).run(sim.P);

// Initialize velocities to semi-random values
  new gFunction1D(
   "vec4 run(float index) {\n"+
   "  vec4 P=Ptex(index);\n"+
   "  float f=10.0; // frequency\n"+
   "  float a=0.01; // amplitude\n"+
   "  vec4 V=vec4(a*sin(f*P.x),a*sin(f*P.y),1.0,0.0);\n"+
   "  V=normalize(V)*0.01;\n"+
   "  V.xy=0.1*(P.xy-vec2(0.5));\n"+
   "  return V;\n"+
   "}\n"
  ) .set("Ptex",sim.P) .run(sim.V);
}

var ground_z=0.01; // up, to prevent Z buffer fighting

// Substepping loop:
if (!sim.target_t) {
   sim.target_t=sim.current_t=0;
}
if (lib.dt<0.2) sim.target_t+=lib.dt;

var dt=0.02; // simulation timestep size
while (sim.current_t<sim.target_t) { // substep loop
  sim.current_t+=dt;

// Splat positions into force field
  sim.F.splat.uniforms.readP.value=sim.P.tex;
  sim.F.uniforms.time.value=sim.current_t;
  sim.F.render();

// This GLSL code updates the V velocities
  sim.fnV // set uniforms and compute V
    .set("dt",dt)
    .set("gravity",new vec3(0,0,-9.8))
    .set("ground_z",ground_z)
    .set("Ptex",sim.P) // read from old P
    .set("Vtex",sim.V) // read from old V
    .set("Ftex",sim.F.tex) // read from old F
    .run(sim.Vnext); // write to Vnext
  sim.Vnext.swap(sim.V); // ping-pong buffers

// This GLSL code updates the P positions
  sim.fnP // set uniforms and compute P
    .set("dt",dt)
    .set("ground_z",ground_z)
    .set("Ptex",sim.P) // read from (old) P
    .set("Vtex",sim.V) // read from (new) V
    .run(sim.Pnext); // write to Pnext
  sim.Pnext.swap(sim.P); // ping-pong buffers

} // end substep loop

// Camera control:

// Update camera coordinate system
var s=camera;
if (!s.X) { // startup: create initial coordinates
    s.X=new vec3(1,0,0);
    s.Y=new vec3(0,0,1);
    s.Z=new vec3(0,-1,-0.2); // camera Z is world Y
    s.P=new vec3(0,-5,1.3); // initial location
}

// Move camera via keyboard
var move=new vec3(0,0,0); // sums current frame motion

// X control via A and D
if (lib.key['a']) move.pe(new vec3(-1,0,0));
if (lib.key['d']) move.pe(new vec3(+1,0,0));
// Y control via W and S
if (lib.key['w']) move.pe(new vec3(0,0,-1));
if (lib.key['s']) move.pe(new vec3(0,0,+1));
// Z control via Q and Z
if (lib.key['q']) move.pe(new vec3(0,+1,0));
if (lib.key['z']) move.pe(new vec3(0,-1,0));

move.te(5.0*lib.dt); // meters/second motion rate
s.P.pe(s.X.t(move.x).p(s.Y.t(move.y)).p(s.Z.t(move.z)));

// Rotate camera via mouse
var speed=0.01; // radians per mouse pixel
if (lib.mouseleft) { // move Z with mouse
    s.Z.pe(s.X.t(-lib.mousedx*speed).p(
           s.Y.t( lib.mousedy*speed)));
}

// Keep level: make sure X is horizontal.
s.X.z=0.0;
s.Y.crossVectors(s.Z,s.X.normalize());

// Orthonormalize
s.X.crossVectors(s.Y,s.Z).normalize();
s.Y.crossVectors(s.Z,s.X).normalize();
s.Z.normalize();

// Write coordinate system into matrix
s.matrixAutoUpdate=false; // don't trash
s.matrixWorldNeedsUpdate=true; // show
  
var m=s.matrix; // the camera's matrix

// Utility function: set a matrix column
function setCol(m,col,vec) {
  m.elements[4*col+0]=vec.x;
  m.elements[4*col+1]=vec.y;
  m.elements[4*col+2]=vec.z;
}

setCol(m,0,s.X);
  setCol(m,1,s.Y);
  setCol(m,2,s.Z);
  setCol(m,3,s.P); // position from sim