sandbox/Antoonvh/distancecube.c
The distance between two point within a cube.
In order to test a double grid iterator implementation (i.e.nested) we aim to find the average distance between all points within a cube. The assumption here is that this number may be approximated using a finite number of points. The exact expression and a decimal expression is provided by mathworld, which will serve as a benchmark for our result.
#include "grid/octree.h"
int main(){
FILE * fp = fopen("distance","w");
int o = BGHOSTS+1;
We check for an increasing number of points the distance to all points, using the cell centered coordinates. We increase the resolution upto 6 levels of refinement. Since we use a rather naive nested iterator, the scaling with increasing resolution is horrible (six dimensional).
for (int g=1;g<6;g++){
int lev=g;
init_grid(1<<lev);
scalar xx[],yy[],zz[];
int n=0;
double dist=0;
Since we are not allowed to nest two foreach()
loops, it was decided to store the locations of the cell centeres in dedicated fields. I imagine persuiing a nested foreach()
loop would be quite difficult, especially when working with a decomposed domain (e.g. when using MPI).
foreach(){
xx[]=x; yy[]=y; zz[]=z;
}
In order to determine the distances, we iterate over each grid cell and store the location as xp,yp,zp.
foreach(){
double xp=x;
double yp=y;
double zp=z;
For each cell, we loop over all cells using a coding friendly Cartesian fashion and determine the corresponding distance. Which is then added to the total distance variable dist. The iterator should also be able to iterate over tree-grid-cells with a variable resolution.
for (int l =1;l<=depth();l++){
int Nl= (1<<l);
for (int i=o; i<Nl+o; i++){
for (int j=o; j<Nl+o; j++){
for (int k=o; k<Nl+o; k++){
Point point;
point.level=l; point.i=i; point.j=j; point.k=k;
if(is_active(cell) && is_local(cell)){ //Does this cell exist on this pid?
if(is_leaf(cell)){ //We only consider the leaf cells.
dist+=pow(sq(xp-xx[])+sq(yy[]-yp)+sq(zp-zz[]),0.5);
n++;
}
}
}
}
}
}
}
the result is written to a file, together with the error based on the analytical result.
fprintf(fp,"%d\t%g %d\t%d\t%g\n",g,dist/(double)n,n,N,fabs(dist/(double)n-0.6617071822));
}
}
Result
set xr [0.5:5.5]
set yr [0.5:0.7]
set ylabel 'Distance'
set xlabel 'Refinement Iteration'
plot 0.6617072 w lines lw 3 t 'Analytical results' ,\
'distance' u 1:2 t 'Grid Approximated average distance'
set xr [32:5000000000]
set yr [0.0001:0.2]
set logscale y
set logscale x
set ylabel 'Error'
set xlabel 'Number of pairs evaluated'
plot 'distance' u 3:5 t 'Error in the grid Approximated average distance' ,\
x**(-0.333) t '2/6-power scaling'
Notice that for each refinement iteration the required effort for the algorithm increases with a factor of (2^6)=64, compared to the previous iteration. For the modestly-sized 32^3 we have evaluated the distances over 1 billion pairs! also keeping in mind that the pow(x,0.5)
function is not particularly cheap. If only we could come up with an adaptive-grid approach…