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

Only check the potential in the unit tests if we compile with the potential-enabled grvity scheme.

parent cc51ca15
No related branches found
No related tags found
No related merge requests found
...@@ -229,16 +229,18 @@ int main(int argc, char *argv[]) { ...@@ -229,16 +229,18 @@ int main(int argc, char *argv[]) {
const struct gpart *gp2 = &ci.gparts[0]; const struct gpart *gp2 = &ci.gparts[0];
const double epsilon = gravity_get_softening(gp, &props); const double epsilon = gravity_get_softening(gp, &props);
#if defined(POTENTIAL_GRAVITY)
double pot_true = double pot_true =
potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr); potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr);
check_value(gp->potential, pot_true, "potential");
#endif
double acc_true = double acc_true =
acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr); acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr);
check_value(gp->a_grav[0], acc_true, "acceleration");
/* message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - gp2->x[0], /* message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - gp2->x[0],
* gp->a_grav[0], acc_true, gp->potential, pot_true); */ * gp->a_grav[0], acc_true, gp->potential, pot_true); */
check_value(gp->potential, pot_true, "potential");
check_value(gp->a_grav[0], acc_true, "acceleration");
} }
message("\n\t\t P-P interactions all good\n"); message("\n\t\t P-P interactions all good\n");
...@@ -271,16 +273,18 @@ int main(int argc, char *argv[]) { ...@@ -271,16 +273,18 @@ int main(int argc, char *argv[]) {
const struct gravity_tensors *mpole = ci.multipole; const struct gravity_tensors *mpole = ci.multipole;
const double epsilon = gravity_get_softening(gp, &props); const double epsilon = gravity_get_softening(gp, &props);
#if defined(POTENTIAL_GRAVITY)
double pot_true = double pot_true =
potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr); potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr);
check_value(gp->potential, pot_true, "potential");
#endif
double acc_true = acceleration(mpole->m_pole.M_000, double acc_true = acceleration(mpole->m_pole.M_000,
gp->x[0] - mpole->CoM[0], epsilon, rlr); gp->x[0] - mpole->CoM[0], epsilon, rlr);
check_value(gp->a_grav[0], acc_true, "acceleration");
/* message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - /* message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] -
* mpole->CoM[0], gp->a_grav[0], acc_true, gp->potential, pot_true); */ * mpole->CoM[0], gp->a_grav[0], acc_true, gp->potential, pot_true); */
check_value(gp->potential, pot_true, "potential");
check_value(gp->a_grav[0], acc_true, "acceleration");
} }
message("\n\t\t basic P-M interactions all good\n"); message("\n\t\t basic P-M interactions all good\n");
...@@ -309,16 +313,18 @@ int main(int argc, char *argv[]) { ...@@ -309,16 +313,18 @@ int main(int argc, char *argv[]) {
const struct gravity_tensors *mpole = ci.multipole; const struct gravity_tensors *mpole = ci.multipole;
const double epsilon = gravity_get_softening(gp, &props); const double epsilon = gravity_get_softening(gp, &props);
#if defined(POTENTIAL_GRAVITY)
double pot_true = double pot_true =
potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr); potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr);
check_value(gp->potential, pot_true, "potential");
#endif
double acc_true = acceleration(mpole->m_pole.M_000, double acc_true = acceleration(mpole->m_pole.M_000,
gp->x[0] - mpole->CoM[0], epsilon, rlr); gp->x[0] - mpole->CoM[0], epsilon, rlr);
message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - mpole->CoM[0],
gp->a_grav[0], acc_true, gp->potential, pot_true);
check_value(gp->potential, pot_true, "potential");
check_value(gp->a_grav[0], acc_true, "acceleration"); check_value(gp->a_grav[0], acc_true, "acceleration");
/* message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - mpole->CoM[0], */
/* gp->a_grav[0], acc_true, gp->potential, pot_true); */
} }
message("\n\t\t truncated P-M interactions all good\n"); message("\n\t\t truncated P-M interactions all good\n");
...@@ -382,7 +388,10 @@ int main(int argc, char *argv[]) { ...@@ -382,7 +388,10 @@ int main(int argc, char *argv[]) {
for (int n = 0; n < num_tests; ++n) { for (int n = 0; n < num_tests; ++n) {
const struct gpart *gp = &cj.gparts[n]; const struct gpart *gp = &cj.gparts[n];
double pot_true = 0, acc_true[3] = {0., 0., 0.}; #if defined(POTENTIAL_GRAVITY)
double pot_true = 0;
#endif
double acc_true[3] = {0., 0., 0.};
for (int i = 0; i < 8; ++i) { for (int i = 0; i < 8; ++i) {
const struct gpart *gp2 = &ci.gparts[i]; const struct gpart *gp2 = &ci.gparts[i];
...@@ -392,20 +401,25 @@ int main(int argc, char *argv[]) { ...@@ -392,20 +401,25 @@ int main(int argc, char *argv[]) {
gp2->x[2] - gp->x[2]}; gp2->x[2] - gp->x[2]};
const double d = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]); const double d = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
#if defined(POTENTIAL_GRAVITY)
pot_true += potential(gp2->mass, d, epsilon, rlr); pot_true += potential(gp2->mass, d, epsilon, rlr);
#endif
acc_true[0] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[0] / d; acc_true[0] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[0] / d;
acc_true[1] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[1] / d; acc_true[1] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[1] / d;
acc_true[2] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[2] / d; acc_true[2] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[2] / d;
} }
#if defined(POTENTIAL_GRAVITY)
check_value_backend(gp->potential, pot_true, "potential", 1e-2, 1e-6);
#endif
check_value_backend(gp->a_grav[0], acc_true[0], "acceleration", 1e-2, 1e-6);
/* const struct gravity_tensors *mpole = ci.multipole; */ /* const struct gravity_tensors *mpole = ci.multipole; */
/* message("x=%e f=%e f_true=%e pot=%e pot_true=%e %e %e", */ /* message("x=%e f=%e f_true=%e pot=%e pot_true=%e %e %e", */
/* gp->x[0] - mpole->CoM[0], gp->a_grav[0], acc_true[0], /* gp->x[0] - mpole->CoM[0], gp->a_grav[0], acc_true[0],
* gp->potential, */ * gp->potential, */
/* pot_true, acc_true[1], acc_true[2]); */ /* pot_true, acc_true[1], acc_true[2]); */
check_value_backend(gp->potential, pot_true, "potential", 1e-2, 1e-6);
check_value_backend(gp->a_grav[0], acc_true[0], "acceleration", 1e-2, 1e-6);
} }
message("\n\t\t high-order P-M interactions all good\n"); message("\n\t\t high-order P-M interactions all good\n");
......
...@@ -184,13 +184,15 @@ int main(int argc, char *argv[]) { ...@@ -184,13 +184,15 @@ int main(int argc, char *argv[]) {
const double epsilon = gravity_get_softening(gp, &props); const double epsilon = gravity_get_softening(gp, &props);
#if defined(POTENTIAL_GRAVITY)
double pot_true = potential(c.gparts[0].mass, gp->x[0], epsilon, rlr); double pot_true = potential(c.gparts[0].mass, gp->x[0], epsilon, rlr);
check_value(gp->potential, pot_true, "potential");
#endif
double acc_true = acceleration(c.gparts[0].mass, gp->x[0], epsilon, rlr); double acc_true = acceleration(c.gparts[0].mass, gp->x[0], epsilon, rlr);
check_value(gp->a_grav[0], acc_true, "acceleration");
// message("x=%e f=%e f_true=%e", gp->x[0], gp->a_grav[0], acc_true); // message("x=%e f=%e f_true=%e", gp->x[0], gp->a_grav[0], acc_true);
check_value(gp->potential, pot_true, "potential");
check_value(gp->a_grav[0], acc_true, "acceleration");
} }
free(c.gparts); free(c.gparts);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment