1 
/* 
/* ========================================== 
2 
running this under valgrind produces: 
## circle.diderot: Mutuallyrepelling 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 (postinitialization, before first iteration) snapshot 





*/ 

3 


4 

This example computes a system of mutuallyrepelling particles constrained to 
5 

the unit circle. The particles are initialized from a given sequence of a 
6 

2vectors, 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 coefficientofvariation of nearest 
13 

interparticle 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 unitlength 2vectors 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 postprocess with: 
43 


44 

export NRRD_STATE_VERBOSE_IO=0 
45 

for PIIN in pos????.nrrd; do 
46 

IIN=${PIIN#*} 
47 

II=${IIN%.*} 
48 

echo "postprocessing 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 [anglesref.png](anglesref.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 interparticle 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 

compiletime analysis, but there is a runtime check on the shape of the 
77 

array. A sequence of scalars should come from a 1D array, and a sequence 
78 

of any other kind of value should come from a 2D 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 nonscalar sequences. */ 
83 
input vec2{} ipos ("initial positions for all particles") = load("circlevec2.nrrd"); 
input vec2{} ipos ("initial positions for all particles") = load("circlevec2.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 coefficientofvariation of distances to nearest neighbors") = 0.05; 
input real eps ("system convergence threshold, computed as the coefficientofvariation 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 interparticle 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 

// Cotangentbased 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 2D 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 


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; 
185 
// save actual step taken 
// save actual step taken 
186 
step = pos  posLast; 
step = pos  posLast; 
187 


188 

/* The ArmijoGoldstein 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; 
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; 
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 }; 