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

Corrected mistake in calculation of multipoles.

parent d97468ea
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 29 additions and 26 deletions
examples/test_fmm.c
+ 29
− 26
View file @ e76ea9ba
......@@ -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)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment