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

SCM Repository

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

Annotation of /tests/lamont-tests/unit-circle.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4640 - (view) (download)

1 : glk 4640 #version 2
2 :    
3 :     // unit-circle
4 :     //
5 :     // Demo of distributing particles on the unit circle
6 :     //
7 :     //real[] initPosns = load_sequence("-");
8 :     real[] initPosns = load_sequence("100.nrrd");
9 :     int numOfParticles = length(initPosns)/2;
10 :     input int iterMax = 300; // maximum number of steps to run
11 :     input real rr = 0.5; // actual particle radius
12 :     real stdev = ∞;
13 :     real threshold = 0.08; // covnergence threshold
14 :     real RR = rr+0.0; // neighbor query radius (MUST be >= rr; should get same results for any RR >= rr)
15 :     real hhInit = 10.0; // initial integration step size; can err too big, will be trimmed down during iterations
16 :     int iter = 1; // which iteration we're on
17 :     real stepMax = 2*π/numOfParticles; // limit on distance to travel per iter
18 :    
19 :     strand Particle (int ii, real posx, real posy) {
20 :     vec2 pos = normalize([posx,posy]);
21 :     real hh = hhInit;
22 :     vec2 posOld1 = pos; // remember last TWO positions
23 :     vec2 posOld2 = pos;
24 :     output vec2 outPos = pos;
25 :     int id = ii;
26 :     real energy = ii; // HEY: can do convergence test based on variance of energies
27 :     vec2 force = [0,0]; // or can test convergence based on sum of |force|
28 :     update {
29 :    
30 :     if( iter > iterMax) {
31 :     stabilize;
32 :     }
33 :    
34 :     // compute energy and forces on us
35 :     energy = 0;
36 :     force = [0,0];
37 :     foreach (Particle p_j in sphere(RR)) {
38 :     vec2 r_ij = (pos - p_j.pos)/rr;
39 :     vec2 d_ij = normalize(r_ij);
40 :     if (|r_ij| < 1) {
41 :     energy += (1 - |r_ij|)^4;
42 :     force += - (-4*(1 - |r_ij|)^3) * d_ij;
43 :     }
44 :     }
45 :     force /= rr; // smaller particles make larger forces
46 :    
47 :     // update position based on force
48 :     posOld2 = posOld1; // shuffle saved positions down
49 :     posOld1 = pos;
50 :     if (energy > 0.0) { // we have neighbors
51 :     tensor[2,2] pten = identity[2] - normalize(pos)⊗normalize(pos);
52 :     force = pten•force; // project force onto tangent surface
53 :     vec2 step = hh*force;
54 :     if (|step| > stepMax) {
55 :     // decrease hh by factor by which step was too big
56 :     hh *= stepMax/|step|;
57 :     // and find smaller step
58 :     step = hh*force;
59 :     }
60 :     // take step and re-find implicit surface
61 :     pos = normalize(pos + step);
62 :     real travel = |pos - posOld1| + |posOld1 - posOld2|;
63 :     if (travel > 0) {
64 :     // if we've moved in the past two steps, but we've moved back
65 :     // to where we were two steps ago, we're oscillating ==>
66 :     // okay = 0. Two steps in the same direction ==> okay = 1.
67 :     real okay = |pos - posOld2|/travel;
68 :     // slow down if oscillating, speed up a little if not
69 :     hh *= lerp(0.8, 1.001, 0, okay, 1);
70 :     }
71 :     }
72 :     if (id == 0 && iter == 140) {
73 :     new Particle(2,pos[0] + 0.001, pos[1] + 0.01);
74 :     }
75 :     outPos = pos;
76 :     }
77 :     }
78 :    
79 :     // can add statements in here to do some global computation
80 :     update {
81 :     // real var = variance{P.energy | P in Particle.all};
82 :     real avg = mean{P.energy | P in Particle.all};
83 :     print("iter ", iter, ": avg ", avg, "\n");
84 :     iter+=1;
85 :     real var = sum{(avg-P.energy)^2 | P in Particle.all};
86 :     stdev = sqrt(var);
87 :     }
88 :    
89 :     create_collection {Particle(ii, initPosns[ii*2], initPosns[ii*2+1])
90 :     | ii in 0 .. numOfParticles-1 }
91 :    

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