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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4917, Wed Feb 15 20:36:49 2017 UTC revision 4918, Wed Feb 15 20:40:37 2017 UTC
# Line 1  Line 1 
1  /*  /* ==========================================
2  running this under valgrind produces:  ## circle.diderot: Mutually-repelling particles on a unit circle
 Syscall param write(buf) points to uninitialised byte(s)" that seems to be about the memory in the nrrd not being set  
 due to saving out the first (post-initialization, before first iteration) snapshot  
   
 */  
3    
4    This example computes a system of mutually-repelling particles constrained to
5    the unit circle. The particles are initialized from a given sequence of a
6    2-vectors, and each update tries to move each particle to a lower potential
7    energy, as determined by looking at nearby particles.  As in the
8    [`life.diderot`](../life) example that introduced strand communication, the
9    special strand variable named `pos` indicates the particle position, and the
10    computation is monitored via a `global` update method. Like
11    [`sieve.diderot`](../sieve), this program computes global reductions in the
12    global update, used here to compute the coefficient-of-variation of nearest
13    inter-particle distances, as the basis of a geometric test of convergence.
14    
15    Like the [`life.diderot`](../life) example, we use snapshots to inspect
16    the system state as it progresses, so we compile with:
17    
18            diderotc --snapshot --double --exec circle.diderot
19    
20    When compiling with `--double`, Diderot represents the `real` type
21    as a `double` in C++, instead of the default `float`, which produces
22    more accurate results.
23    
24    The following creates a list of `N` randomly oriented unit-length 2-vectors in
25    `vec2.nrrd` which will be the initial input.  Setting the random number seed
26    of `unu 1op nrand` with `-s RNG` ensures reproducible results, or you can
27    leave that out to give different results each time.
28    
29            N=60
30            RNG=42
31            eval echo {1..$[2*N]} | unu reshape -s 2 $N | unu 1op nrand -s $RNG -o vec2.nrrd
32            unu project -i vec2.nrrd -a 0 -m l2 | unu axinsert -a 0 -s 2 | unu 2op / vec2.nrrd - -o vec2.nrrd
33    
34    Then to run with snapshots saved every iteration (`-s 1`), but limiting the program
35    to 500 iterations (`-l 500`), as well as cleaning results from previous run:
36    
37            rm -f pos-????.nrrd pos.nrrd
38            ./circle -s 1 -l 500
39    
40    If this fails with `ERROR: unexpected arg (or unrecognized flag): "-s"`, it means that
41    `circle` wasn't compiled with `--snapshot`, as above.  It should converge in under
42    500 iterations, producing many `pos-????.nrrd` files, which we can post-process with:
43    
44            export NRRD_STATE_VERBOSE_IO=0
45            for PIIN in pos-????.nrrd; do
46               IIN=${PIIN#*-}
47               II=${IIN%.*}
48               echo "post-processing snapshot $II ... "
49               unu dice -i $PIIN -a 0 -o ./
50               unu 2op atan2 0.nrrd 1.nrrd | unu histax -a 0 -min -pi -max pi -b 800 -t float -o angle-$II.nrrd
51            done
52            rm -f 0.nrrd 1.nrrd
53            unu join -i angle-*.nrrd -a 1 |
54               unu quantize -b 8 -min 0 -max 1 -o angles.png
55            rm -f angle-????.nrrd
56    
57    This produces `angles.png` by unrolling (with atan2) positions along the circle,
58    and laying these out along horizontal scanlines, one per iteration. The result
59    should look something like [angles-ref.png](angles-ref.png). It is clear that the
60    particle interactions made a roughly uniform distribution early on, but
61    subsequent refinements took longer. Scrutinizing the top of the image shows
62    where strands were stationary either because they repeatedly backtracked in the
63    Armijo line search (until they found a step size that led to a sufficient
64    energy decrease), or because they had no other particles to interact with.
65    
66    Some experiments that can be tried with this example:
67    * If `--double` is removed from the compilation, the system may not converge to the default `eps` before the iteration limit, due to the decrease in numerical accuracy.
68    * There are different possible inter-particle potentials `phi` (uncomment one at a time), which can lead to different dynamics and different convergence speed towards a uniform distribution.
69    * If the interaction radius `rad` is decreased while maintining the number of particles, gaps may form.
70    * Increasing radius `rad` will tend to force convergence faster, at the expense of computing more particle interactions.
71    * There are constants associated with Armijo line search, changing them may lead to faster convergence.
72    
73    ========================================== */
74    
75    /* The file containing a sequence (read via "load") currently undergoes no
76       compile-time analysis, but there is a run-time check on the shape of the
77       array.  A sequence of scalars should come from a 1-D array, and a sequence
78       of any other kind of value should come from a 2-D array: the faster axis
79       should contain the components per value, and the sequence of values is
80       along the slower axis. The filename "-" may be used to say that the
81       sequence should be read from standard in.  A text file (one value per line)
82       may be used for non-scalar sequences. */
83  input vec2{} ipos ("initial positions for all particles") = load("circle-vec2.nrrd");  input vec2{} ipos ("initial positions for all particles") = load("circle-vec2.nrrd");
84  input real rad ("radius of particle's potential energy support") = 0.25;  input real rad ("radius of particle's potential energy support") = 0.25;
85  input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;  input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.01;
86    
87  real hhInit = 1;        // initial step size  real hhInit = 1;        // initial step size
88  real stepMax = rad;     // limit on distance to travel per iter  real stepMax = rad;     // limit on distance to travel per iter
89  int iter = 0;           // which iteration we're on  int iter = 0;           // which iteration we're on
90    
91    /* phi defines the potential energy around a particle, as a function of
92       distance from particle r > 0. By design, the potential function smoothly
93       falls to 0 at r == 1, so that the potential has compact support, which
94       simplifies computing the inter-particle interactions. phi'(r) is the
95       derivative of phi(r). The ' in phi' is merely suggestive; ' is not acting
96       as an operator. Below are some different possibilities for phi and phi';
97       only one should be uncommented.
98    
99       In a conservative force field, force is the *negative* gradient of some
100       potential. While phi(r) is the potential at r > 0, due to a particle at 0,
101       it can also be considered as the potential at 0, due to a particle at
102       r. Moving closer to the particle at r by a small dr > 0, the change in
103       potential is dphi = -phi'(r)*dr, which is positive for a monotonically
104       decreasing phi(r). The force at 0, or negative gradient of the potential,
105       is then phi'(r), which justifies how no negation is used in the frc()
106       definition below, based on phi'(). */
107    
108  // phi(0) and phi'(0) are bounded  // phi(0) and phi'(0) are bounded
109  function real  phi(real r) = (1 - r)^4;  function real  phi(real r) = (1 - r)^4;
110  function real phi'(real r) = -4*(1 - r)^3;  function real phi'(real r) = -4*(1 - r)^3;
111    
112    /*
113    // electrostatic potential 1/r, scaled to be C^2 continuous with 0 at r==1
114    function real  phi(real r) =  (1/r)*(1 - r)^3;
115    function real phi'(real r) = 3 - 1/(r^2) - 2*r;
116    */
117    /*
118    // Cotangent-based potential from Meyer et al. SMI'05
119    function real  phi(real r) = 1/tan(r*π/2) + r*π/2 - π/2;
120    function real phi'(real r) = (π/2)*(1 - (1/sin(r*π/2))^2);
121    */
122    
123    /* The enr and frc functions use phi and phi' to define the potential due to,
124       and force from, a particle at 2-D offset x with a potential field extending
125       to radius rad. */
126  function real enr(vec2 x) = phi(|x|/rad);  function real enr(vec2 x) = phi(|x|/rad);
127  function vec2 frc(vec2 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule  function vec2 frc(vec2 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule
128    
# Line 24  Line 130 
130  strand particle (vec2 pos0) {  strand particle (vec2 pos0) {
131     /* These strand variables are visible to the global update and to other     /* These strand variables are visible to the global update and to other
132        strands making spatial queries. Any variable inside the scope of        strands making spatial queries. Any variable inside the scope of
133        update{} will not be visible in this way. NOTE: "pos" is a special        update{} will not be visible in this way. NOTE: "pos" is the special
134        variable name that *must* be used to enable spatial queries of        variable name that *must* be used to enable spatial queries of
135        neighboring particles via sphere(). */        neighboring particles via sphere(). */
136     output vec2 pos = pos0/|pos0|;     output vec2 pos = pos0/|pos0|;
# Line 79  Line 185 
185        // save actual step taken        // save actual step taken
186        step = pos - posLast;        step = pos - posLast;
187    
188          /* The Armijo-Goldstein sufficient decrease condition: The potential
189             gradient, -force, dotted with change in position (pos - posLast),
190             should predict the change in potential (energy - energyLast), which
191             should be negative. If the potential change is too large (not
192             negative enough), decrease step size and try again on the next
193             iteration. */
194        if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {        if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
195           hh *= 0.5; // backtrack           hh *= 0.5; // backtrack
196           pos = posLast;           pos = posLast;
# Line 88  Line 200 
200     }     }
201  }  }
202    
203    /* In some programs, like ../heron/heron.diderot, strands can individually
204       determine when they have converged, and then execute "stabilize".  In this
205       kind of particle system, the judgement of convergence has to come from a
206       more global view of the state of computation.  This global update monitors
207       the particle system by measuring how uniformly the points are distributed
208       along the circle, based on the coefficient of variation of distances to
209       nearest neighbors. If the result has converged sufficiently, the
210       "stabilize" here results in all active strands being stabilized (and all
211       their stabilize{} methods would be called, if they had them). On the other
212       hand, if it is clear that the computation cannot converge, global update
213       can also call "die". */
214  global {  global {
215       // Convergence test could be based on movement
216       //real mvmt = max { |P.step|/rad | P in particle.all};
217       // This tests convergence based on distances to closest neighbor
218     real meancl = mean { P.closest | P in particle.all};     real meancl = mean { P.closest | P in particle.all};
219     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
220     real covcl = sqrt(varicl) / meancl;     real covcl = sqrt(varicl) / meancl;
# Line 101  Line 227 
227     iter += 1;     iter += 1;
228  }  }
229    
230    /* In this example, there is no "population control" (with "new" or "die"), so
231       the set of strands is constant throughout the program execution. It thus
232       doesn't matter if the strands are created via "initially {}" (a collection
233       of strands) or "initially []" (an array of strands). We create a collection
234       of strands here to be consistent with future particle system examples.*/
235  initially { particle(ipos[ii]) | ii in 0 .. length(ipos)-1 };  initially { particle(ipos[ii]) | ii in 0 .. length(ipos)-1 };

Legend:
Removed from v.4917  
changed lines
  Added in v.4918

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