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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4640 - (view) (download)

1 : glk 4640 // unit-circle
2 :     //
3 :     // Demo of distributing particles on the unit circle
4 :     //
5 :     input vec2{} initPos = load("vec2.nrrd");
6 :     // BUG: input vec2{} initPos = {[0.15606558, -0.56803983]};
7 :     int numOfParticles = length(initPos);
8 :     input int iterMax = 1; // maximum number of steps to run
9 :     input real rad = 0.2; // particle radius
10 :     real hhInit = 1; // initial step size, can be big
11 :     int iter = 1; // which iteration we're on
12 :     real stepMax = 2*rad; // limit on distance to travel per iter
13 :    
14 :     // inter-particle energy, and its derivative
15 :     function real φ(real r) = (1 - r)^4;
16 :     function real φ'(real r) = -4*(1 - r)^3;
17 :     //function real φ(real r) = (1/r)*(1-r)^3;
18 :     //function real φ'(real r) = 3 - 1/(r^2) - 2*r;
19 :     function real enr(vec2 x) = φ(|x|/rad);
20 :     function vec2 frc(vec2 x) = φ'(|x|/rad) * (1/rad) * x/|x|; // chain rule
21 :    
22 :     strand Particle (int parent_, int ii, vec2 pos0, real hh0) {
23 :     int parent = parent_;
24 :     vec2 pos = pos0/|pos0|;
25 :     real hh = hh0;
26 :     output vec2 outPos = pos;
27 :     vec2 dpos = [0,0];
28 :     int id = ii;
29 :     int born = 0;
30 :     real energyLast = ii; // HEY: can do convergence test based on variance of energies
31 :     vec2 force = [0,0]; // or can test convergence based on sum of |force|
32 :     bool first = true;
33 :     update {
34 :     if (first) {
35 :     first = false;
36 :     }
37 :     // compute energy and forces on us
38 :     energyLast = 0;
39 :     force = [0,0];
40 :     int ncount = 0;
41 :     foreach (Particle p_j in sphere(rad)) {
42 :     ncount += 1;
43 :     vec2 r_ij = p_j.pos - pos;
44 :     energyLast += enr(r_ij);
45 :     force += frc(r_ij);
46 :     }
47 :     vec2 norm = normalize(pos);
48 :     // update position based on force
49 :     force = (identity[2] - norm⊗norm)•force; // project force onto tangent surface
50 :     dpos = hh*force;
51 :     if (|dpos| > stepMax) {
52 :     // decrease hh by factor by which step was too big
53 :     hh *= stepMax/|dpos|;
54 :     // and find smaller step
55 :     dpos = hh*force; // == stepMax*normalize(force);
56 :     }
57 :     vec2 posLast = pos;
58 :     // take step and re-apply constraint
59 :     pos = normalize(pos + dpos);
60 :     real energy = 0;
61 :     foreach (Particle p_j in sphere(rad)) {
62 :     energy += enr(p_j.pos - pos);
63 :     }
64 :     if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
65 :     hh *= 0.5; // backtrack
66 :     pos = posLast;
67 :     } else {
68 :     hh *= 1.1;
69 :     /* putting new particle creation here, to avoid possibility of
70 :     starting multiple particles at same position while, while not
71 :     moving due to backtracking */
72 :     if (ncount < 2) {
73 :     vec2 npos = pos + 0.2*rad*[norm[1],-norm[0]];
74 :     new Particle(ii, born, npos, hh);
75 :     born += 1;
76 :     }
77 :     }
78 :     outPos = pos;
79 :     }
80 :     }
81 :    
82 :     global {
83 :     /*
84 :     real mvmt = mean { |P.dpos|/rad | P in Particle.all};
85 :     if (mvmt < 0.0001) {
86 :     print("stabilized at iter ", iter, "\n");
87 :     stabilize;
88 :     }
89 :     */
90 :     iter+=1;
91 :     if (iter > iterMax) {
92 :     print("numActive() = ", numActive(), "\n");
93 :     stabilize; // HEY really want "die" here: didn't converge
94 :     }
95 :     // real var = variance{P.energy | P in Particle.all};
96 :     // real var = sum{(avg-P.energy)^2 | P in Particle.all};
97 :     }
98 :    
99 :     initially {Particle(-1, ii, initPos[ii], hhInit)
100 :     | ii in 0 .. numOfParticles-1 };

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