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

SCM Repository

[diderot] View of /tests/vis15-bugs/circle.diderot
ViewVC logotype

View of /tests/vis15-bugs/circle.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4729 - (download) (annotate)
Tue Oct 11 21:32:41 2016 UTC (2 years, 10 months ago) by glk
File size: 3880 byte(s)
adding script for (fixed) bug with saving first snapshot
/*
running this under valgrind produces:
Syscall param write(buf) points to uninitialised byte(s)" that seems to be about the memory in the nrrd not being set
due to saving out the first (post-initialization, before first iteration) snapshot

*/

input vec2{} ipos ("initial positions for all particles") = load("vec2.nrrd");
input real rad ("radius of particle's potential energy support") = 0.25;
input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;

real hhInit = 1;        // initial step size
real stepMax = rad;     // limit on distance to travel per iter
int iter = 0;           // which iteration we're on

// phi(0) and phi'(0) are bounded
function real  phi(real r) = (1 - r)^4;
function real phi'(real r) = -4*(1 - r)^3;

function real enr(vec2 x) = phi(|x|/rad);
function vec2 frc(vec2 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule

/* The particle is initialized at position pos0 */
strand particle (vec2 pos0) {
   /* These strand variables are visible to the global update and to other
      strands making spatial queries. Any variable inside the scope of
      update{} will not be visible in this way. NOTE: "pos" is a special
      variable name that *must* be used to enable spatial queries of
      neighboring particles via sphere(). */
   output vec2 pos = pos0/|pos0|;
   real hh = hhInit;
   vec2 step = [0,0];  // step along force
   real closest = 0;   // distance to closest neighbor
   int ncount = 0;     // how many neighbors did we have

   update {
      // Compute energy and forces on us from neighbors
      real energyLast = 0;
      vec2 force = [0,0];
      ncount = 0;
      foreach (particle P in sphere(rad)) {
         energyLast += enr(P.pos - pos);
         force += frc(P.pos - pos);
         ncount += 1;
      }
      if (0 == ncount) {
         /* no neighbors; do nothing to do but set closest to big value */
         closest = rad;
         continue;
      }

      // Else we have interacting neighbors
      vec2 norm = normalize(pos); // surface normal for unit circle
      // project force onto tangent surface
      force = (identity[2] - norm⊗norm)•force;

      /* Limiting the step size (even before testing for sufficient decrease,
         below) helps keep particles near the surface they are supposed to be
         sampling; this precaution matters more with a non-trivial
         (data-driven) constraint surface */
      step = hh*force;           // compute step along force
      if (|step| > stepMax) {
         // decrease hh by factor by which step was too big
         hh *= stepMax/|step|;
         // and find smaller step (of length stepMax)
         step = hh*force;
      }

      // take step and re-apply constraint
      vec2 posLast = pos;
      pos = normalize(pos + step);
      // find energy at new location, and distance to closest neighbor
      real energy = 0;
      closest = rad;
      foreach (particle P in sphere(rad)) {
         energy += enr(P.pos - pos);
         closest = min(closest, |P.pos - pos|);
      }
      // save actual step taken
      step = pos - posLast;

      if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
         hh *= 0.5; // backtrack
         pos = posLast;
      } else {
         hh *= 1.02; // opportunistically increase step size
      }
   }
}

global {
   real meancl = mean { P.closest | P in particle.all};
   real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
   real covcl = sqrt(varicl) / meancl;
   print("(iter ", iter, ") mean(hh)=", mean { P.hh | P in particle.all},
         "; COV(closest) = ", covcl, "\n");
   if (covcl < eps) {
      print("Stabilizing with COV(closest) ", covcl, " < ", eps, " (iter ", iter, ")\n");
      stabilize;
   }
   iter += 1;
}

initially { particle(ipos[ii]) | ii in 0 .. length(ipos)-1 };

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