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

# SCM Repository

[diderot] Diff of /tests/examples/sphere/sphere.diderot
 [diderot] / tests / examples / sphere / sphere.diderot # Diff of /tests/examples/sphere/sphere.diderot

revision 4682, Mon Oct 3 16:25:37 2016 UTC revision 4683, Mon Oct 3 16:27:34 2016 UTC
# Line 8  Line 8
8  input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd");  input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd");
9  input real rad ("radius of particle's energy support") = 0.1;  input real rad ("radius of particle's energy support") = 0.1;
10  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;
11    input int pcp ("periodicity of population control (if non-zero)") = 5;
12    input real hhInit ("initial step size for potential energy gradient descent") = 1;
13
14  real hhInit = 1;        // initial step size  real newDist = 0.5*rad; // how far away to birth new particles
15  real stepMax = rad;     // limit on distance to travel per iter  real stepMax = rad;     // limit on distance to travel per iter
16  int iter = 0;           // which iteration we're on  int iter = 0;           // which iteration we're on
17
18  /* phi defines the potential energy around a particle, as a function of  /* energy functions */
distance from particle r > 0. By design, the potential function smoothly
falls to 0 at r == 1, so that the potential has compact support, which
simplifies computing the inter-particle interactions. phi'(r) is the
derivative of phi(r). The ' in phi' is merely suggestive; ' is not acting
as an operator. Below are some different possibilities for phi and phi';
only one should be uncommented.

In a conservative force field, force is the *negative* gradient of some
potential. phi(r) is potential at r > 0, due to a particle at 0. It is also
the potential at 0, due to a particle at r. Moving closer to the particle
at r by a small dr > 0, the change in potential is dphi = -phi'(r)*dr,
which is positive for a monotonically decreasing phi(r). The force, or
negative gradient of the potential, at 0 is then phi'(r), which justifies
how no negation is used in the frc() definition below. */
19
20  // phi(0) and phi'(0) are bounded  // phi(0) and phi'(0) are bounded
21  function real  phi(real r) = (1 - r)^4;  function real  phi(real r) = (1 - r)^4;
# Line 44  Line 32
32  function real phi'(real r) = (π/2)*(1 - (1/sin(r*π/2))^2);  function real phi'(real r) = (π/2)*(1 - (1/sin(r*π/2))^2);
33  */  */
34
35  /* The enr and frc functions use phi and phi' to define the potential  // Energy and force from particle at vec3 x
due to, and force from, a particle at 2-D offset x with a potential
36  function real enr(vec3 x) = phi(|x|/rad);  function real enr(vec3 x) = phi(|x|/rad);
37  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
38
39  /* The particle is initialized at position pos0, with stepsize hh0 */  // Returns a non-zero vector perpendicular to given non-zero vector v
40    function vec3 perp3(vec3 v) {
41       int c = 0;
42       if (|v| < |v|) {
43          c = 1;
44       }
45       // not v[c] because tensors can only be indexed by constants
46       if (|v if 1==c else v| < |v|) {
47          c = 2;
48       }
49       // now c is index of v component with largest absolute value
50       vec3 ret = ([v - v, -v, v] if (c == 0) else
51                   [-v, v - v, v] if (c == 1) else
52                   [-v, v, v - v]);
53       return ret;
54    }
55
56    // From a given vec3, find a random-ish value uniformly sampling [0,1)
57    function real posrnd(vec3 v) {
58       vec3 p = 10000*v/rad;
59       return fmod(|fmod(p,1)| + |fmod(p,1)| + |fmod(p,1)|, 1);
60    }
61
62    // Is this an iteration in which to do population control?
63    function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp);
64
65    /* The particle is initialized at position pos0, with initial stepsize hh0.
66       The initial set of particles is passed hhInit for initial stepsize, but
67       new particles created by population control benefit from the stepsize that
68       was adaptively learned by the parent */
69  strand particle (vec3 pos0, real hh0) {  strand particle (vec3 pos0, real hh0) {
/* "pos" is a special variable name that must be used to enable
spatial queries of neighboring particles via sphere() */
70     output vec3 pos = pos0/|pos0|;     output vec3 pos = pos0/|pos0|;
71     real hh = hh0;     real hh = hh0;
72     vec3 step = [0,0,0]; // step along force     vec3 step = [0,0,0]; // step along force
73     real closest = 0;   // distance to closest neighbor     real closest = 0;   // distance to closest neighbor
74       int ncount = 0;
75       /* This "done" variable is a simple way to signal to the global
76          update that something has just changed in a way that should
77          prevent convergence; in this program it is reset to 0 when new
78          particles are created, and otherwise incrementally increased */
79       real done = 0;
80     update {     update {
81        // compute energy and forces on us from neighbors        // compute energy and forces on us from neighbors
82        real energyLast = 0;        real energyLast = 0;
83        vec3 force = [0,0,0];        vec3 force = [0,0,0];
84        int ncount = 0;        ncount = 0;
85        foreach (particle P in sphere(rad)) {        foreach (particle P in sphere(rad)) {
86           energyLast += enr(P.pos - pos);           energyLast += enr(P.pos - pos);
87           force += frc(P.pos - pos);           force += frc(P.pos - pos);
88           ncount += 1;           ncount += 1;
89        }        }
90        if (0 == ncount) {        vec3 norm = normalize(pos); // surface normal for unit circle
91           /* no neighbors; do nothing to do but set closest to big value */        if (0 == ncount && pcIter()) {
92           closest = rad;           // no neighbors, so let's make one
93             vec3 npos = normalize(pos + newDist*normalize(perp3(norm)));
94             new particle(npos, hh);
95             // have to set closest to something; used in global update
96             closest = |npos - pos|;
97             done = 0;
98           continue;           continue;
99        }        }
100
101        /* Else we have interacting neighbors */        /* Else we have interacting neighbors; project force onto
102        vec3 norm = normalize(pos); // surface normal for unit circle           tangent surface, find step, limit step size */
// projection onto tangent surface
103        force = (identity - norm⊗norm)•force;        force = (identity - norm⊗norm)•force;

// Limit step size
104        step = hh*force;        step = hh*force;
105        if (|step| > stepMax) {        if (|step| > stepMax) {
106           hh *= stepMax/|step|;           hh *= stepMax/|step|;
107           step = hh*force;           step = hh*force;
108        }        }
109
110        // take step, re-apply constraint, find new energy        // Take step, re-apply constraint, find new energy
111        vec3 posLast = pos;        vec3 posLast = pos;
112        pos = normalize(pos + step);        pos = normalize(pos + step);
113        real energy = 0;        real energy = 0;
114        closest = rad;        closest = rad;
115          ncount = 0;
116        foreach (particle P in sphere(rad)) {        foreach (particle P in sphere(rad)) {
117           energy += enr(P.pos - pos);           energy += enr(P.pos - pos);
118           closest = min(closest, |P.pos - pos|);           closest = min(closest, |P.pos - pos|);
119             ncount += 1;
120        }        }
121        step = pos - posLast;        done = lerp(done, 1, 0.5);  // may be reset below
122
123        /* The Armijo-Goldstein sufficient decrease condition */        // Armijo-Goldstein sufficient decrease condition
124        if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {        if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
125           hh *= 0.5; // backtrack           hh *= 0.5;  // energy didn't decrease as expected, backtrack
126           pos = posLast;           pos = posLast;
127        } else {        } else {
128           hh *= 1.02; // opportunistically increase step size           hh *= 1.02; // opportunistically increase step size
129             /* try to have between 5 and 8 neighbors */
130             if (pcIter()) {
131                if (ncount < 5) {
132                   new particle(pos + newDist*normalize(force), hh);
133                   done = 0;
134                } else if (ncount > 8) {
135                   /* If this particles has ncount neighbors, then all of those
136                      neighbors probably have a similar number of neighbors. To
137                      get down to having about 8 neighbors, all of them should
138                      die with a chance of ncount-8 out of ncount. */
139                   if (posrnd(pos) < (ncount - 8.0)/ncount) {
140                      die;
141                   }
142                   done = 0;
143                }
144        }        }
145     }     }
146
147          // Record actual step taken, in case used in global update
148          step = pos - posLast;
149       }
150  }  }
151
152  global {  global {
//print("(iter ", iter, ") hello from global\n");
153     /* could test convergence based on movement */     /* could test convergence based on movement */
154     //real mvmt = max { |P.step|/rad | P in particle.all};     //real mvmt = max { |P.step|/rad | P in particle.all};
155     /* testing convergence based on distance to closest neighbor */     /* This tests convergence based on distance to closest neighbor */
156     real meancl = mean { P.closest | P in particle.all};     real meancl = mean { P.closest | P in particle.all};
157     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
158     real covcl = sqrt(varicl) / meancl;     real covcl = sqrt(varicl) / meancl; // coefficient of variation
159     print("   mean hh = ", mean { P.hh | P in particle.all}, "; covcl = ", covcl, "\n");     real meanncount = mean  { real(P.ncount) | P in particle.all};
160     if (covcl < eps) {     real mindone = min  { P.done | P in particle.all};
161        print("Stabilizing with COV(closest) ", covcl, " < ", eps, " (iter ", iter, ")\n");     print("(iter ", iter, ") COV(cl)=", covcl,
162             "; mean(ncount)=", meanncount, "; min(done)=", mindone, "\n");
163       if (covcl < eps          // seem to be geometrically uniform
164           && mindone > 0.9) {  // and no particle recently set done=0
165          print("Stabilizing ", numActive(), " points with COV(closest) ", covcl,
166                " < ", eps, " and mindone > ", mindone, " (iter ", iter, ")\n");
167        stabilize;        stabilize;
168     }     }
169     iter += 1;     iter += 1;

Legend:
 Removed from v.4682 changed lines Added in v.4683

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