SCM Repository
[diderot] / tests / examples / sphere / sphere.diderot |
View of /tests/examples/sphere/sphere.diderot
Parent Directory | Revision Log
Revision 4675 -
(download)
(annotate)
Sat Oct 1 12:36:25 2016 UTC (2 years, 11 months ago) by glk
File size: 4920 byte(s)
Sat Oct 1 12:36:25 2016 UTC (2 years, 11 months ago) by glk
File size: 4920 byte(s)
another particle example, starting with minor modification of circle
/* ========================================== Mutually-repelling particles on a unit sphere This is based on the unit circle example; see that program for more detailed and explanatory comments. ========================================== */ 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; 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 defines the potential energy around a particle, as a function of 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. */ // 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); */ /* The enr and frc functions use phi and phi' to define the potential due to, and force from, a particle at 2-D offset x with a potential field extending to radius rad. */ function real enr(vec3 x) = phi(|x|/rad); function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule /* The particle is initialized at position pos0, with stepsize 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() */ output vec3 pos = pos0/|pos0|; real hh = hh0; vec3 step = [0,0,0]; // step along force real closest = 0; // distance to closest neighbor update { // compute energy and forces on us from neighbors real energyLast = 0; vec3 force = [0,0,0]; int 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 */ vec3 norm = normalize(pos); // surface normal for unit circle // projection onto tangent surface force = (identity[3] - norm⊗norm)•force; // Limit step size 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; foreach (particle P in sphere(rad)) { energy += enr(P.pos - pos); closest = min(closest, |P.pos - pos|); } step = pos - posLast; /* The Armijo-Goldstein sufficient decrease condition */ if (energy - energyLast > 0.5*(pos - posLast)•(-force)) { hh *= 0.5; // backtrack pos = posLast; } else { hh *= 1.02; // opportunistically increase step size } } } global { //print("(iter ", iter, ") hello from global\n"); /* could test convergence based on movement */ //real mvmt = max { |P.step|/rad | P in particle.all}; /* testing convergence based on 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; print(" mean hh = ", mean { P.hh | P in particle.all}, "; covcl = ", covcl, "\n"); if (covcl < eps) { print("Stabilizing with COV(closest) ", covcl, " < ", eps, " (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 |