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

SCM Repository

[diderot] Diff of /tests/examples/sphere/sphere.diderot
ViewVC logotype

Diff of /tests/examples/sphere/sphere.diderot

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4691, Mon Oct 3 21:36:42 2016 UTC revision 4692, Tue Oct 4 17:04:27 2016 UTC
# Line 1  Line 1 
1  /* ==========================================  /* ==========================================
2  Mutually-repelling particles on a unit sphere  Mutually-repelling particles populating a unit sphere
3    
4  This is based on the unit circle example; see that program for more detailed  This is heavily based on the [`circle.diderot`](../circle) example; see that
5  and explanatory comments. The new things added with this example are  program for more detailed and explanatory comments. The new things added with
6  documented in comments below. The most significant new capability is  this example are documented in comments below. A significant capability
7  "population control", whereby particles create new particles if there are too  demonstrated here is "population control", whereby particles create new
8  few (using "new"), or die if there are too may (using "die").  particles (using `new`) if there seem to be too few, or die if there are too
9    may (using `die`).  Along with this there is a new variable `undone` that acts
10    as a signal (to the global update method monitoring the computation) that
11    things are unfinished because the system population is changing.
12    
13  ... TODO: example of running program and looking at output  ... TODO: example of running program and looking at output
14    
# Line 14  Line 17 
17  input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd");  input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd");
18  input real rad ("radius of particle's potential energy support") = 0.1;  input real rad ("radius of particle's potential energy support") = 0.1;
19  input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;  input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;
20  input int pcp ("periodicity of population control (if non-zero)") = 5;  input int pcp ("periodicity of population control (if greater than zero)") = 2;
21  input real hhInit ("initial step size for potential energy gradient descent") = 1;  input real hhInit ("initial step size for potential energy gradient descent") = 1;
22    
23  real newDist = 0.5*rad; // how far away to birth new particles  real newDist = 0.5*rad; // how far away to birth new particles
24  real stepMax = rad;     // limit on distance to travel per iter  real stepMax = rad;     // limit on distance to travel per iter
25  int iter = 0;           // which iteration we're on  int iter = 0;           // which iteration we're on
26    
27  // Univariate energy functions; see ../circle for alternatives  // Univariate energy functions; see ../circle/circle.diderot for alternatives
28  function real  phi(real r) = (1 - r)^4;  function real  phi(real r) = (1 - r)^4;
29  function real phi'(real r) = -4*(1 - r)^3;  function real phi'(real r) = -4*(1 - r)^3;
30  // Energy and force from particle at vec3 x  // Energy and force from particle (with radius rad) at vec3 x
31  function real enr(vec3 x) = phi(|x|/rad);  function real enr(vec3 x) = phi(|x|/rad);
32  function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule  function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule
33    
# Line 62  Line 65 
65     output vec3 pos = pos0/|pos0|;     output vec3 pos = pos0/|pos0|;
66     real hh = hh0;     real hh = hh0;
67     vec3 step = [0,0,0]; // step along force     vec3 step = [0,0,0]; // step along force
68     real closest = 0;    // distance to closest neighbor     real closest = rad;  // distance to closest neighbor
69     int ncount = 0;      // how many neighbors did we have     int ncount = 0;      // how many neighbors did we have
70     /* This "done" variable is a simple way to signal to the global update that     /* This "undone" variable signals to global update that something is
71        something is happening or just changed in a way that should prevent        happening or just changed that should delay convergence. In this program
72        convergence. In this program it is reset to 0 when new particles are        it is reset to 1 when new particles are created and when there are too
73        created and when there are too many neighbors; otherwise it is        many neighbors; otherwise it is slowly decreased towards 0. */
74        incrementally increased towards 1. */     real undone = 1;
    real done = 0;  
75     update {     update {
76        // compute energy and forces on us from neighbors        // compute energy and forces on us from neighbors
77        real energyLast = 0;        real energyLast = 0;
78        vec3 force = [0,0,0];        vec3 force = [0,0,0];
79        ncount = 0;        ncount = 0;
80        foreach (particle P in sphere(rad)) {        foreach (particle P in sphere(rad)) {
81           energyLast += enr(P.pos - pos);           vec3 x = P.pos - pos;
82           force += frc(P.pos - pos);           if (|x| == 0) {
83                /* we're exactly overlapping with another particle; would be
84                   nice to have exactly one strand persist and kill the others;
85                   but simpler to have all overlap-ees die here and let
86                   population control fill in the hole as needed later */
87                die;
88             }
89             energyLast += enr(x);
90             force += frc(x);
91           ncount += 1;           ncount += 1;
92        }        }
93        vec3 norm = normalize(pos); // surface normal for unit circle        vec3 norm = normalize(pos); // surface normal for unit circle
94        if (0 == ncount && pcIter()) {        if (0 == ncount) {
95           // no neighbors, so let's make one           if (pcIter()) {  // no neighbors, so let's make one
96           vec3 npos = normalize(pos + newDist*normalize(perp3(norm)));              vec3 npos = pos + newDist*normalize(perp3(norm));
97           new particle(npos, hh);           new particle(npos, hh);
98                undone = 1;
99             }
100           // set closest to something in case used in global update           // set closest to something in case used in global update
101           closest = |npos - pos|;           closest = newDist;
102           done = 0;           undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
103           continue;           continue;
104        }        }
105    
# Line 111  Line 123 
123           closest = min(closest, |P.pos - pos|);           closest = min(closest, |P.pos - pos|);
124           ncount += 1;           ncount += 1;
125        }        }
       done = lerp(done, 1, 0.5);  // may be reset below  
126    
127        // Armijo-Goldstein sufficient decrease condition        // Armijo-Goldstein sufficient decrease condition
128        if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {        if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
129           hh *= 0.5;  // energy didn't decrease as expected, so backtrack           hh *= 0.5;  // energy didn't decrease as expected: backtrack
130           pos = posLast;           pos = posLast;  // try again next iteration
131             // no progress, so decrease of undone
132        } else {        } else {
133           hh *= 1.02; // opportunistically increase step size           hh *= 1.02; // opportunistically increase step size
134           /* try to have between 5 and 8 neighbors */           // indicate progress; may be over-written below
135             undone *= pow(0.5, 1.0/(2*(pcp if pcp > 0 else 1)));
136             // try to have between 5 and 8 neighbors
137           if (pcIter()) {           if (pcIter()) {
138              if (ncount < 5) {              if (ncount < 5) {
139                 new particle(pos + newDist*normalize(force), hh);                 new particle(pos + newDist*normalize(force), hh);
140                 done = 0;                 undone = 1;
141              } else if (ncount > 8) {              } else if (ncount > 8) {
142                 /* If this particle has ncount neighbors, then all of those                 /* If this particle has ncount neighbors, then all of those
143                    neighbors probably have a similar number of neighbors. To                    neighbors probably have a similar number of neighbors. To
# Line 132  Line 146 
146                 if (posrnd(pos) < (ncount - 8.0)/ncount) {                 if (posrnd(pos) < (ncount - 8.0)/ncount) {
147                    die;                    die;
148                 }                 }
149                 done = 0; // else not done if too many neighbors                 // else not done if too many neighbors, w/ population control
150                   undone = 1;
151              }              }
152           }           }
153        }        }
# Line 148  Line 163 
163     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
164     real covcl = sqrt(varicl) / meancl;     real covcl = sqrt(varicl) / meancl;
165     real meanncount = mean  { real(P.ncount) | P in particle.all};     real meanncount = mean  { real(P.ncount) | P in particle.all};
166     real mindone = min  { P.done | P in particle.all};     real maxundone = max  { P.undone | P in particle.all};
167     print("(iter ", iter, ") COV(cl)=", covcl,     print("(iter ", iter, ") COV(cl)=", covcl,
168           "; mean(ncount)=", meanncount, "; min(done)=", mindone, "\n");           "; mean(ncount)=", meanncount, "; max(undone)=", maxundone, "\n");
169    
170     if (covcl < eps          // seem to be geometrically uniform     if (covcl < eps          // seem to be geometrically uniform
171         && mindone > 0.9) {  // and no particle recently set done=0         && maxundone < 0.5) { // and no particle recently set undone=1
172        print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,        print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,
173              " < ", eps, " and mindone > ", mindone, " (iter ", iter, ")\n");              " < ", eps, " and maxundone > ", maxundone, " (iter ", iter, ")\n");
174        stabilize;        stabilize;
175     }     }
176     iter += 1;     iter += 1;

Legend:
Removed from v.4691  
changed lines
  Added in v.4692

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