Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] View of /tests/lamont-tests/unit-sphere-bug1.diderot
 [diderot] / tests / lamont-tests / unit-sphere-bug1.diderot

# View of /tests/lamont-tests/unit-sphere-bug1.diderot

Tue Oct 11 17:28:35 2016 UTC (2 years, 10 months ago) by glk
File size: 4994 byte(s)
`these used to crash the compiler; running the code is not what revealed the bug`
```// unit-sphere
//
// Demo of distributing particles on the unit sphere
//
input vec3{} initPos ("initial positions for all particles") = load("vec3.nrrd");
input int iterMax ("maximum number of iterations to run") = 100;
input int pcp ("periodicity of population control (if non-zero)") = 5;
input real mvmtEps ("convergence threshold on particle system movement") = 0.001;
input bool noDie ("all strands stabilize instead of die at iterMax") = false;

real hhInit = 10;       // initial step size, can be big
real newDist = 0.9*rad; // how far away to birth new particles
real stepMax = 2*rad;   // limit on distance to travel per iter

int iter = 0;           // which iteration we're on

// inter-particle energy, and its derivative
//function real  φ(real r) = (1 - r)^4;
//function real φ'(real r) = -4*(1 - r)^3;
function real  φ(real r)  =  (1/r)*(1-r)^3;
function real φ'(real r) = 3 - 1/(r^2) - 2*r;
function real enr(vec3 x) = φ(|x|/rad);
function vec3 frc(vec3 x) = φ'(|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;
}
// can't say 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;
}

function real prnd(vec3 v) {
real ret = fmod(v[0]*1000, 1);
print(v[0], " --> ", ret, "\n");
return ret;
}

strand Particle (vec3 pos0, real hh0) {
vec3 pos = pos0/|pos0|;
real hh = hh0;
output vec3 outPos = pos;
/* putting this here allows it to be accessed by the global update,
where it is used in a convergence test */
vec3 delta = [0,0,0];
real dummy = 0;
update {
// compute energy and forces on us
int ncount = 0;
real energyLast = 0;
vec3 force = [0,0,0];
foreach (Particle p_j in sphere(rad)) {
ncount += 1;
vec3 r_ij = p_j.pos - pos;
energyLast += enr(r_ij);
force += frc(r_ij);
}
vec3 norm = normalize(pos); // surface normal (for unit circle)
if (ncount == 0 && (pcp > 0 && 0 == iter % pcp)) {
vec3 npos = pos + newDist*rad*normalize(perp3(norm));
new Particle(npos, hh);
continue;
}
// else update position based on force from neighbors.
// first, project force onto tangent surface
force = (identity[3] - norm⊗norm)•force;
delta = hh*force;
if (|delta| > stepMax) {
// decrease hh by factor by which step was too big
hh *= stepMax/|delta|;
// and find smaller step
delta = hh*force; // == stepMax*normalize(force);
}
vec3 posLast = pos;
// take step and re-apply constraint
pos = normalize(pos + delta);
real energy = 0;
ncount = 0;
real minz = 2; // like infinity
foreach (Particle p_j in sphere(rad)) {
energy += enr(p_j.pos - pos);
ncount += 1;
minz = min(minz, p_j.pos[2]);
dummy = prnd(p_j.pos);
}
if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
hh *= 0.5; // backtrack
pos = posLast;
} else {
hh *= 1.1;
/* putting new particle creation here avoids possibility of
starting multiple particles at same position, while this
particle is not moving due to backtracking (above) */
if (pcp > 0 && 0 == iter % pcp) {
if (ncount < 5) {
vec3 npos = pos + newDist*rad*normalize(force);
new Particle(npos, hh);
} else if (ncount > 8) {
if (pos[2] < minz) { die; }
}
}
}
outPos = pos;
}
}

global {
print("(iter ", iter, ") hello from global\n");
real dmm = sum { P.dummy| P in Particle.all};
print("    ", dmm, "\n");
real mvmt = max { |P.delta|/rad | P in Particle.all};
if (mvmt < mvmtEps          // there is very little movement
&& (pcp == 0            // and, either there is no population control
|| iter > 2*pcp)) {  // or new particles should be moving by now
print("Stabilizing ", numActive(), " points with movement ", mvmt,
" < ", mvmtEps, " (iter ", iter, ")\n");
stabilize;
}
iter += 1;
if (iter >= iterMax && iterMax > 0)  {
print("Stabilizing (due to noDie) " if noDie else "STOPPING ",
numActive(), " points at iter ", iter, " >= ", iterMax,
" (movement ", mvmt, ")\n");
if (noDie) {
stabilize;
} else {
die;
}
}
}

initially {Particle(initPos[ii], hhInit)
| ii in 0 .. length(initPos)-1 };

```