diff --git a/examples/test_fmm.c b/examples/test_fmm.c index 5a0fb2cda47d17dfaf254b7fcac11ac56c4fc258..8c73cc8187fa7d55dd59ed14afb7d212d824dc5c 100644 --- a/examples/test_fmm.c +++ b/examples/test_fmm.c @@ -47,7 +47,7 @@ #define SANITY_CHECKS #define COM_AS_TASK #define COUNTERS -#define NO_QUADRUPOLES +#define QUADRUPOLES #include "multipoles.h" @@ -191,9 +191,9 @@ void comp_com(struct cell *c) { com[2] += cp->new.M_001 * cp_mass; mass += cp_mass; #ifdef QUADRUPOLES - M_200 += cp->new.M_200 + 0.5f * cp_mass * cp->new.M_100 * cp->new.M_100; - M_020 += cp->new.M_020 + 0.5f * cp_mass * cp->new.M_010 * cp->new.M_010; - M_002 += cp->new.M_002 + 0.5f * cp_mass * cp->new.M_001 * cp->new.M_001; + M_200 += 2. * cp->new.M_200 + cp_mass * cp->new.M_100 * cp->new.M_100; + M_020 += 2. * cp->new.M_020 + cp_mass * cp->new.M_010 * cp->new.M_010; + M_002 += 2. * cp->new.M_002 + cp_mass * cp->new.M_001 * cp->new.M_001; M_110 += cp->new.M_110 + cp_mass * cp->new.M_100 * cp->new.M_010; M_101 += cp->new.M_101 + cp_mass * cp->new.M_100 * cp->new.M_001; M_011 += cp->new.M_011 + cp_mass * cp->new.M_010 * cp->new.M_001; @@ -1420,28 +1420,31 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) { 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); - /* 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); */ - - /* double M_200 = 0., M_020 = 0., M_002 = 0., M_110 = 0., M_101 = 0., M_011 = 0.; */ - /* for (i = 0; i < N; ++i) { */ - /* M_200 += root->parts[i].x[0] * root->parts[i].x[0] * root->parts[i].mass; */ - /* M_020 += root->parts[i].x[1] * root->parts[i].x[1] * root->parts[i].mass; */ - /* M_002 += root->parts[i].x[2] * root->parts[i].x[2] * root->parts[i].mass; */ - /* M_110 += root->parts[i].x[0] * root->parts[i].x[1] * root->parts[i].mass; */ - /* M_101 += root->parts[i].x[0] * root->parts[i].x[2] * root->parts[i].mass; */ - /* M_011 += root->parts[i].x[1] * root->parts[i].x[2] * root->parts[i].mass; */ - /* } */ - /* M_200 = M_200 - root->new.M_000 * root->new.M_100 * root->new.M_100; */ - /* M_020 = M_020 - root->new.M_000 * root->new.M_010 * root->new.M_010; */ - /* M_002 = M_002 - root->new.M_000 * root->new.M_001 * root->new.M_001; */ - /* M_110 = M_110 - root->new.M_000 * root->new.M_100 * root->new.M_010; */ - /* M_101 = M_101 - root->new.M_000 * root->new.M_100 * root->new.M_001; */ - /* M_011 = M_011 - root->new.M_000 * root->new.M_010 * root->new.M_001; */ - /* 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); */ + 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); + + double M_200 = 0., M_020 = 0., M_002 = 0., M_110 = 0., M_101 = 0., M_011 = 0.; + for (i = 0; i < N; ++i) { + M_200 += root->parts[i].x[0] * root->parts[i].x[0] * root->parts[i].mass; + M_020 += root->parts[i].x[1] * root->parts[i].x[1] * root->parts[i].mass; + M_002 += root->parts[i].x[2] * root->parts[i].x[2] * root->parts[i].mass; + M_110 += root->parts[i].x[0] * root->parts[i].x[1] * root->parts[i].mass; + M_101 += root->parts[i].x[0] * root->parts[i].x[2] * root->parts[i].mass; + M_011 += root->parts[i].x[1] * root->parts[i].x[2] * root->parts[i].mass; + } + M_200 = M_200 - root->new.M_000 * root->new.M_100 * root->new.M_100; + M_020 = M_020 - root->new.M_000 * root->new.M_010 * root->new.M_010; + M_002 = M_002 - root->new.M_000 * root->new.M_001 * root->new.M_001; + M_110 = M_110 - root->new.M_000 * root->new.M_100 * root->new.M_010; + M_101 = M_101 - root->new.M_000 * root->new.M_100 * root->new.M_001; + M_011 = M_011 - root->new.M_000 * root->new.M_010 * root->new.M_001; + M_200 *= 0.5f; + M_020 *= 0.5f; + M_002 *= 0.5f; + 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); #if ICHECK >= 0 for (i = 0; i < N; ++i)