Skip to content
Snippets Groups Projects
Commit 5802a4e7 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Correct version of the task-based tree walk. Works for all numbers of... 

Correct version of the task-based tree walk. Works for all numbers of particles per cell and particles per task.
parent 6ec0120a
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
examples/plot.py
+ 18
− 17
View file @ 5802a4e7
...@@ -93,21 +93,20 @@ errz_new = (accz_new - accz_e )/abs(accz_e) ...@@ -93,21 +93,20 @@ errz_new = (accz_new - accz_e )/abs(accz_e)
# erry_g = (accy_g - accy_e )/abs(accy_e) # erry_g = (accy_g - accy_e )/abs(accy_e)
# errz_g = (accz_g - accz_e )/abs(accz_e) # errz_g = (accz_g - accz_e )/abs(accz_e)
# Statistics # Statistics
meanx_bh = mean(errx_bh) meanx_bh = mean(errx_bh[abs(errx_bh) < 0.1])
stdx_bh = sqrt(mean(errx_bh**2) - meanx_bh**2) stdx_bh = std(errx_bh[abs(errx_bh) < 0.1])
meany_bh = mean(erry_bh) meany_bh = mean(erry_bh[abs(erry_bh) < 0.1])
stdy_bh = std(erry_bh) stdy_bh = std(erry_bh[abs(erry_bh) < 0.1])
meanz_bh = mean(errz_bh) meanz_bh = mean(errz_bh[abs(errz_bh) < 0.1])
stdz_bh = std(errz_bh) stdz_bh = std(errz_bh[abs(errz_bh) < 0.1])
meanx_new = mean(errx_new) meanx_new = mean(errx_new[abs(errx_new) < 0.1])
stdx_new = std(errx_new) stdx_new = std(errx_new[abs(errx_new) < 0.1])
meany_new = mean(erry_new) meany_new = mean(erry_new[abs(erry_new) < 0.1])
stdy_new = std(erry_new) stdy_new = std(erry_new[abs(erry_new) < 0.1])
meanz_new = mean(errz_new) meanz_new = mean(errz_new[abs(errz_new) < 0.1])
stdz_new = std(errz_new) stdz_new = std(errz_new[abs(errz_new) < 0.1])
# meanx_g = mean(errx_g) # meanx_g = mean(errx_g)
# stdx_g = std(errx_g) # stdx_g = std(errx_g)
...@@ -123,7 +122,7 @@ subplot(311, title="Acceleration along X") ...@@ -123,7 +122,7 @@ subplot(311, title="Acceleration along X")
#plot(id, errx_g , 'gs') #plot(id, errx_g , 'gs')
plot(id, errx_bh , 'rx') plot(id, errx_bh , 'rx')
plot(id, errx_new , 'b.') plot(id, errx_new , 'b.')
#text(id[-1], 0.18, "Gadget: $%5.3f\\pm%5.3f$ \n B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanx_g, stdx_g, meanx_bh, stdx_bh, meanx_new, stdx_new), backgroundcolor="w", va="top", ha="right" ) text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanx_bh, stdx_bh, meanx_new, stdx_new), backgroundcolor="w", va="top", ha="right" )
ylim(-0.2, 0.2) ylim(-0.2, 0.2)
...@@ -134,7 +133,8 @@ subplot(312, title="Acceleration along Y") ...@@ -134,7 +133,8 @@ subplot(312, title="Acceleration along Y")
#plot(id, erry_g , 'gs') #plot(id, erry_g , 'gs')
plot(id, erry_bh , 'rx') plot(id, erry_bh , 'rx')
plot(id, erry_new , 'b.') plot(id, erry_new , 'b.')
#text(id[-1], 0.18, "Gadget: $%5.3f\\pm%5.3f$ \n B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_g, stdy_g, meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" ) text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" )
ylim(-0.2, 0.2) ylim(-0.2, 0.2)
xlim(0,id[-1]) xlim(0,id[-1])
...@@ -145,7 +145,8 @@ subplot(313, title="Acceleration along Z") ...@@ -145,7 +145,8 @@ subplot(313, title="Acceleration along Z")
#plot(id, errz_g , 'gs') #plot(id, errz_g , 'gs')
plot(id, errz_bh , 'rx') plot(id, errz_bh , 'rx')
plot(id, errz_new , 'b.') plot(id, errz_new , 'b.')
#text(id[-1], 0.18, "Gadget: $%5.3f\\pm%5.3f$ \n B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_g, stdy_g, meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" ) text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanz_bh, stdz_bh, meanz_new, stdz_new), backgroundcolor="w", va="top", ha="right" )
ylim(-0.2, 0.2) ylim(-0.2, 0.2)
xlim(0,id[-1]) xlim(0,id[-1])
......
examples/test_bh_3.c
+ 32
− 17
View file @ 5802a4e7
...@@ -36,11 +36,10 @@ ...@@ -36,11 +36,10 @@
/* Some local constants. */ /* Some local constants. */
#define cell_pool_grow 1000 #define cell_pool_grow 1000
#define cell_maxparts 1 #define cell_maxparts 1
#define task_limit 500 #define task_limit 5000
#define const_G 6.6738e-8 #define const_G 6.6738e-8
#define dist_min 0.5 // 0.5 #define dist_min 0.5 // 0.5
#define ICHECK -1 #define ICHECK -1
/** Data structure for the particles. */ /** Data structure for the particles. */
...@@ -65,18 +64,25 @@ struct cell { ...@@ -65,18 +64,25 @@ struct cell {
struct cell* firstchild; /* Next node if opening */ struct cell* firstchild; /* Next node if opening */
struct cell* sibling; /* Next node */ struct cell* sibling; /* Next node */
union { /* We keep both CoMs and masses to make sure the comp_com calculation is correct */
union {
/* Information for the legacy walk */
struct { struct {
double com[3]; double com[3];
double mass; double mass;
} legacy; } legacy;
/* Information for the QuickShed walk */
struct { struct {
double com[3]; double com[3];
double mass; double mass;
} new; } new;
} u; } u;
int res, com_tid; int res, com_tid;
} __attribute__((aligned (128))); } __attribute__((aligned (128)));
...@@ -404,7 +410,7 @@ void iact_pair_pc ( struct cell *ci , struct cell *cj ) { ...@@ -404,7 +410,7 @@ void iact_pair_pc ( struct cell *ci , struct cell *cj ) {
printf("%e %e %e %d %p %d %p\n", cj->u.new.com[0], cj->u.new.com[1], cj->u.new.com[2], cj->count, cj , cj->split, cj->sibling); printf("%e %e %e %d %p %d %p\n", cj->u.new.com[0], cj->u.new.com[1], cj->u.new.com[2], cj->count, cj , cj->split, cj->sibling);
for ( j = 0 ; j < cj->count ; ++j ) for ( j = 0 ; j < cj->count ; ++j )
printf("part %d mass= %e id= %d\n", j, cj->parts[j].mass , cj->parts[j].id ); printf("part %d mass= %e id= %d\n", j, cj->parts[j].mass , cj->parts[j].id );
error( "com does not seem to have been set." ); error( "com does not seem to have been set." );
} }
...@@ -456,7 +462,7 @@ void iact_pair ( struct cell *ci , struct cell *cj ) { ...@@ -456,7 +462,7 @@ void iact_pair ( struct cell *ci , struct cell *cj ) {
int i, j, k; int i, j, k;
int count_i = ci->count, count_j = cj->count; int count_i = ci->count, count_j = cj->count;
double dx[3], xi[3], ai[3], mi, mj, r2, w, ir; double dx[3], xi[3], ai[3], mi, mj, r2, r2_i, r2_j, w, ir;
struct part *parts_i = ci->parts, *parts_j = cj->parts; struct part *parts_i = ci->parts, *parts_j = cj->parts;
struct cell *cp; struct cell *cp;
...@@ -470,15 +476,17 @@ void iact_pair ( struct cell *ci , struct cell *cj ) { ...@@ -470,15 +476,17 @@ void iact_pair ( struct cell *ci , struct cell *cj ) {
/* Distance between the CoMs */ /* Distance between the CoMs */
for ( r2 = 0.0, k = 0 ; k < 3 ; k++ ) { for ( r2 = 0.0, r2_i = 0.0, r2_j = 0.0, k = 0 ; k < 3 ; k++ ) {
dx[k] = fabs( ci->u.new.com[k] - cj->u.new.com[k] ); // dx[k] = fabs( ci->u.new.com[k] - cj->u.new.com[k] );
r2 += dx[k]*dx[k]; dx[k] = fabs( ci->loc[k] - cj->loc[k] );
r2 += dx[k] * dx[k];
r2_i += ( dx[k] - 0.5*ci->h ) * ( dx[k] - 0.5*ci->h );
r2_j += ( dx[k] - 0.5*cj->h ) * ( dx[k] - 0.5*cj->h );
} }
double s_max_i = ci->h; if ( ( dist_min * dist_min * r2_j > ci->h * ci->h ) && ( dist_min * dist_min * r2_i > cj->h * cj->h ) )
double s_max_j = cj->h;
if ( ( dist_min * dist_min * r2 > s_max_i * s_max_i ) && ( dist_min * dist_min * r2 > s_max_j * s_max_j ) )
{ {
iact_pair_pc( ci, cj ); iact_pair_pc( ci, cj );
iact_pair_pc( cj, ci ); iact_pair_pc( cj, ci );
...@@ -585,6 +593,7 @@ void iact_pair ( struct cell *ci , struct cell *cj ) { ...@@ -585,6 +593,7 @@ void iact_pair ( struct cell *ci , struct cell *cj ) {
cp = cp->sibling; cp = cp->sibling;
} }
} }
else else
error("Want to split unpslitable cells !\n"); error("Want to split unpslitable cells !\n");
...@@ -703,7 +712,7 @@ void create_tasks ( struct qsched *s , struct cell *ci , struct cell *cj ) { ...@@ -703,7 +712,7 @@ void create_tasks ( struct qsched *s , struct cell *ci , struct cell *cj ) {
int k; int k;
qsched_task_t tid; qsched_task_t tid;
struct cell *data[2], *cp, *cps; struct cell *data[2], *cp, *cps;
double dx, r2; double dx, r2, r2_i, r2_j;
/* If either cell is empty, stop. */ /* If either cell is empty, stop. */
...@@ -761,13 +770,17 @@ void create_tasks ( struct qsched *s , struct cell *ci , struct cell *cj ) { ...@@ -761,13 +770,17 @@ void create_tasks ( struct qsched *s , struct cell *ci , struct cell *cj ) {
else { else {
/* Distance between the cells */ /* Distance between the cells */
for ( r2 = 0.0, k = 0 ; k < 3 ; k++ ) { for ( r2 = 0.0, r2_i = 0.0, r2_j = 0.0, k = 0 ; k < 3 ; k++ ) {
dx = fabs( ci->loc[k] - cj->loc[k] ); dx = fabs( ci->loc[k] - cj->loc[k] );
r2 += dx*dx; r2 += dx*dx;
r2_i += ( dx - 0.5 * ci->h ) * ( dx - 0.5 * ci->h );
r2_j += ( dx - 0.5 * cj->h ) * ( dx - 0.5 * cj->h );
} }
/* Check whether we can use the multipoles. */ /* Check whether we can use the multipoles. */
if ( ( dist_min * dist_min * r2 > ci->h * ci->h ) && ( dist_min * dist_min * r2 > cj->h * cj->h ) ) if ( ( dist_min * dist_min * r2_j > ci->h * ci->h ) && ( dist_min * dist_min * r2_i > cj->h * cj->h ) )
{ {
data[0] = ci; data[1] = cj; data[0] = ci; data[1] = cj;
tid = qsched_addtask( s , task_type_pair_pc , task_flag_none , data , sizeof(struct cell *) * 2 , ci->count ); tid = qsched_addtask( s , task_type_pair_pc , task_flag_none , data , sizeof(struct cell *) * 2 , ci->count );
...@@ -819,6 +832,7 @@ void create_tasks ( struct qsched *s , struct cell *ci , struct cell *cj ) { ...@@ -819,6 +832,7 @@ void create_tasks ( struct qsched *s , struct cell *ci , struct cell *cj ) {
cp = cp->sibling; cp = cp->sibling;
} }
} }
else else
error("Want to split unpslitable cells !\n"); error("Want to split unpslitable cells !\n");
...@@ -840,6 +854,7 @@ void create_tasks ( struct qsched *s , struct cell *ci , struct cell *cj ) { ...@@ -840,6 +854,7 @@ void create_tasks ( struct qsched *s , struct cell *ci , struct cell *cj ) {
cp = cp->sibling; cp = cp->sibling;
} }
} }
else else
error("Want to split unpslitable cells !\n"); error("Want to split unpslitable cells !\n");
...@@ -1289,7 +1304,7 @@ void test_bh ( int N , int nr_threads , int runs , char* fileName ) { ...@@ -1289,7 +1304,7 @@ void test_bh ( int N , int nr_threads , int runs , char* fileName ) {
/* Do a N^2 interactions calculation */ /* Do a N^2 interactions calculation */
tic_exact = getticks(); tic_exact = getticks();
interact_exact( N , parts , ICHECK ); //interact_exact( N , parts , ICHECK );
toc_exact = getticks(); toc_exact = getticks();
printf( "Exact calculation (1 thread) took %lli (= %e) ticks\n", toc_exact - tic_exact , (float)(toc_exact - tic_exact) ); printf( "Exact calculation (1 thread) took %lli (= %e) ticks\n", toc_exact - tic_exact , (float)(toc_exact - tic_exact) );
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment