Commit d9d8742a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Correct P-M forces for multipole order < 3

parent 297ed246
......@@ -152,7 +152,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
/* In the case where the order is < 3, then there is only a monopole term left.
* We can default to the normal P-P interaction with the mass of the multipole
* and its CoM as the "particle" property */
#if 1 //SELF_GRAVITY_MULTIPOLE_ORDER < 3
#if 1 // SELF_GRAVITY_MULTIPOLE_ORDER < 3
float f_ij, pot_ij;
runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
......@@ -240,9 +240,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
m->M_000, r_s_inv, &f_ij, &pot_ij);
*f_x = -f_ij * r_x;
*f_y = -f_ij * r_y;
*f_z = -f_ij * r_z;
*f_x = f_ij * r_x;
*f_y = f_ij * r_y;
*f_z = f_ij * r_z;
*pot = pot_ij;
#else
......
......@@ -31,7 +31,7 @@
#include <float.h>
#include <math.h>
#define GADGET2_LONG_RANGE_CORRECTION
//#define GADGET2_LONG_RANGE_CORRECTION
#ifdef GADGET2_LONG_RANGE_CORRECTION
#define kernel_long_gravity_truncation_name "Gadget-like (using erfc())"
......@@ -121,7 +121,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_derivatives(
const float x = c1 * r;
/* e^(2r / r_s) */
const float exp_x = good_approx_expf(x);
const float exp_x = expf(x); // good_approx_expf(x);
/* 1 / alpha(w) */
const float a_inv = 1.f + exp_x;
......@@ -164,7 +164,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval(
#else
const float x = 2.f * u;
const float exp_x = good_approx_expf(x);
const float exp_x = expf(x); // good_approx_expf(x);
const float alpha = 1.f / (1.f + exp_x);
/* We want 2 - 2 exp(x) * alpha */
......@@ -197,7 +197,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
#else
const float x = 2.f * u;
const float exp_x = good_approx_expf(x);
const float exp_x = expf(x); // good_approx_expf(x);
const float alpha = 1.f / (1.f + exp_x);
/* We want 2*(x*alpha - x*alpha^2 - exp(x)*alpha + 1) */
......
......@@ -751,10 +751,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
ci->gparts, cj->gparts);
/* Then the M2P */
if (0)
runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
multi_j, dim, r_s_inv, e, ci->gparts,
gcount_i, cj);
runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
multi_j, dim, r_s_inv, e, ci->gparts,
gcount_i, cj);
}
if (cj_active && symmetric) {
......@@ -764,10 +763,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
cj->gparts, ci->gparts);
/* Then the M2P */
if (0)
runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
multi_i, dim, r_s_inv, e, cj->gparts,
gcount_j, ci);
runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
multi_i, dim, r_s_inv, e, cj->gparts,
gcount_j, ci);
}
} else {
......@@ -783,10 +781,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
ci->gparts, cj->gparts);
/* Then the M2P */
if (0)
runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
periodic, dim, e, ci->gparts, gcount_i,
cj);
runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
periodic, dim, e, ci->gparts, gcount_i, cj);
}
if (cj_active && symmetric) {
......@@ -796,10 +792,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
cj->gparts, ci->gparts);
/* Then the M2P */
if (0)
runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
periodic, dim, e, cj->gparts, gcount_j,
ci);
runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
periodic, dim, e, cj->gparts, gcount_j, ci);
}
}
}
......
......@@ -42,7 +42,7 @@ const double eps = 0.1;
*/
void check_value(double a, double b, const char *s) {
if (fabs(a - b) / fabs(a + b) > 2e-6 && fabs(a - b) > 1.e-6)
error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s);
error("Values are inconsistent: SWIFT:%12.15e true:%12.15e (%s)!", a, b, s);
}
/* Definitions of the potential and force that match
......@@ -92,25 +92,34 @@ int main(int argc, char *argv[]) {
#ifdef HAVE_FE_ENABLE_EXCEPT
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
/* Initialise a few things to get us going */
/* Non-truncated forces first */
double rlr = FLT_MAX;
struct engine e;
e.max_active_bin = num_time_bins;
e.time = 0.1f;
e.ti_current = 8;
e.time_base = 1e-10;
struct space s;
s.periodic = 0;
e.s = &s;
struct pm_mesh mesh;
mesh.periodic = 0;
mesh.dim[0] = 10.;
mesh.dim[1] = 10.;
mesh.dim[2] = 10.;
mesh.r_s_inv = 0.;
mesh.r_s = rlr;
mesh.r_s_inv = 1. / rlr;
mesh.r_cut_min = 0.;
mesh.r_cut_max = FLT_MAX;
e.mesh = &mesh;
struct gravity_props props;
props.a_smooth = 1.25;
props.theta_crit2 = 0.;
props.epsilon_cur = eps;
e.gravity_properties = &props;
......@@ -119,8 +128,6 @@ int main(int argc, char *argv[]) {
bzero(&r, sizeof(struct runner));
r.e = &e;
const double rlr = FLT_MAX;
/* Init the cache for gravity interaction */
gravity_cache_init(&r.ci_gravity_cache, num_tests);
gravity_cache_init(&r.cj_gravity_cache, num_tests);
......@@ -225,8 +232,8 @@ int main(int argc, char *argv[]) {
double acc_true =
acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr);
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);
/* 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); */
check_value(gp->potential, pot_true, "potential");
check_value(gp->a_grav[0], acc_true, "acceleration");
......@@ -262,13 +269,13 @@ int main(int argc, char *argv[]) {
const struct gravity_tensors *mpole = ci.multipole;
const double epsilon = gravity_get_softening(gp, &props);
double pot_true = potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0],
epsilon, rlr * FLT_MAX);
double acc_true = acceleration(
mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr * FLT_MAX);
double pot_true =
potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr);
double acc_true = acceleration(mpole->m_pole.M_000,
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);
/* 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");
......@@ -279,11 +286,47 @@ int main(int argc, char *argv[]) {
/* Reset the accelerations */
for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]);
/***************************************/
/* Test the high-order PM interactions */
/***************************************/
/***************************************/
/* Test the truncated PM interactions */
/***************************************/
rlr = 2.;
mesh.r_s = rlr;
mesh.r_s_inv = 1. / rlr;
mesh.periodic = 1;
s.periodic = 1;
props.epsilon_cur = FLT_MIN; /* No softening */
/* Now compute the forces */
runner_dopair_grav_pp(&r, &ci, &cj, 1);
/* Verify everything */
for (int n = 0; n < num_tests; ++n) {
const struct gpart *gp = &cj.gparts[n];
const struct gravity_tensors *mpole = ci.multipole;
const double epsilon = gravity_get_softening(gp, &props);
double pot_true =
potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr);
double acc_true = acceleration(mpole->m_pole.M_000,
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");
}
message("\n\t\t truncated P-M interactions all good\n");
/***************************************/
/* Test the high-order PM interactions */
/***************************************/
/* Reset the accelerations */
for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]);
#if SELF_GRAVITY_MULTIPOLE_ORDER >= 3
#if 0 // SELF_GRAVITY_MULTIPOLE_ORDER >= 3
/* Let's make ci more interesting */
free(ci.gparts);
......@@ -346,13 +389,10 @@ int main(int argc, char *argv[]) {
gp2->x[2] - gp->x[2]};
const double d = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
pot_true += potential(gp2->mass, d, epsilon, rlr * FLT_MAX);
acc_true[0] -=
acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[0] / d;
acc_true[1] -=
acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[1] / d;
acc_true[2] -=
acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[2] / d;
pot_true += potential(gp2->mass, d, epsilon, rlr);
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[2] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[2] / d;
}
message("x=%e f=%e f_true=%e pot=%e pot_true=%e %e %e",
......
......@@ -95,7 +95,7 @@ int main(int argc, char *argv[]) {
#ifdef HAVE_FE_ENABLE_EXCEPT
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
/* Initialise a few things to get us going */
struct engine e;
e.max_active_bin = num_time_bins;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment