Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] View of /tests/examples/sphere/sphere.diderot
 [diderot] / tests / examples / sphere / sphere.diderot

# View of /tests/examples/sphere/sphere.diderot

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
========================================== */

input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd");
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
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 */
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;
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 };
```