michael@0: //
michael@0: // NBody adapted from Intel's nbody benchmark.
michael@0: //
michael@0: 
michael@0: load(libdir + "util.js");
michael@0: load(libdir + "seedrandom.js");
michael@0: 
michael@0: var NBody = {
michael@0:   Constant: {
michael@0:     "deltaTime": 1,     // 0.005 in their code.
michael@0:     "epsSqr": 50,       // softening factor, when they compute, set to 50.
michael@0:     "initialVelocity": 8 // set to 0 to turn off
michael@0:   },
michael@0: 
michael@0:   init: function init(mode, numBodies) {
michael@0:     var initPos = new Array(numBodies);
michael@0:     var initVel = new Array(numBodies);
michael@0: 
michael@0:     // initialization of inputs
michael@0:     for (var i = 0; i < numBodies; i++) {
michael@0:       // [x,y,z]
michael@0:       initPos[i] = [Math.floor((Math.random()) * 40000),
michael@0:                     Math.floor((Math.random()) * 20000),
michael@0:                     Math.floor((Math.random() - .25) * 50000)];
michael@0: 
michael@0:       // [x,y,z,x,y,z]
michael@0:       initVel[i] = [(Math.random() - 0.5) * NBody.Constant.initialVelocity,
michael@0:                     (Math.random() - 0.5) * NBody.Constant.initialVelocity,
michael@0:                     (Math.random()) * NBody.Constant.initialVelocity + 10,
michael@0: 
michael@0:                     (Math.random() - 0.5) * NBody.Constant.initialVelocity,
michael@0:                     (Math.random() - 0.5) * NBody.Constant.initialVelocity,
michael@0:                     (Math.random()) * NBody.Constant.initialVelocity];
michael@0:     }
michael@0: 
michael@0:     NBody.private = {};
michael@0: 
michael@0:     if (mode === "par") {
michael@0:       NBody.private.pos = new ParallelArray(initPos);
michael@0:       NBody.private.vel = new ParallelArray(initVel);
michael@0:     } else {
michael@0:       NBody.private.pos = initPos;
michael@0:       NBody.private.vel = initVel;
michael@0:     }
michael@0: 
michael@0:     NBody.numBodies = numBodies;
michael@0:     NBody.time = 0;
michael@0:   },
michael@0: 
michael@0:   // Parallel
michael@0: 
michael@0:   tickPar: function tickPar() {
michael@0:     NBody.private.vel = new ParallelArray([NBody.numBodies], NBody.velocityPar);
michael@0:     NBody.private.pos = new ParallelArray([NBody.numBodies], NBody.positionPar);
michael@0:     NBody.time++;
michael@0:   },
michael@0: 
michael@0:   velocityPar: function velocityPar(index) {
michael@0:     var pos = NBody.private.pos;
michael@0:     var vel = NBody.private.vel;
michael@0: 
michael@0:     var deltaTime = NBody.Constant.deltaTime;
michael@0:     var epsSqr = NBody.Constant.epsSqr;
michael@0:     var time = NBody.time;
michael@0: 
michael@0:     var shape = vel.shape[0];
michael@0: 
michael@0:     var newVel;
michael@0:     var newX, newY, newZ;
michael@0:     var newX2, newY2, newZ2;
michael@0: 
michael@0:     var cX = Math.cos(time / 22) * -4200;
michael@0:     var cY = Math.sin(time / 14) * 9200;
michael@0:     var cZ = Math.sin(time / 27) * 6000;
michael@0: 
michael@0:     // pull to center
michael@0:     var maxDistance = 3400;
michael@0:     var pullStrength = .042;
michael@0: 
michael@0:     var speedLimit = 8;
michael@0: 
michael@0:     // zones
michael@0:     var zone = 400;
michael@0:     var repel = 100;
michael@0:     var align = 300;
michael@0:     var attract = 100;
michael@0: 
michael@0:     if (time < 500) {
michael@0:       speedLimit = 2000;
michael@0:       var attractPower = 100.9;
michael@0:     } else {
michael@0:       speedLimit = .2;
michael@0:       attractPower = 20.9;
michael@0:     }
michael@0: 
michael@0:     var zoneSqrd = zone * zone + zone * zone + zone * zone;
michael@0: 
michael@0:     var accX = 0, accY = 0, accZ = 0;
michael@0:     var accX2 = 0, accY2 = 0, accZ2 = 0;
michael@0:     var i;
michael@0: 
michael@0:     // define particle 1 center distance
michael@0:     var dirToCenterX = cX - pos.get(index)[0];
michael@0:     var dirToCenterY = cY - pos.get(index)[1];
michael@0:     var dirToCenterZ = cZ - pos.get(index)[2];
michael@0: 
michael@0:     var distanceSquaredTo = dirToCenterX * dirToCenterX + dirToCenterY * dirToCenterY + dirToCenterZ * dirToCenterZ;
michael@0:     var distToCenter = Math.sqrt(distanceSquaredTo);
michael@0: 
michael@0:     // orient to center
michael@0:     if (distToCenter > maxDistance) {
michael@0:       var velc = (distToCenter - maxDistance) * pullStrength;
michael@0:       if (time < 200)
michael@0:         velc = .2;
michael@0:       else velc = (distToCenter - maxDistance) * pullStrength;
michael@0: 
michael@0:       accX += (dirToCenterX / distToCenter) * velc;
michael@0:       accY += (dirToCenterY / distToCenter) * velc;
michael@0:       accZ += (dirToCenterZ / distToCenter) * velc;
michael@0:     }
michael@0: 
michael@0:     for (i = 0; i < shape; i = i + 1) {
michael@0:       var rx = pos.get(i)[0] - pos.get(index)[0];
michael@0:       var ry = pos.get(i)[1] - pos.get(index)[1];
michael@0:       var rz = pos.get(i)[2] - pos.get(index)[2];
michael@0: 
michael@0:       // make sure we are not testing the particle against its own position
michael@0:       var areSame = 0;
michael@0:       if (pos.get(i)[0] == pos.get(index)[0] && pos.get(i)[1] == pos.get(index)[1] && pos.get(i)[2] == pos.get(index)[2])
michael@0:         areSame += 1;
michael@0: 
michael@0:       var distSqrd = rx * rx + ry * ry + rz * rz;
michael@0: 
michael@0:       // cant use eqals to test, only <= or >= WTF
michael@0:       if (distSqrd < zoneSqrd && areSame <= 0) {
michael@0:         var length = Math.sqrt(distSqrd);
michael@0:         var percent = distSqrd / zoneSqrd;
michael@0: 
michael@0:         if (distSqrd < repel) {
michael@0:           var F = (repel / percent - 1) * .025;
michael@0: 
michael@0:           var normalRx = (rx / length) * F;
michael@0:           var normalRy = (ry / length) * F;
michael@0:           var normalRz = (rz / length) * F;
michael@0: 
michael@0:           accX = accX + normalRx;
michael@0:           accY = accY + normalRy;
michael@0:           accZ = accZ + normalRz;
michael@0: 
michael@0:           accX2 = accX2 - normalRx;
michael@0:           accY2 = accY2 - normalRy;
michael@0:           accZ2 = accZ2 - normalRz;
michael@0:         } else if (distSqrd < align) { //align
michael@0:           var threshDelta = align - repel;
michael@0:           var adjustedPercent = (percent - repel) / threshDelta;
michael@0:           var Q = (.5 - Math.cos(adjustedPercent * 3.14159265 * 2) * .5 + .5) * 100.9;
michael@0: 
michael@0:           // get velocity 2
michael@0:           var velX2 = vel.get(i)[4];
michael@0:           var velY2 = vel.get(i)[5];
michael@0:           var velZ2 = vel.get(i)[6];
michael@0: 
michael@0:           var velLength2 = Math.sqrt(velX2 * velX2 + velY2 * velY2 + velZ2 * velZ2);
michael@0: 
michael@0:           // normalize vel2 and multiply by factor
michael@0:           velX2 = (velX2 / velLength2) * Q;
michael@0:           velY2 = (velY2 / velLength2) * Q;
michael@0:           velZ2 = (velZ2 / velLength2) * Q;
michael@0: 
michael@0:           // get own velocity
michael@0:           var velX = vel.get(i)[0];
michael@0:           var velY = vel.get(i)[1];
michael@0:           var velZ = vel.get(i)[2];
michael@0: 
michael@0:           var velLength = Math.sqrt(velX * velX + velY * velY + velZ * velZ);
michael@0: 
michael@0:           // normalize own velocity
michael@0:           velX = (velX / velLength) * Q;
michael@0:           velY = (velY / velLength) * Q;
michael@0:           velZ = (velZ / velLength) * Q;
michael@0: 
michael@0:           accX += velX2;
michael@0:           accY += velY2;
michael@0:           accZ += velZ2;
michael@0: 
michael@0:           accX2 += velX;
michael@0:           accY2 += velY;
michael@0:           accZ2 += velZ;
michael@0:         }
michael@0: 
michael@0:         if (distSqrd > attract) { // attract
michael@0:           var threshDelta2 = 1 - attract;
michael@0:           var adjustedPercent2 = (percent - attract) / threshDelta2;
michael@0:           var C = (1 - (Math.cos(adjustedPercent2 * 3.14159265 * 2) * 0.5 + 0.5)) * attractPower;
michael@0: 
michael@0:           // normalize the distance vector
michael@0:           var dx = (rx / (length)) * C;
michael@0:           var dy = (ry / (length)) * C;
michael@0:           var dz = (rz / (length)) * C;
michael@0: 
michael@0:           accX += dx;
michael@0:           accY += dy;
michael@0:           accZ += dz;
michael@0: 
michael@0:           accX2 -= dx;
michael@0:           accY2 -= dy;
michael@0:           accZ2 -= dz;
michael@0:         }
michael@0:       }
michael@0:     }
michael@0: 
michael@0:     // Speed limits
michael@0:     if (time > 500) {
michael@0:       var accSquared = accX * accX + accY * accY + accZ * accZ;
michael@0:       if (accSquared > speedLimit) {
michael@0:         accX = accX * .015;
michael@0:         accY = accY * .015;
michael@0:         accZ = accZ * .015;
michael@0:       }
michael@0: 
michael@0:       var accSquared2 = accX2 * accX2 + accY2 * accY2 + accZ2 * accZ2;
michael@0:       if (accSquared2 > speedLimit) {
michael@0:         accX2 = accX2 * .015;
michael@0:         accY2 = accY2 * .015;
michael@0:         accZ2 = accZ2 * .015;
michael@0:       }
michael@0:     }
michael@0: 
michael@0:     // Caclulate new velocity
michael@0:     newX = (vel.get(index)[0]) + accX;
michael@0:     newY = (vel.get(index)[1]) + accY;
michael@0:     newZ = (vel.get(index)[2]) + accZ;
michael@0: 
michael@0:     newX2 = (vel.get(index)[3]) + accX2;
michael@0:     newY2 = (vel.get(index)[4]) + accY2;
michael@0:     newZ2 = (vel.get(index)[5]) + accZ2;
michael@0: 
michael@0:     if (time < 500) {
michael@0:       var acs = newX2 * newX2 + newY2 * newY2 + newZ2 * newZ2;
michael@0:       if (acs > speedLimit) {
michael@0:         newX2 = newX2 * .15;
michael@0:         newY2 = newY2 * .15;
michael@0:         newZ2 = newZ2 * .15;
michael@0:       }
michael@0: 
michael@0:       var acs2 = newX * newX + newY * newY + newZ * newZ;
michael@0:       if (acs2 > speedLimit) {
michael@0:         newX = newX * .15;
michael@0:         newY = newY * .15;
michael@0:         newZ = newZ * .15;
michael@0:       }
michael@0:     }
michael@0: 
michael@0:     return [newX, newY, newZ, newX2, newY2, newZ2];
michael@0:   },
michael@0: 
michael@0:   positionPar: function positionPar(index) {
michael@0:     var vel = NBody.private.vel;
michael@0:     var pos = NBody.private.pos;
michael@0: 
michael@0:     var x = 0;
michael@0:     var y = 0;
michael@0:     var z = 0;
michael@0: 
michael@0:     var velX = vel.get(index)[0];
michael@0:     var velY = vel.get(index)[1];
michael@0:     var velZ = vel.get(index)[2];
michael@0: 
michael@0:     var velX2 = vel.get(index)[3];
michael@0:     var velY2 = vel.get(index)[4];
michael@0:     var velZ2 = vel.get(index)[5];
michael@0: 
michael@0:     var netVelX = (velX - velX2);
michael@0:     var netVelY = (velY - velY2);
michael@0:     var netVelZ = (velZ - velZ2);
michael@0: 
michael@0:     x = pos.get(index)[0] + (netVelX);
michael@0:     y = pos.get(index)[1] + (netVelY);
michael@0:     z = pos.get(index)[2] + (netVelZ);
michael@0: 
michael@0:     return [x, y, z];
michael@0:   },
michael@0: 
michael@0:   // Sequential
michael@0: 
michael@0:   tickSeq: function tickSeq() {
michael@0:     var numBodies = NBody.numBodies;
michael@0:     var newVel = new Array(numBodies);
michael@0:     var newPos = new Array(numBodies);
michael@0: 
michael@0:     for (var i = 0; i < numBodies; i++)
michael@0:       newVel[i] = NBody.velocitySeq(i);
michael@0: 
michael@0:     for (var i = 0; i < numBodies; i++)
michael@0:       newPos[i] = NBody.positionSeq(i);
michael@0: 
michael@0:     NBody.private.vel = newVel;
michael@0:     NBody.private.pos = newPos;
michael@0: 
michael@0:     NBody.time++;
michael@0:   },
michael@0: 
michael@0:   velocitySeq: function velocitySeq(index) {
michael@0:     var vel = NBody.private.vel;
michael@0:     var pos = NBody.private.pos;
michael@0: 
michael@0:     var deltaTime = NBody.Constant.deltaTime;
michael@0:     var epsSqr = NBody.Constant.epsSqr;
michael@0:     var time = NBody.time;
michael@0: 
michael@0:     var shape = pos.length;
michael@0: 
michael@0:     var newVel;
michael@0:     var newX, newY, newZ;
michael@0:     var newX2, newY2, newZ2;
michael@0: 
michael@0:     var cX = Math.cos(time / 22) * -4200;
michael@0:     var cY = Math.sin(time / 14) * 9200;
michael@0:     var cZ = Math.sin(time / 27) * 6000;
michael@0: 
michael@0:     // pull to center
michael@0:     var maxDistance = 3400;
michael@0:     var pullStrength = .042;
michael@0: 
michael@0:     var speedLimit = 8;
michael@0: 
michael@0:     // zones
michael@0:     var zone = 400;
michael@0:     var repel = 100;
michael@0:     var align = 300;
michael@0:     var attract = 100;
michael@0: 
michael@0: 
michael@0:     if (time < 500) {
michael@0:       speedLimit = 2000;
michael@0:       var attractPower = 100.9;
michael@0:     } else {
michael@0:       speedLimit = .2;
michael@0:       attractPower = 20.9;
michael@0:     }
michael@0: 
michael@0:     var zoneSqrd = zone * zone + zone * zone + zone * zone;
michael@0: 
michael@0:     var accX = 0, accY = 0, accZ = 0;
michael@0:     var accX2 = 0, accY2 = 0, accZ2 = 0;
michael@0:     var i;
michael@0: 
michael@0:     // define particle 1 center distance
michael@0:     var dirToCenterX = cX - pos[index][0];
michael@0:     var dirToCenterY = cY - pos[index][1];
michael@0:     var dirToCenterZ = cZ - pos[index][2];
michael@0: 
michael@0:     var distanceSquaredTo = dirToCenterX * dirToCenterX + dirToCenterY * dirToCenterY + dirToCenterZ * dirToCenterZ;
michael@0:     var distToCenter = Math.sqrt(distanceSquaredTo);
michael@0: 
michael@0:     // orient to center
michael@0:     if (distToCenter > maxDistance) {
michael@0:       var velc = (distToCenter - maxDistance) * pullStrength;
michael@0:       if (time < 200)
michael@0:         velc = .2;
michael@0:       else velc = (distToCenter - maxDistance) * pullStrength;
michael@0: 
michael@0:       accX += (dirToCenterX / distToCenter) * velc;
michael@0:       accY += (dirToCenterY / distToCenter) * velc;
michael@0:       accZ += (dirToCenterZ / distToCenter) * velc;
michael@0:     }
michael@0: 
michael@0:     for (i = 0; i < shape; i = i + 1) {
michael@0:       var rx = pos[i][0] - pos[index][0];
michael@0:       var ry = pos[i][1] - pos[index][1];
michael@0:       var rz = pos[i][2] - pos[index][2];
michael@0: 
michael@0:       var areSame = 0;
michael@0:       if (pos[i][0] == pos[index][0] && pos[i][1] == pos[index][1] && pos[i][2] == pos[index][2])
michael@0:         areSame += 1;
michael@0: 
michael@0:       var distSqrd = rx * rx + ry * ry + rz * rz;
michael@0: 
michael@0:       // cant use eqals to test, only <= or >= WTF
michael@0:       if (distSqrd < zoneSqrd && areSame <= 0) {
michael@0:         var length = Math.sqrt(distSqrd);
michael@0:         var percent = distSqrd / zoneSqrd;
michael@0: 
michael@0: 
michael@0:         if (distSqrd < repel) {
michael@0:           var F = (repel / percent - 1) * .025;
michael@0: 
michael@0:           var normalRx = (rx / length) * F;
michael@0:           var normalRy = (ry / length) * F;
michael@0:           var normalRz = (rz / length) * F;
michael@0: 
michael@0:           accX = accX + normalRx;
michael@0:           accY = accY + normalRy;
michael@0:           accZ = accZ + normalRz;
michael@0: 
michael@0:           accX2 = accX2 - normalRx;
michael@0:           accY2 = accY2 - normalRy;
michael@0:           accZ2 = accZ2 - normalRz;
michael@0:         } else if (distSqrd < align) { //align
michael@0:           var threshDelta = align - repel;
michael@0:           var adjustedPercent = (percent - repel) / threshDelta;
michael@0:           var Q = (.5 - Math.cos(adjustedPercent * 3.14159265 * 2) * .5 + .5) * 100;
michael@0: 
michael@0:           // get velocity 2
michael@0:           var velX2 = vel[i][3];
michael@0:           var velY2 = vel[i][4];
michael@0:           var velZ2 = vel[i][5];
michael@0: 
michael@0:           var velLength2 = Math.sqrt(velX2 * velX2 + velY2 * velY2 + velZ2 * velZ2);
michael@0: 
michael@0:           // normalize vel2 and multiply by factor
michael@0:           velX2 = (velX2 / velLength2) * Q;
michael@0:           velY2 = (velY2 / velLength2) * Q;
michael@0:           velZ2 = (velZ2 / velLength2) * Q;
michael@0: 
michael@0:           // get own velocity
michael@0:           var velX = vel[i][0];
michael@0:           var velY = vel[i][1];
michael@0:           var velZ = vel[i][2];
michael@0: 
michael@0:           var velLength = Math.sqrt(velX * velX + velY * velY + velZ * velZ);
michael@0: 
michael@0:           // normalize own velocity
michael@0:           velX = (velX / velLength) * Q;
michael@0:           velY = (velY / velLength) * Q;
michael@0:           velZ = (velZ / velLength) * Q;
michael@0: 
michael@0:           accX += velX2;
michael@0:           accY += velY2;
michael@0:           accZ += velZ2;
michael@0: 
michael@0:           accX2 += velX;
michael@0:           accY2 += velY;
michael@0:           accZ2 += velZ;
michael@0:         }
michael@0: 
michael@0:         if (distSqrd > attract) {        //attract
michael@0:           var threshDelta2 = 1 - align;
michael@0:           var adjustedPercent2 = (percent - align) / threshDelta2;
michael@0:           var C = (1 - (Math.cos(adjustedPercent2 * 3.14159265 * 2) * 0.5 + 0.5)) * attractPower;
michael@0: 
michael@0:           // normalize the distance vector
michael@0:           var dx = (rx / (length)) * C;
michael@0:           var dy = (ry / (length)) * C;
michael@0:           var dz = (rz / (length)) * C;
michael@0: 
michael@0:           debug = 1.1;
michael@0: 
michael@0:           accX += dx;
michael@0:           accY += dy;
michael@0:           accZ += dz;
michael@0: 
michael@0:           accX2 -= dx;
michael@0:           accY2 -= dy;
michael@0:           accZ2 -= dz;
michael@0:         }
michael@0:       }
michael@0:     }
michael@0: 
michael@0:     // enforce speed limits
michael@0:     if (time > 500) {
michael@0:       var accSquared = accX * accX + accY * accY + accZ * accZ;
michael@0:       if (accSquared > speedLimit) {
michael@0:         accX = accX * .015;
michael@0:         accY = accY * .015;
michael@0:         accZ = accZ * .015;
michael@0:       }
michael@0: 
michael@0:       var accSquared2 = accX2 * accX2 + accY2 * accY2 + accZ2 * accZ2;
michael@0:       if (accSquared2 > speedLimit) {
michael@0:         accX2 = accX2 * .015;
michael@0:         accY2 = accY2 * .015;
michael@0:         accZ2 = accZ2 * .015;
michael@0:       }
michael@0:     }
michael@0: 
michael@0:     // Caclulate new velocity
michael@0:     newX = vel[index][0] + accX;
michael@0:     newY = vel[index][1] + accY;
michael@0:     newZ = vel[index][2] + accZ;
michael@0: 
michael@0:     newX2 = vel[index][3] + accX2;
michael@0:     newY2 = vel[index][4] + accY2;
michael@0:     newZ2 = vel[index][5] + accZ2;
michael@0: 
michael@0:     if (time < 500) {
michael@0:       var acs = newX2 * newX2 + newY2 * newY2 + newZ2 * newZ2;
michael@0:       if (acs > speedLimit) {
michael@0:         newX2 = newX2 * .15;
michael@0:         newY2 = newY2 * .15;
michael@0:         newZ2 = newZ2 * .15;
michael@0:       }
michael@0: 
michael@0:       var acs2 = newX * newX + newY * newY + newZ * newZ;
michael@0:       if (acs2 > speedLimit) {
michael@0:         newX = newX * .15;
michael@0:         newY = newY * .15;
michael@0:         newZ = newZ * .15;
michael@0:       }
michael@0:     }
michael@0: 
michael@0:     return [newX, newY, newZ, newX2, newY2, newZ2];
michael@0:   },
michael@0: 
michael@0:   positionSeq: function positionSeq(index) {
michael@0:     var vel = NBody.private.vel;
michael@0:     var pos = NBody.private.pos;
michael@0: 
michael@0:     var x = 0;
michael@0:     var y = 0;
michael@0:     var z = 0;
michael@0: 
michael@0:     var velX = vel[index][0];
michael@0:     var velY = vel[index][1];
michael@0:     var velZ = vel[index][2];
michael@0: 
michael@0:     var velX2 = vel[index][3];
michael@0:     var velY2 = vel[index][4];
michael@0:     var velZ2 = vel[index][5];
michael@0: 
michael@0:     var netX = (velX - velX2);
michael@0:     var netY = (velY - velY2);
michael@0:     var netZ = (velZ - velZ2);
michael@0: 
michael@0:     x = pos[index][0] + netX;
michael@0:     y = pos[index][1] + netY;
michael@0:     z = pos[index][2] + netZ;
michael@0: 
michael@0:     return [x, y, z];
michael@0:   }
michael@0: };
michael@0: 
michael@0: function emulateNBody(mode, numBodies, ticks) {
michael@0:   NBody.init(mode, numBodies);
michael@0:   for (var i = 0; i < ticks; i++) {
michael@0:     var start = Date.now();
michael@0:     if (mode === "par")
michael@0:       NBody.tickPar();
michael@0:     else
michael@0:       NBody.tickSeq();
michael@0:     //print(NBody.private.pos);
michael@0:     print(mode + " bodies=" + numBodies + " tick=" + (i+1) + "/" + ticks + ": " + (Date.now() - start) + " ms");
michael@0:   }
michael@0: }
michael@0: 
michael@0: // Using 4000 bodies going off Rick's comment as 4000 being a typical workload.
michael@0: const NUMBODIES = 4000;
michael@0: const TICKS = 10;
michael@0: 
michael@0: Math.seedrandom("seed");
michael@0: 
michael@0: benchmark("NBODY", 1, DEFAULT_MEASURE,
michael@0:           function () { emulateNBody("seq", NUMBODIES, TICKS); },
michael@0:           function () { emulateNBody("par", NUMBODIES, TICKS); });