SCM Repository
[diderot] / tests / examples / sphere / sphere.diderot |
View of /tests/examples/sphere/sphere.diderot
Parent Directory | Revision Log
Revision 4686 -
(download)
(annotate)
Mon Oct 3 16:46:39 2016 UTC (2 years, 9 months ago) by glk
File size: 6880 byte(s)
Mon Oct 3 16:46:39 2016 UTC (2 years, 9 months ago) by glk
File size: 6880 byte(s)
comment tweaks
/* ========================================== Mutually-repelling particles on a unit sphere This is based on the unit circle example; see that program for more detailed and explanatory comments. The new things added with this example are documented in comments below. The most significant new capability is "population control", whereby particles create new particles if there are too few (using "new"), or die if there are too may (using "die"). ========================================== */ input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd"); input real rad ("radius of particle's energy support") = 0.1; input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05; input int pcp ("periodicity of population control (if non-zero)") = 5; input real hhInit ("initial step size for potential energy gradient descent") = 1; real newDist = 0.5*rad; // how far away to birth new particles real stepMax = rad; // limit on distance to travel per iter int iter = 0; // which iteration we're on /* energy functions */ // phi(0) and phi'(0) are bounded function real phi(real r) = (1 - r)^4; function real phi'(real r) = -4*(1 - r)^3; /* // electrostatic potential 1/r, scaled to be C^2 continuous with 0 at r==1 function real phi(real r) = (1/r)*(1 - r)^3; function real phi'(real r) = 3 - 1/(r^2) - 2*r; */ /* // Cotangent-based potential from Meyer et al. SMI'05 function real phi(real r) = 1/tan(r*π/2) + r*π/2 - π/2; function real phi'(real r) = (π/2)*(1 - (1/sin(r*π/2))^2); */ // Energy and force from particle at vec3 x function real enr(vec3 x) = phi(|x|/rad); function vec3 frc(vec3 x) = phi'(|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; } // not 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; } // From a given vec3, find a random-ish value uniformly sampling [0,1) function real posrnd(vec3 v) { vec3 p = 10000*v/rad; return fmod(|fmod(p[0],1)| + |fmod(p[1],1)| + |fmod(p[2],1)|, 1); } // Is this an iteration in which to do population control? function bool pcIter() = (pcp > 0 && iter > 0 && 0 == iter % pcp); /* The particle is initialized at position pos0, with initial stepsize hh0. The first set of particles gets hhInit for the initial stepsize, but new particles created by population control benefit from getting the stepsize that was adaptively learned by the parent. */ strand particle (vec3 pos0, real hh0) { output vec3 pos = pos0/|pos0|; real hh = hh0; vec3 step = [0,0,0]; // step along force real closest = 0; // distance to closest neighbor int ncount = 0; // how many neighbors did we have /* This "done" variable is a simple way to signal to the global update that something is happening or just changed in a way that should prevent convergence. In this program it is reset to 0 when new particles are created and when there are too many neighbors; otherwise it is incrementally increased towards 1. */ real done = 0; update { // compute energy and forces on us from neighbors real energyLast = 0; vec3 force = [0,0,0]; ncount = 0; foreach (particle P in sphere(rad)) { energyLast += enr(P.pos - pos); force += frc(P.pos - pos); ncount += 1; } vec3 norm = normalize(pos); // surface normal for unit circle if (0 == ncount && pcIter()) { // no neighbors, so let's make one vec3 npos = normalize(pos + newDist*normalize(perp3(norm))); new particle(npos, hh); // set closest to something in case used in global update closest = |npos - pos|; done = 0; continue; } /* Else we have interacting neighbors; project force onto tangent surface, find step, limit step size */ force = (identity[3] - norm⊗norm)•force; step = hh*force; if (|step| > stepMax) { hh *= stepMax/|step|; step = hh*force; } // Take step, re-apply constraint, find new energy vec3 posLast = pos; pos = normalize(pos + step); real energy = 0; closest = rad; ncount = 0; foreach (particle P in sphere(rad)) { energy += enr(P.pos - pos); closest = min(closest, |P.pos - pos|); ncount += 1; } done = lerp(done, 1, 0.5); // may be reset below // Armijo-Goldstein sufficient decrease condition if (energy - energyLast > 0.5*(pos - posLast)•(-force)) { hh *= 0.5; // energy didn't decrease as expected, so backtrack pos = posLast; } else { hh *= 1.02; // opportunistically increase step size /* try to have between 5 and 8 neighbors */ if (pcIter()) { if (ncount < 5) { new particle(pos + newDist*normalize(force), hh); done = 0; } else if (ncount > 8) { /* If this particle has ncount neighbors, then all of those neighbors probably have a similar number of neighbors. To get down to having about 8 neighbors, all of them should die with a chance of ncount-8 out of ncount. */ if (posrnd(pos) < (ncount - 8.0)/ncount) { die; } done = 0; // else not done if too many neighbors } } } // Record actual step taken, in case used in global update step = pos - posLast; } } global { /* Compute coefficient-of-variation of distance to closest neighbor 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; real meanncount = mean { real(P.ncount) | P in particle.all}; real mindone = min { P.done | P in particle.all}; print("(iter ", iter, ") COV(cl)=", covcl, "; mean(ncount)=", meanncount, "; min(done)=", mindone, "\n"); if (covcl < eps // seem to be geometrically uniform && mindone > 0.9) { // and no particle recently set done=0 print("Stabilizing ", numActive(), " points with COV(closest) ", covcl, " < ", eps, " and mindone > ", mindone, " (iter ", iter, ")\n"); stabilize; } iter += 1; } initially { particle(ipos[ii], hhInit) | ii in 0 .. length(ipos)-1 };
root@smlnj-gforge.cs.uchicago.edu | ViewVC Help |
Powered by ViewVC 1.0.0 |