Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /tests/lamont-tests/unit-sphere.diderot
ViewVC logotype

View of /tests/lamont-tests/unit-sphere.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4666 - (download) (annotate)
Thu Sep 29 18:28:48 2016 UTC (2 years, 9 months ago) by glk
File size: 4803 byte(s)
tweaking
// unit-sphere
//
// Demo of distributing particles on the unit sphere
//
input vec3{} initPos ("initial positions for all particles") = load("vec3.nrrd");
input int iterMax ("maximum number of iterations to run") = 100;
input real rad ("radius of particle's energy support") = 0.1;
input int pcp ("periodicity of population control (if non-zero)") = 5;
input real mvmtEps ("convergence threshold on particle system movement") = 0.001;
input bool noDie ("all strands stabilize instead of die at iterMax") = false;

real hhInit = 10;       // initial step size, can be big
real newDist = 0.9*rad; // how far away to birth new particles
real stepMax = 2*rad;   // limit on distance to travel per iter

int iter = 0;           // which iteration we're on

// inter-particle energy, and its derivative
//function real  φ(real r) = (1 - r)^4;
//function real φ'(real r) = -4*(1 - r)^3;
function real  φ(real r)  =  (1/r)*(1-r)^3;
function real φ'(real r) = 3 - 1/(r^2) - 2*r;
function real enr(vec3 x) = φ(|x|/rad);
function vec3 frc(vec3 x) = φ'(|x|/rad) * (1/rad) * x/|x|; // chain rule

// returns a non-zero vector perpendicular to given non-zero vector v
function vec3 perp3(vec3 v) {
   int c = 0;
   if (|v[0]| < |v[1]|) {
      c = 1;
   }
   // can't say v[c] because tensors can only be indexed by constants
   if (|v[1] if 1==c else v[0]| < |v[2]|) {
      c = 2;
   }
   // now c is index of v component with largest absolute value
   vec3 ret = ([v[1] - v[2], -v[0], v[0]] if (c == 0) else
               [-v[1], v[0] - v[2], v[1]] if (c == 1) else
               [-v[2], v[2], v[0] - v[1]]);
   return ret;
}

strand Particle (int parent_, int ii, vec3 pos0, real hh0) {
   int parent = parent_;
   vec3 pos = pos0/|pos0|;
   real hh = hh0;
   output vec3 outPos = pos;
   vec3 delta = [0,0,0];
   int id = ii;
   int born = 0;
   real energyLast = ii;     // HEY: can do convergence test based on variance of energies
   vec3 force = [0,0,0];   // or can test convergence based on sum of |force|
   bool first = true;
   int ncount = 0;
   update {
      if (first) {
         //print("(iter ", iter, ") hello from ", ii, "\n");
         first = false; // HEY remove me
      }
      // compute energy and forces on us
      energyLast = 0;
      force = [0,0,0];
      foreach (Particle p_j in sphere(rad)) {
         ncount += 1;
         vec3 r_ij = p_j.pos - pos;
         energyLast += enr(r_ij);
         force += frc(r_ij);
      }
      vec3 norm = normalize(pos); // surface normal (for unit circle)
      if (ncount == 0 && (pcp > 0 && 0 == iter % pcp)) {
         vec3 npos = pos + newDist*rad*normalize(perp3(norm));
         new Particle(ii, born, npos, hh);
         born += 1;
         continue;
      }
      // else update position based on force from neighbors.
      // first, project force onto tangent surface
      force = (identity[3] - norm⊗norm)•force;
      delta = hh*force;
      if (|delta| > stepMax) {
         // decrease hh by factor by which step was too big
         hh *= stepMax/|delta|;
         // and find smaller step
         delta = hh*force; // == stepMax*normalize(force);
      }
      vec3 posLast = pos;
      // take step and re-apply constraint
      pos = normalize(pos + delta);
      real energy = 0;
      ncount = 0;
      foreach (Particle p_j in sphere(rad)) {
         energy += enr(p_j.pos - pos);
         ncount += 1;
      }
      if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
         hh *= 0.5; // backtrack
         pos = posLast;
      } else {
         hh *= 1.1;
         /* putting new particle creation here avoids possibility of
            starting multiple particles at same position while not
            moving due to backtracking */
         if (ncount < 5 && (pcp > 0 && 0 == iter % pcp)) {
            vec3 npos = pos + newDist*rad*normalize(force);
            new Particle(ii, born, npos, hh);
            born += 1;
         }
      }
      outPos = pos;
   }
}

global {
   //print("(iter ", iter, ") hello from global\n");
   real mvmt = max { |P.delta|/rad | P in Particle.all};
   if (mvmt < mvmtEps          // there is very little movement
       && (pcp == 0            // and, either there is no population control
          || iter > 2*pcp)) {  // or new particles should be moving by now
      print("Stabilizing with movement ", mvmt, " < ", mvmtEps, " (iter ", iter, ")\n");
      stabilize;
   }
   iter += 1;
   if (iter >= iterMax && iterMax > 0)  {
      print("Stabilizing (due to noDie)" if noDie else "STOPPING",
            " at iter ", iter, " >= ", iterMax, " (movement ", mvmt, ")\n");
      if (noDie) {
         stabilize;
      } else {
         die;
      }
   }
}

initially {Particle(-1, ii, initPos[ii], hhInit)
           | ii in 0 .. length(initPos)-1 };


root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0