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

SCM Repository

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

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

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