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

SCM Repository

[diderot] Annotation of /tests/lamont-tests/unit-sphere-bug1.diderot
ViewVC logotype

Annotation of /tests/lamont-tests/unit-sphere-bug1.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4712 - (view) (download)

1 : glk 4712 // unit-sphere
2 :     //
3 :     // Demo of distributing particles on the unit sphere
4 :     //
5 :     input vec3{} initPos ("initial positions for all particles") = load("vec3.nrrd");
6 :     input int iterMax ("maximum number of iterations to run") = 100;
7 :     input real rad ("radius of particle's energy support") = 0.1;
8 :     input int pcp ("periodicity of population control (if non-zero)") = 5;
9 :     input real mvmtEps ("convergence threshold on particle system movement") = 0.001;
10 :     input bool noDie ("all strands stabilize instead of die at iterMax") = false;
11 :    
12 :     real hhInit = 10; // initial step size, can be big
13 :     real newDist = 0.9*rad; // how far away to birth new particles
14 :     real stepMax = 2*rad; // limit on distance to travel per iter
15 :    
16 :     int iter = 0; // which iteration we're on
17 :    
18 :     // inter-particle energy, and its derivative
19 :     //function real φ(real r) = (1 - r)^4;
20 :     //function real φ'(real r) = -4*(1 - r)^3;
21 :     function real φ(real r) = (1/r)*(1-r)^3;
22 :     function real φ'(real r) = 3 - 1/(r^2) - 2*r;
23 :     function real enr(vec3 x) = φ(|x|/rad);
24 :     function vec3 frc(vec3 x) = φ'(|x|/rad) * (1/rad) * x/|x|; // chain rule
25 :    
26 :     // returns a non-zero vector perpendicular to given non-zero vector v
27 :     function vec3 perp3(vec3 v) {
28 :     int c = 0;
29 :     if (|v[0]| < |v[1]|) {
30 :     c = 1;
31 :     }
32 :     // can't say v[c] because tensors can only be indexed by constants
33 :     if (|v[1] if 1==c else v[0]| < |v[2]|) {
34 :     c = 2;
35 :     }
36 :     // now c is index of v component with largest absolute value
37 :     vec3 ret = ([v[1] - v[2], -v[0], v[0]] if (c == 0) else
38 :     [-v[1], v[0] - v[2], v[1]] if (c == 1) else
39 :     [-v[2], v[2], v[0] - v[1]]);
40 :     return ret;
41 :     }
42 :    
43 :     function real prnd(vec3 v) {
44 :     real ret = fmod(v[0]*1000, 1);
45 :     print(v[0], " --> ", ret, "\n");
46 :     return ret;
47 :     }
48 :    
49 :     strand Particle (vec3 pos0, real hh0) {
50 :     vec3 pos = pos0/|pos0|;
51 :     real hh = hh0;
52 :     output vec3 outPos = pos;
53 :     /* putting this here allows it to be accessed by the global update,
54 :     where it is used in a convergence test */
55 :     vec3 delta = [0,0,0];
56 :     real dummy = 0;
57 :     update {
58 :     // compute energy and forces on us
59 :     int ncount = 0;
60 :     real energyLast = 0;
61 :     vec3 force = [0,0,0];
62 :     foreach (Particle p_j in sphere(rad)) {
63 :     ncount += 1;
64 :     vec3 r_ij = p_j.pos - pos;
65 :     energyLast += enr(r_ij);
66 :     force += frc(r_ij);
67 :     }
68 :     vec3 norm = normalize(pos); // surface normal (for unit circle)
69 :     if (ncount == 0 && (pcp > 0 && 0 == iter % pcp)) {
70 :     vec3 npos = pos + newDist*rad*normalize(perp3(norm));
71 :     new Particle(npos, hh);
72 :     continue;
73 :     }
74 :     // else update position based on force from neighbors.
75 :     // first, project force onto tangent surface
76 :     force = (identity[3] - norm⊗norm)•force;
77 :     delta = hh*force;
78 :     if (|delta| > stepMax) {
79 :     // decrease hh by factor by which step was too big
80 :     hh *= stepMax/|delta|;
81 :     // and find smaller step
82 :     delta = hh*force; // == stepMax*normalize(force);
83 :     }
84 :     vec3 posLast = pos;
85 :     // take step and re-apply constraint
86 :     pos = normalize(pos + delta);
87 :     real energy = 0;
88 :     ncount = 0;
89 :     real minz = 2; // like infinity
90 :     foreach (Particle p_j in sphere(rad)) {
91 :     energy += enr(p_j.pos - pos);
92 :     ncount += 1;
93 :     minz = min(minz, p_j.pos[2]);
94 :     dummy = prnd(p_j.pos);
95 :     }
96 :     if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
97 :     hh *= 0.5; // backtrack
98 :     pos = posLast;
99 :     } else {
100 :     hh *= 1.1;
101 :     /* putting new particle creation here avoids possibility of
102 :     starting multiple particles at same position, while this
103 :     particle is not moving due to backtracking (above) */
104 :     if (pcp > 0 && 0 == iter % pcp) {
105 :     if (ncount < 5) {
106 :     vec3 npos = pos + newDist*rad*normalize(force);
107 :     new Particle(npos, hh);
108 :     } else if (ncount > 8) {
109 :     if (pos[2] < minz) { die; }
110 :     }
111 :     }
112 :     }
113 :     outPos = pos;
114 :     }
115 :     }
116 :    
117 :     global {
118 :     print("(iter ", iter, ") hello from global\n");
119 :     real dmm = sum { P.dummy| P in Particle.all};
120 :     print(" ", dmm, "\n");
121 :     real mvmt = max { |P.delta|/rad | P in Particle.all};
122 :     if (mvmt < mvmtEps // there is very little movement
123 :     && (pcp == 0 // and, either there is no population control
124 :     || iter > 2*pcp)) { // or new particles should be moving by now
125 :     print("Stabilizing ", numActive(), " points with movement ", mvmt,
126 :     " < ", mvmtEps, " (iter ", iter, ")\n");
127 :     stabilize;
128 :     }
129 :     iter += 1;
130 :     if (iter >= iterMax && iterMax > 0) {
131 :     print("Stabilizing (due to noDie) " if noDie else "STOPPING ",
132 :     numActive(), " points at iter ", iter, " >= ", iterMax,
133 :     " (movement ", mvmt, ")\n");
134 :     if (noDie) {
135 :     stabilize;
136 :     } else {
137 :     die;
138 :     }
139 :     }
140 :     }
141 :    
142 :     initially {Particle(initPos[ii], hhInit)
143 :     | ii in 0 .. length(initPos)-1 };
144 :    

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