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

FMM is now correct. Downpass still not task-based.

parent d6f4208a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
examples/plot.py
+ 3
− 3
View file @ d9cba760
......@@ -127,7 +127,7 @@ plot(id, errx_new , 'b.')
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.02, 0.02)
xlim(0,id[-1])
grid()
......@@ -138,7 +138,7 @@ plot(id, erry_new , 'b.')
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.02, 0.02)
xlim(0,id[-1])
grid()
......@@ -150,7 +150,7 @@ plot(id, errz_bh , 'rx', label="Legacy")
#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" )
legend(loc="upper right")
ylim(-0.2, 0.2)
ylim(-0.02, 0.02)
xlim(0,id[-1])
grid()
......
examples/test_fmm.c
+ 33
− 46
View file @ d9cba760
......@@ -37,7 +37,7 @@
/* Some local constants. */
#define cell_pool_grow 1000
#define cell_maxparts 1
#define cell_maxparts 128
#define task_limit 1
#define const_G 1 // 6.6738e-8
#define dist_min 0.5 /* Used for legacy walk only */
......@@ -255,19 +255,15 @@ void comp_down(struct cell *c) {
double dx=0.f, dy=0.f, dz=0.f;
struct part *parts = c->parts;
if (c->split) {
//message("Downpass: %f", c->field_x.F_000);
/* Loop over the siblings and propagate accelerations. */
for (cp = c->firstchild; cp != c->sibling; cp = cp->sibling) {
/* Distance over which to shift */
dx = cp->new.M_100 - c->new.M_100;
dy = cp->new.M_100 - c->new.M_100;
dz = cp->new.M_100 - c->new.M_100;
dy = cp->new.M_010 - c->new.M_010;
dz = cp->new.M_001 - c->new.M_001;
/* Propagate the tensors downwards */
shiftAndAddTensor(c->field_x, &cp->field_x, dx, dy, dz);
......@@ -294,10 +290,6 @@ void comp_down(struct cell *c) {
ay = const_G * applyFieldAcceleration(c->field_y, dx, dy, dz);
az = const_G * applyFieldAcceleration(c->field_z, dx, dy, dz);
//message("Apply to particle: %f %f", ax, c->field_x.F_000);
parts[k].a[0] += ax;
parts[k].a[1] += ay;
parts[k].a[2] += az;
......@@ -569,12 +561,6 @@ static inline void iact_pair_pc(struct cell *ci, struct cell *cj) {
computeFieldTensors_x( -r_x, -r_y, -r_z, inv_r, cp->new, &cps->field_x );
computeFieldTensors_y( -r_x, -r_y, -r_z, inv_r, cp->new, &cps->field_y );
computeFieldTensors_z( -r_x, -r_y, -r_z, inv_r, cp->new, &cps->field_z );
/* message("r= %f %f %f %f M_000= %f %f F_000= %f", r_x, r_y, r_z, sqrt(r2), */
/* cps->new.M_000, */
/* cps->new.M_000 * D_100( r_x, r_y, r_z, inv_r), */
/* cp->field_x.F_000); */
}
}
}
......@@ -1260,36 +1246,36 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
/* If no input file was specified, generate random particle positions. */
if (fileName == NULL || fileName[0] == 0) {
/* for (k = 0; k < N; k++) { */
/* parts[k].id = k; */
/* parts[k].x[0] = ((double)rand()) / RAND_MAX; */
/* parts[k].x[1] = ((double)rand()) / RAND_MAX; */
/* parts[k].x[2] = ((double)rand()) / RAND_MAX; */
/* parts[k].mass = 1.f;//((double)rand()) / RAND_MAX; */
/* parts[k].a_legacy[0] = 0.0; */
/* parts[k].a_legacy[1] = 0.0; */
/* parts[k].a_legacy[2] = 0.0; */
/* } */
for (k = 0; k < N; k++) {
parts[k].id = k;
parts[k].x[0] = ((double)rand()) / RAND_MAX;
parts[k].x[1] = ((double)rand()) / RAND_MAX;
parts[k].x[2] = ((double)rand()) / RAND_MAX;
parts[k].mass = 1.f;//((double)rand()) / RAND_MAX;
parts[k].a_legacy[0] = 0.0;
parts[k].a_legacy[1] = 0.0;
parts[k].a_legacy[2] = 0.0;
}
int ii, jj, kk;
for (ii = 0; ii < 8; ++ii) {
for (jj = 0; jj < 8; ++jj) {
for (kk = 0; kk < 8; ++kk) {
/* int ii, jj, kk; */
/* for (ii = 0; ii < 8; ++ii) { */
/* for (jj = 0; jj < 8; ++jj) { */
/* for (kk = 0; kk < 8; ++kk) { */
int index = 64*ii + 8*jj + kk;
parts[index].id = index;
parts[index].x[0] = ii*0.125 + 0.0625;
parts[index].x[1] = jj*0.125 + 0.0625;
parts[index].x[2] = kk*0.125 + 0.0625;
parts[index].mass = 1.;
parts[index].a_legacy[0] = 0.f;
parts[index].a_legacy[1] = 0.f;
parts[index].a_legacy[2] = 0.f;
}
}
}
/* int index = 64*ii + 8*jj + kk; */
/* parts[index].id = index; */
/* parts[index].x[0] = ii*0.125 + 0.0625; */
/* parts[index].x[1] = jj*0.125 + 0.0625; */
/* parts[index].x[2] = kk*0.125 + 0.0625; */
/* parts[index].mass = 1.; */
/* parts[index].a_legacy[0] = 0.f; */
/* parts[index].a_legacy[1] = 0.f; */
/* parts[index].a_legacy[2] = 0.f; */
/* } */
/* } */
/* } */
/* Otherwise, read them from a file. */
} else {
......@@ -1419,7 +1405,7 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
message("[check] root CoMx= %f %f", root->legacy.com[0], root->new.M_100);
message("[check] root CoMy= %f %f", root->legacy.com[1], root->new.M_010);
message("[check] root CoMz= %f %f", root->legacy.com[2], root->new.M_001);
#ifdef QUADRUPOLES
message("[check] | %f %f %f |", root->new.M_200, root->new.M_110, root->new.M_101);
message("[check] I = | %f %f %f |", root->new.M_110, root->new.M_020, root->new.M_011);
message("[check] | %f %f %f |", root->new.M_101, root->new.M_011, root->new.M_002);
......@@ -1445,6 +1431,7 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
message("[check] | %f %f %f |", M_200, M_110, M_101);
message("[check] I = | %f %f %f |", M_110, M_020, M_011);
message("[check] | %f %f %f |", M_101, M_011, M_002);
#endif
#if ICHECK >= 0
for (i = 0; i < N; ++i)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment