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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4701 - (view) (download)

1 : glk 4701 #version 2
2 : glk 4640 // unit-circle
3 :     //
4 :     // Demo of distributing particles on the unit circle
5 :     //
6 : glk 4701 input vec2[] initPos = load_sequence("vec2.nrrd");
7 : glk 4640 //input vec2[] initPos = {[0.15606558, -0.56803983]}; // BUG
8 :     input int iterMax ("maximum number of iterations to run") = 1;
9 :     input real rad ("radius of particle's energy support") = 0.2;
10 :     input int pcp ("periodicity of population control") = 5;
11 :     input real mvmtEps ("system convergence criterion on particle movement ") = 0.001;
12 :    
13 :     real hhInit = 10; // initial step size, can be big
14 :     real newDist = 0.5*rad; // how far away to birth new particles
15 :     real stepMax = 2*rad; // limit on distance to travel per iter
16 :    
17 :     int iter = 0; // which iteration we're on
18 :    
19 :     // inter-particle energy, and its derivative
20 :     function real φ(real r) = (1 - r)^4;
21 :     function real φ'(real r) = -4*(1 - r)^3;
22 :     //function real φ(real r) = (1/r)*(1-r)^3;
23 :     //function real φ'(real r) = 3 - 1/(r^2) - 2*r;
24 :     function real enr(vec2 x) = φ(|x|/rad);
25 :     function vec2 frc(vec2 x) = φ'(|x|/rad) * (1/rad) * x/|x|; // chain rule
26 :    
27 :     strand Particle (int parent_, int ii, vec2 pos0, real hh0) {
28 :     int parent = parent_;
29 :     vec2 pos = pos0/|pos0|;
30 :     real hh = hh0;
31 :     output vec2 outPos = pos;
32 :     vec2 delta = [0,0];
33 :     int id = ii;
34 :     int born = 0;
35 :     real energyLast = ii; // HEY: can do convergence test based on variance of energies
36 :     vec2 force = [0,0]; // or can test convergence based on sum of |force|
37 :     bool first = true;
38 :     update {
39 :     if (first) {
40 :     //print("(iter ", iter, ") hello from ", ii, "\n");
41 :     first = false; // HEY remove me
42 :     }
43 :     // compute energy and forces on us
44 :     energyLast = 0;
45 :     force = [0,0];
46 :     int ncount = 0;
47 :     foreach (Particle p_j in sphere(rad)) {
48 :     ncount += 1;
49 :     vec2 r_ij = p_j.pos - pos;
50 :     energyLast += enr(r_ij);
51 :     force += frc(r_ij);
52 :     }
53 :     vec2 norm = normalize(pos); // surface normal (for unit circle)
54 :     if (ncount == 0 && (pcp > 0 && 0 == iter % pcp)) {
55 :     vec2 npos = pos + newDist*rad*[norm[1],-norm[0]];
56 :     new Particle(ii, born, npos, hh);
57 :     born += 1;
58 :     continue;
59 :     }
60 :     // else update position based on force from neighbors.
61 :     // first, project force onto tangent surface
62 :     force = (identity[2] - norm⊗norm)•force;
63 :     delta = hh*force;
64 :     if (|delta| > stepMax) {
65 :     // decrease hh by factor by which step was too big
66 :     hh *= stepMax/|delta|;
67 :     // and find smaller step
68 :     delta = hh*force; // == stepMax*normalize(force);
69 :     }
70 :     vec2 posLast = pos;
71 :     // take step and re-apply constraint
72 :     pos = normalize(pos + delta);
73 :     real energy = 0;
74 :     ncount = 0;
75 :     foreach (Particle p_j in sphere(rad)) {
76 :     energy += enr(p_j.pos - pos);
77 :     ncount += 1;
78 :     }
79 :     if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
80 :     hh *= 0.5; // backtrack
81 :     pos = posLast;
82 :     } else {
83 :     hh *= 1.1;
84 :     /* putting new particle creation here avoids possibility of
85 :     starting multiple particles at same position while not
86 :     moving due to backtracking */
87 :     if (ncount < 2 && (pcp > 0 && 0 == iter % pcp)) {
88 :     vec2 npos = pos + newDist*rad*normalize(force);
89 :     new Particle(ii, born, npos, hh);
90 :     born += 1;
91 :     }
92 :     }
93 :     outPos = pos;
94 :     }
95 :     }
96 :    
97 : glk 4701 update {
98 : glk 4640 //print("(iter ", iter, ") hello from global\n");
99 :     real mvmt = max { |P.delta|/rad | P in Particle.all};
100 : glk 4648 if (numActive() > 1 && mvmt < mvmtEps) {
101 : glk 4640 print("stabilizing with movment ", mvmt, " at iter ", iter, "\n");
102 :     stabilize;
103 :     }
104 :     iter += 1;
105 :     if (iter >= iterMax) {
106 :     print("stopping with movment ", mvmt, " at iter ", iter, "\n");
107 :     die;
108 :     }
109 :     }
110 :    
111 : glk 4701 create_collection {Particle(-1, ii, initPos[ii], hhInit)
112 :     | ii in 0 .. length(initPos)-1 }

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