Skip to content
GitLab
Commit 4434bc1d authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

fixed bug in cell split. 


Former-commit-id: 2f3e816cab93990d8d69e44ba478bd61ce67976e
parent 22d1a432
Hide whitespace changes
Inline Side-by-side
runner.c
View file @ 4434bc1d
......@@ -371,7 +371,7 @@ void runner_dopair ( struct runner_thread *rt , struct cell *ci , struct cell *c
ri_max = ci->r_max - rshift; rj_max = cj->r_max - rshift;
count_i = ci->count; count_j = cj->count;
parts_i = ci->parts; parts_j = cj->parts;
di_max = sort_i[count_i-1].d;
di_max = sort_i[count_i-1].d - rshift;
dj_min = sort_j[0].d;
/* if ( ci->split && cj->split && sid == 4 )
......
space.c
View file @ 4434bc1d
......@@ -410,7 +410,8 @@ void space_splittasks ( struct space *s ) {
break;
case 1: /* ( 1 , 1 , 0 ) */
if ( !ci->progeny[6]->split && !ci->progeny[7]->split &&
if ( space_dosub &&
!ci->progeny[6]->split && !ci->progeny[7]->split &&
!cj->progeny[0]->split && !cj->progeny[1]->split &&
ci->progeny[6]->count + ci->progeny[7]->count + cj->progeny[0]->count + cj->progeny[1]->count < space_splitsize ) {
t->type = tid_sub; t->flags = 1;
......@@ -443,7 +444,8 @@ void space_splittasks ( struct space *s ) {
break;
case 3: /* ( 1 , 0 , 1 ) */
if ( !ci->progeny[5]->split && !ci->progeny[7]->split &&
if ( space_dosub &&
!ci->progeny[5]->split && !ci->progeny[7]->split &&
!cj->progeny[0]->split && !cj->progeny[2]->split &&
ci->progeny[5]->count + ci->progeny[7]->count + cj->progeny[0]->count + cj->progeny[2]->count < space_splitsize ) {
t->type = tid_sub; t->flags = 3;
......@@ -469,7 +471,8 @@ void space_splittasks ( struct space *s ) {
break;
case 4: /* ( 1 , 0 , 0 ) */
if ( !ci->progeny[4]->split && !ci->progeny[5]->split && !ci->progeny[6]->split && !ci->progeny[7]->split &&
if ( space_dosub &&
!ci->progeny[4]->split && !ci->progeny[5]->split && !ci->progeny[6]->split && !ci->progeny[7]->split &&
!cj->progeny[0]->split && !cj->progeny[1]->split && !cj->progeny[2]->split && !cj->progeny[3]->split &&
ci->progeny[4]->count + ci->progeny[5]->count + ci->progeny[6]->count + ci->progeny[7]->count + cj->progeny[0]->count + cj->progeny[1]->count + cj->progeny[2]->count + cj->progeny[3]->count < space_splitsize ) {
t->type = tid_sub; t->flags = 4;
......@@ -535,7 +538,8 @@ void space_splittasks ( struct space *s ) {
break;
case 5: /* ( 1 , 0 , -1 ) */
if ( !ci->progeny[4]->split && !ci->progeny[6]->split &&
if ( space_dosub &&
!ci->progeny[4]->split && !ci->progeny[6]->split &&
!cj->progeny[1]->split && !cj->progeny[3]->split &&
ci->progeny[4]->count + ci->progeny[6]->count + cj->progeny[1]->count + cj->progeny[3]->count < space_splitsize ) {
t->type = tid_sub; t->flags = 5;
......@@ -568,7 +572,8 @@ void space_splittasks ( struct space *s ) {
break;
case 7: /* ( 1 , -1 , 0 ) */
if ( !ci->progeny[4]->split && !ci->progeny[5]->split &&
if ( space_dosub &&
!ci->progeny[4]->split && !ci->progeny[5]->split &&
!cj->progeny[2]->split && !cj->progeny[3]->split &&
ci->progeny[4]->count + ci->progeny[5]->count + cj->progeny[2]->count + cj->progeny[3]->count < space_splitsize ) {
t->type = tid_sub; t->flags = 7;
......@@ -601,7 +606,8 @@ void space_splittasks ( struct space *s ) {
break;
case 9: /* ( 0 , 1 , 1 ) */
if ( !ci->progeny[3]->split && !ci->progeny[7]->split &&
if ( space_dosub &&
!ci->progeny[3]->split && !ci->progeny[7]->split &&
!cj->progeny[0]->split && !cj->progeny[4]->split &&
ci->progeny[3]->count + ci->progeny[7]->count + cj->progeny[0]->count + cj->progeny[4]->count < space_splitsize ) {
t->type = tid_sub; t->flags = 9;
......@@ -693,7 +699,8 @@ void space_splittasks ( struct space *s ) {
break;
case 11: /* ( 0 , 1 , -1 ) */
if ( !ci->progeny[2]->split && !ci->progeny[6]->split &&
if ( space_dosub &&
!ci->progeny[2]->split && !ci->progeny[6]->split &&
!cj->progeny[1]->split && !cj->progeny[5]->split &&
ci->progeny[2]->count + ci->progeny[6]->count + cj->progeny[1]->count + cj->progeny[5]->count < space_splitsize ) {
t->type = tid_sub; t->flags = 11;
......@@ -719,7 +726,8 @@ void space_splittasks ( struct space *s ) {
break;
case 12: /* ( 0 , 0 , 1 ) */
if ( !ci->progeny[1]->split && !ci->progeny[3]->split && !ci->progeny[5]->split && !ci->progeny[7]->split &&
if ( space_dosub &&
!ci->progeny[1]->split && !ci->progeny[3]->split && !ci->progeny[5]->split && !ci->progeny[7]->split &&
!cj->progeny[0]->split && !cj->progeny[2]->split && !cj->progeny[4]->split && !cj->progeny[6]->split &&
ci->progeny[1]->count + ci->progeny[3]->count + ci->progeny[5]->count + ci->progeny[7]->count + cj->progeny[0]->count + cj->progeny[2]->count + cj->progeny[4]->count + cj->progeny[6]->count < space_splitsize ) {
t->type = tid_sub; t->flags = 12;
......@@ -883,7 +891,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
maketasks_rec( c->progeny[k] , sort , nr_sort , c );
/* Worth splitting into several tasks? */
if ( c->count > 1.5*space_splitsize ) {
if ( !space_dosub || c->count > 1.5*space_splitsize ) {
/* Make a task for eac pair of progeny. */
for ( j = 0 ; j < 8 ; j++ )
......@@ -1012,7 +1020,7 @@ void space_maketasks ( struct space *s , int do_sort ) {
void cell_split ( struct cell *c ) {
int i, j, k;
int i, j, k, kk;
struct part temp, *parts = c->parts;
int left[8], right[8];
double pivot[3];
......@@ -1023,30 +1031,44 @@ void cell_split ( struct cell *c ) {
/* Split along the x-axis. */
i = 0; j = c->count - 1;
while ( i < j ) {
while ( i < c->count-1 && parts[i].x[0] <= pivot[0] )
while ( i <= j ) {
while ( i <= c->count-1 && parts[i].x[0] <= pivot[0] )
i += 1;
while ( j > 0 && parts[j].x[0] > pivot[0] )
while ( j >= 0 && parts[j].x[0] > pivot[0] )
j -= 1;
if ( i < j ) {
temp = parts[i]; parts[i] = parts[j]; parts[j] = temp;
}
}
for ( k = 0 ; k <= j ; k++ )
if ( parts[k].x[0] > pivot[0] )
error( "cell_split: sorting failed." );
for ( k = i ; k < c->count ; k++ )
if ( parts[k].x[0] < pivot[0] )
error( "cell_split: sorting failed." );
left[1] = i; right[1] = c->count - 1;
left[0] = 0; right[0] = j;
/* Split along the y axis, twice. */
for ( k = 1 ; k >= 0 ; k-- ) {
i = left[k]; j = right[k];
while ( i < j ) {
while ( i < right[k] && parts[i].x[1] <= pivot[1] )
while ( i <= j ) {
while ( i <= right[k] && parts[i].x[1] <= pivot[1] )
i += 1;
while ( j > left[k] && parts[j].x[1] > pivot[1] )
while ( j >= left[k] && parts[j].x[1] > pivot[1] )
j -= 1;
if ( i < j ) {
temp = parts[i]; parts[i] = parts[j]; parts[j] = temp;
}
}
for ( kk = left[k] ; kk <= j ; kk++ )
if ( parts[kk].x[1] > pivot[1] ) {
printf( "cell_split: ival=[%i,%i], i=%i, j=%i.\n" , left[k] , right[k] , i , j );
error( "sorting failed (left)." );
}
for ( kk = i ; kk <= right[k] ; kk++ )
if ( parts[kk].x[1] < pivot[1] )
error( "sorting failed (right)." );
left[2*k+1] = i; right[2*k+1] = right[k];
left[2*k] = left[k]; right[2*k] = j;
}
......@@ -1054,15 +1076,25 @@ void cell_split ( struct cell *c ) {
/* Split along the z axis, four times. */
for ( k = 3 ; k >= 0 ; k-- ) {
i = left[k]; j = right[k];
while ( i < j ) {
while ( i < right[k] && parts[i].x[2] <= pivot[2] )
while ( i <= j ) {
while ( i <= right[k] && parts[i].x[2] <= pivot[2] )
i += 1;
while ( j > left[k] && parts[j].x[2] > pivot[2] )
while ( j >= left[k] && parts[j].x[2] > pivot[2] )
j -= 1;
if ( i < j ) {
temp = parts[i]; parts[i] = parts[j]; parts[j] = temp;
}
}
for ( kk = left[k] ; kk <= j ; kk++ )
if ( parts[kk].x[2] > pivot[2] ) {
printf( "cell_split: ival=[%i,%i], i=%i, j=%i.\n" , left[k] , right[k] , i , j );
error( "sorting failed (left)." );
}
for ( kk = i ; kk <= right[k] ; kk++ )
if ( parts[kk].x[2] < pivot[2] ) {
printf( "cell_split: ival=[%i,%i], i=%i, j=%i.\n" , left[k] , right[k] , i , j );
error( "sorting failed (right)." );
}
left[2*k+1] = i; right[2*k+1] = right[k];
left[2*k] = left[k]; right[2*k] = j;
}
......@@ -1315,11 +1347,6 @@ void space_init ( struct space *s , double dim[3] , struct part *parts , int N ,
c->count += 1;
}
/* Loop over the cells and split them. */
/* for ( k = 0 ; k < nr_cells ; k++ )
if ( cells[k].count > space_maxppc )
space_split( s , &cells[k] ); */
/* Store eveything in the space. */
s->r_min = r_min; s->r_max = r_max;
s->dim[0] = dim[0]; s->dim[1] = dim[1]; s->dim[2] = dim[2];
......@@ -1338,6 +1365,6 @@ void space_init ( struct space *s , double dim[3] , struct part *parts , int N ,
/* Loop over the cells and split them. */
for ( k = 0 ; k < nr_cells ; k++ )
space_split( s , &cells[k] );
}
space.h
View file @ 4434bc1d
......@@ -25,6 +25,7 @@
#define space_cellallocchunk 1000
#define space_splitratio 0.875
#define space_splitsize_default 800
#define space_dosub 0
#define task_maxwait 3
#define task_maxunlock 39
......
test.c
View file @ 4434bc1d
......@@ -117,6 +117,32 @@ void map_cells_plot ( struct cell *c , void *data ) {
* @brief Mapping function for neighbour count.
*/
void map_cellcheck ( struct cell *c , void *data ) {
int k, *count = (int *)data;
struct part *p;
/* Loop over all parts and check if they are in the cell. */
for ( k = 0 ; k < c->count ; k++ ) {
*count += 1;
p = &c->parts[k];
if ( p->x[0] < c->loc[0] || p->x[1] < c->loc[1] || p->x[2] < c->loc[2] ||
p->x[0] > c->loc[0] + c->h[0] || p->x[1] > c->loc[1] + c->h[1] || p->x[2] > c->loc[2] + c->h[2] ) {
printf( "map_cellcheck: particle at [ %e %e %e ] outside of cell [ %e %e %e ] - [ %e %e %e ].\n" ,
p->x[0] , p->x[1] , p->x[2] ,
c->loc[0] , c->loc[1] , c->loc[2] ,
c->loc[0] + c->h[0] , c->loc[1] + c->h[1] , c->loc[2] + c->h[2] );
error( "particle out of bounds!" );
}
}
}
/**
* @brief Mapping function for maxdepth cell count.
*/
void map_maxdepth ( struct cell *c , void *data ) {
int maxdepth = ((int *)data)[0];
......@@ -374,11 +400,11 @@ void read_dt ( char *fname , struct part *parts , int N ) {
void pairs_n2 ( double *dim , struct part *__restrict__ parts , int N , int periodic ) {
int i, j, k, count = 0, mj, mk;
double r, dx[3], dcount = 0.0, maxratio = 1.0;
double r2, dx[3], dcount = 0.0, maxratio = 1.0;
float f1, f2;
/* Loop over all particle pairs. */
#pragma omp parallel for schedule(dynamic), default(none), private(k,i,dx,r,f1,f2), shared(maxratio,mj,mk,periodic,parts,dim,N,stdout), reduction(+:count,dcount)
#pragma omp parallel for schedule(dynamic), default(none), private(k,i,dx,r2,f1,f2), shared(maxratio,mj,mk,periodic,parts,dim,N,stdout), reduction(+:count,dcount)
for ( j = 0 ; j < N ; j++ ) {
if ( j % 1000 == 0 ) {
printf( "pairs_n2: j=%i, count=%i.\n" , j , count );
......@@ -394,10 +420,10 @@ void pairs_n2 ( double *dim , struct part *__restrict__ parts , int N , int peri
dx[i] -= dim[i];
}
}
r = sqrt( dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2] );
if ( r < parts[j].r || r < parts[k].r ) {
r2 = dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2];
if ( r2 < parts[j].r*parts[j].r || r2 < parts[k].r*parts[k].r ) {
f1 = 0.0; f2 = 0.0;
iact( r*r , parts[j].r , parts[k].r , &f1 , &f2 , &count , &count );
iact( r2 , parts[j].r , parts[k].r , &f1 , &f2 , &count , &count );
dcount += f1 + f2;
if ( parts[j].r / parts[k].r > maxratio )
#pragma omp critical
......@@ -579,7 +605,7 @@ int main ( int argc , char *argv[] ) {
}
/* Get the brute-force number of pairs. */
// pairs_n2( dim , parts , N , periodic );
pairs_n2( dim , parts , N , periodic );
// pairs_single( dim , parts , N , periodic , 63628 );
fflush( stdout );
......@@ -600,6 +626,11 @@ int main ( int argc , char *argv[] ) {
printf( "main: maximum depth is %d.\n" , s.maxdepth );
printf( "main: cutoffs in [ %.3f %.3f ].\n" , s.r_min , s.r_max ); fflush(stdout);
/* Verify that each particle is in it's propper cell. */
icount = 0;
space_map_cells( &s , &map_cellcheck , &icount );
printf( "main: map_cellcheck picked up %i parts.\n" , icount );
data[0] = s.maxdepth; data[1] = 0;
space_map_cells( &s , &map_maxdepth , data );
printf( "main: nr of cells at depth %i is %i.\n" , data[0] , data[1] );
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment