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

SCM Repository

[diderot] Annotation of /tests/vis15-bugs/circle.diderot
ViewVC logotype

Annotation of /tests/vis15-bugs/circle.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4917 - (view) (download)

1 : glk 4729 /*
2 :     running this under valgrind produces:
3 :     Syscall param write(buf) points to uninitialised byte(s)" that seems to be about the memory in the nrrd not being set
4 :     due to saving out the first (post-initialization, before first iteration) snapshot
5 :    
6 :     */
7 :    
8 : glk 4917 input vec2{} ipos ("initial positions for all particles") = load("circle-vec2.nrrd");
9 : glk 4729 input real rad ("radius of particle's potential energy support") = 0.25;
10 :     input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;
11 :    
12 :     real hhInit = 1; // initial step size
13 :     real stepMax = rad; // limit on distance to travel per iter
14 :     int iter = 0; // which iteration we're on
15 :    
16 :     // phi(0) and phi'(0) are bounded
17 :     function real phi(real r) = (1 - r)^4;
18 :     function real phi'(real r) = -4*(1 - r)^3;
19 :    
20 :     function real enr(vec2 x) = phi(|x|/rad);
21 :     function vec2 frc(vec2 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule
22 :    
23 :     /* The particle is initialized at position pos0 */
24 :     strand particle (vec2 pos0) {
25 :     /* These strand variables are visible to the global update and to other
26 :     strands making spatial queries. Any variable inside the scope of
27 :     update{} will not be visible in this way. NOTE: "pos" is a special
28 :     variable name that *must* be used to enable spatial queries of
29 :     neighboring particles via sphere(). */
30 :     output vec2 pos = pos0/|pos0|;
31 :     real hh = hhInit;
32 :     vec2 step = [0,0]; // step along force
33 :     real closest = 0; // distance to closest neighbor
34 :     int ncount = 0; // how many neighbors did we have
35 :    
36 :     update {
37 :     // Compute energy and forces on us from neighbors
38 :     real energyLast = 0;
39 :     vec2 force = [0,0];
40 :     ncount = 0;
41 :     foreach (particle P in sphere(rad)) {
42 :     energyLast += enr(P.pos - pos);
43 :     force += frc(P.pos - pos);
44 :     ncount += 1;
45 :     }
46 :     if (0 == ncount) {
47 :     /* no neighbors; do nothing to do but set closest to big value */
48 :     closest = rad;
49 :     continue;
50 :     }
51 :    
52 :     // Else we have interacting neighbors
53 :     vec2 norm = normalize(pos); // surface normal for unit circle
54 :     // project force onto tangent surface
55 :     force = (identity[2] - norm⊗norm)•force;
56 :    
57 :     /* Limiting the step size (even before testing for sufficient decrease,
58 :     below) helps keep particles near the surface they are supposed to be
59 :     sampling; this precaution matters more with a non-trivial
60 :     (data-driven) constraint surface */
61 :     step = hh*force; // compute step along force
62 :     if (|step| > stepMax) {
63 :     // decrease hh by factor by which step was too big
64 :     hh *= stepMax/|step|;
65 :     // and find smaller step (of length stepMax)
66 :     step = hh*force;
67 :     }
68 :    
69 :     // take step and re-apply constraint
70 :     vec2 posLast = pos;
71 :     pos = normalize(pos + step);
72 :     // find energy at new location, and distance to closest neighbor
73 :     real energy = 0;
74 :     closest = rad;
75 :     foreach (particle P in sphere(rad)) {
76 :     energy += enr(P.pos - pos);
77 :     closest = min(closest, |P.pos - pos|);
78 :     }
79 :     // save actual step taken
80 :     step = pos - posLast;
81 :    
82 :     if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
83 :     hh *= 0.5; // backtrack
84 :     pos = posLast;
85 :     } else {
86 :     hh *= 1.02; // opportunistically increase step size
87 :     }
88 :     }
89 :     }
90 :    
91 :     global {
92 :     real meancl = mean { P.closest | P in particle.all};
93 :     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
94 :     real covcl = sqrt(varicl) / meancl;
95 :     print("(iter ", iter, ") mean(hh)=", mean { P.hh | P in particle.all},
96 :     "; COV(closest) = ", covcl, "\n");
97 :     if (covcl < eps) {
98 :     print("Stabilizing with COV(closest) ", covcl, " < ", eps, " (iter ", iter, ")\n");
99 :     stabilize;
100 :     }
101 :     iter += 1;
102 :     }
103 :    
104 :     initially { particle(ipos[ii]) | ii in 0 .. length(ipos)-1 };

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