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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4640 - (download) (annotate)
Tue Sep 27 20:54:47 2016 UTC (2 years, 10 months ago) by glk
File size: 2972 byte(s)
initial result of svn export --username anonsvn --password=anonsvn https://svn.smlnj-gforge.cs.uchicago.edu/svn/diderot/branches/vis15/src/tests/
// unit-circle
//
// Demo of distributing particles on the unit circle
//
vec2{} initPosns = load("-");
int numOfParticles = length(initPosns);
input int iterMax = 500;  // maximum number of steps to run
input real rad = 0.2;       // particle radius
real hh0 = 1000; // initial step size, can be big
int iter = 1;             // which iteration we're on
real stepMax = 2*rad;  // limit on distance to travel per iter

// 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(vec2 x) = φ(|x|/rad);
function vec2 frc(vec2 x) = φ'(|x|/rad) * (1/rad) * x/|x|; // chain rule

strand Particle (int ii, vec2 pos0) {
   vec2 pos = pos0/|pos0|;
   real hh = hh0;
   output vec2 outPos = pos;
   vec2 dpos = [0,0];
   int id = ii;
   real energyLast = ii;     // HEY: can do convergence test based on variance of energies
   vec2 force = [0,0];   // or can test convergence based on sum of |force|
   update {
      // compute energy and forces on us
      if (!( pos[0] <= 0 || pos[0] > 0 )) {
         print(ii"HEY pos = ", pos, "\n");
      }
      energyLast = 0;
      force = [0,0];
      int ncount = 0;
      foreach (Particle p_j in sphere(rad)) {
         ncount += 1;
         vec2 r_ij = p_j.pos - pos;
         energyLast += enr(r_ij);
         force += frc(r_ij);
      }
      vec2 norm = normalize(pos);
      if (ncount < 2) {
        if (iter % 3 == 0) {
           vec2 npos = pos + 0.2*rad*[norm[1],-norm[0]];
           new Particle(0, npos);
        }
      }
      // update position based on force
      force = (identity[2] - norm⊗norm)•force;  // project force onto tangent surface
      dpos = hh*force;
      if (|dpos| > stepMax) {
         // decrease hh by factor by which step was too big
         hh *= stepMax/|dpos|;
         // and find smaller step
         dpos = hh*force; // == stepMax*normalize(force);
      }
      vec2 posLast = pos;
      // take step and re-apply constraint
      pos = normalize(pos + dpos);
      real energy = 0;
      foreach (Particle p_j in sphere(rad)) {
         energy += enr(p_j.pos - pos);
      }
      if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
         hh *= 0.5; // backtrack
         pos = posLast;
      } else {
         hh *= 1.1;
      }
      outPos = pos;
   }
}

global {
   /*
   real mvmt = mean { |P.dpos|/rad | P in Particle.all};
   if (mvmt < 0.0001) {
      print("stabilized at iter ", iter, "\n");
      stabilize;
   }
   */
   iter+=1;
   if (iter > iterMax)  {
      print("numActive() = ", numActive(), "\n");
      stabilize; // HEY really want "die" here: didn't converge
   }
//  real var = variance{P.energy | P in Particle.all};
//  real var = sum{(avg-P.energy)^2 | P in Particle.all};
}

initially {Particle(ii, initPosns[ii])
           | ii in 1 .. numOfParticles };

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