From 96d88ca08b743e9b1adf035117cc295ddd2c32b0 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Sun, 9 Sep 2018 15:35:34 +0200 Subject: [PATCH] Added a debugging check to verify that the gparts have been initialised before receiving any force. --- src/gravity/Default/gravity.h | 5 +++ src/gravity/Default/gravity_part.h | 3 ++ src/gravity/Potential/gravity.h | 5 +++ src/gravity/Potential/gravity_part.h | 3 ++ src/mesh_gravity.c | 14 ++++++++- src/runner_doiact_grav.h | 46 +++++++++++++++++++++++++--- 6 files changed, 70 insertions(+), 6 deletions(-) diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h index 2713c9ee7a..83decb0cfb 100644 --- a/src/gravity/Default/gravity.h +++ b/src/gravity/Default/gravity.h @@ -155,6 +155,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart( #ifdef SWIFT_DEBUG_CHECKS gp->num_interacted = 0; + gp->initialised = 1; #endif } @@ -187,6 +188,10 @@ __attribute__((always_inline)) INLINE static void gravity_end_force( gp->a_grav_PM[1] *= const_G; gp->a_grav_PM[2] *= const_G; #endif + +#ifdef SWIFT_DEBUG_CHECKS + gp->initialised = 0; /* Ready for next step */ +#endif } /** diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h index bd73c56da8..f065e6d3a2 100644 --- a/src/gravity/Default/gravity_part.h +++ b/src/gravity/Default/gravity_part.h @@ -55,6 +55,9 @@ struct gpart { /* Time of the last kick */ integertime_t ti_kick; + /* Has this particle been initialised? */ + int initialised; + #endif #ifdef SWIFT_GRAVITY_FORCE_CHECKS diff --git a/src/gravity/Potential/gravity.h b/src/gravity/Potential/gravity.h index 3a6c0fba18..10628dcea9 100644 --- a/src/gravity/Potential/gravity.h +++ b/src/gravity/Potential/gravity.h @@ -151,6 +151,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart( #ifdef SWIFT_DEBUG_CHECKS gp->num_interacted = 0; + gp->initialised = 1; #endif } @@ -183,6 +184,10 @@ __attribute__((always_inline)) INLINE static void gravity_end_force( gp->a_grav_PM[1] *= const_G; gp->a_grav_PM[2] *= const_G; #endif + +#ifdef SWIFT_DEBUG_CHECKS + gp->initialised = 0; /* Ready for next step */ +#endif } /** diff --git a/src/gravity/Potential/gravity_part.h b/src/gravity/Potential/gravity_part.h index 252c18a4dc..229d801108 100644 --- a/src/gravity/Potential/gravity_part.h +++ b/src/gravity/Potential/gravity_part.h @@ -58,6 +58,9 @@ struct gpart { /* Time of the last kick */ integertime_t ti_kick; + /* Has this particle been initialised? */ + int initialised; + #endif #ifdef SWIFT_GRAVITY_FORCE_CHECKS diff --git a/src/mesh_gravity.c b/src/mesh_gravity.c index 7bf1dbf07e..d193f19082 100644 --- a/src/mesh_gravity.c +++ b/src/mesh_gravity.c @@ -467,8 +467,20 @@ void pm_mesh_interpolate_forces(const struct pm_mesh* mesh, for (int i = 0; i < gcount; ++i) { struct gpart* gp = &gparts[i]; - if (gpart_is_active(gp, e)) + if (gpart_is_active(gp, e)) { + +#ifdef SWIFT_DEBUG_CHECKS + /* Check that particles have been drifted to the current time */ + if (gp->ti_drift != e->ti_current) + error("gpart not drifted to current time"); + + /* Check that the particle was initialised */ + if (gp->initialised == 0) + error("Adding forces to an un-initialised gpart."); +#endif + mesh_to_gparts_CIC(gp, potential, N, cell_fac, dim); + } } #else error("No FFTW library found. Cannot compute periodic long-range forces."); diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index 34d3e690dd..961e64006f 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -119,6 +119,10 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c, error("gpart not drifted to current time"); if (c->multipole->pot.ti_init != e->ti_current) error("c->field tensor not initialised"); + + /* Check that the particle was initialised */ + if (gp->initialised == 0) + error("Adding forces to an un-initialised gpart."); #endif /* Apply the kernel */ gravity_L2P(pot, CoM, gp); @@ -223,6 +227,10 @@ static INLINE void runner_dopair_grav_pp_full( error("gpi not drifted to current time"); if (pjd < gcount_j && gparts_j[pjd].ti_drift != e->ti_current) error("gpj not drifted to current time"); + + /* Check that the particle was initialised */ + if (gparts_i[pid].initialised == 0) + error("Adding forces to an un-initialised gpart."); #endif /* Interact! */ @@ -349,6 +357,10 @@ static INLINE void runner_dopair_grav_pp_truncated( error("gpi not drifted to current time"); if (pjd < gcount_j && gparts_j[pjd].ti_drift != e->ti_current) error("gpj not drifted to current time"); + + /* Check that the particle was initialised */ + if (gparts_i[pid].initialised == 0) + error("Adding forces to an un-initialised gpart."); #endif /* Interact! */ @@ -433,6 +445,14 @@ static INLINE void runner_dopair_grav_pm_full( if (pid < gcount_i && !gpart_is_active(&gparts_i[pid], e)) error("Active particle went through the cache"); + /* Check that particles have been drifted to the current time */ + if (gparts_i[pid].ti_drift != e->ti_current) + error("gpi not drifted to current time"); + + /* Check that the particle was initialised */ + if (gparts_i[pid].initialised == 0) + error("Adding forces to an un-initialised gpart."); + if (pid >= gcount_i) error("Adding forces to padded particle"); #endif @@ -458,13 +478,13 @@ static INLINE void runner_dopair_grav_pm_full( const float r2 = dx * dx + dy * dy + dz * dz; -#ifdef SWIFT_DEBUG_CHECKSa +#ifdef SWIFT_DEBUG_CHECKS const float r_max_j = cj->multipole->r_max; const float r_max2 = r_max_j * r_max_j; const float theta_crit2 = e->gravity_properties->theta_crit2; - /* 1.01 to avoid FP rounding false-positives */ - if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.01, r2)) + /* Note: 1.1 to avoid FP rounding false-positives */ + if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2)) error( "use_mpole[i] set when M2P accept fails CoM=[%e %e %e] pos=[%e %e " "%e], rmax=%e", @@ -554,6 +574,14 @@ static INLINE void runner_dopair_grav_pm_truncated( if (pid < gcount_i && !gpart_is_active(&gparts_i[pid], e)) error("Active particle went through the cache"); + /* Check that particles have been drifted to the current time */ + if (gparts_i[pid].ti_drift != e->ti_current) + error("gpi not drifted to current time"); + + /* Check that the particle was initialised */ + if (gparts_i[pid].initialised == 0) + error("Adding forces to an un-initialised gpart."); + if (pid >= gcount_i) error("Adding forces to padded particle"); #endif @@ -582,8 +610,8 @@ static INLINE void runner_dopair_grav_pm_truncated( const float r_max2 = r_max_j * r_max_j; const float theta_crit2 = e->gravity_properties->theta_crit2; - /* 1.01 to avoid FP rounding false-positives */ - if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.01, r2)) + /* 1.1 to avoid FP rounding false-positives */ + if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2)) error( "use_mpole[i] set when M2P accept fails CoM=[%e %e %e] pos=[%e %e " "%e], rmax=%e", @@ -894,6 +922,10 @@ static INLINE void runner_doself_grav_pp_full( error("gpi not drifted to current time"); if (pjd < gcount && gparts[pjd].ti_drift != e->ti_current) error("gpj not drifted to current time"); + + /* Check that the particle was initialised */ + if (gparts[pid].initialised == 0) + error("Adding forces to an un-initialised gpart."); #endif /* Interact! */ @@ -1004,6 +1036,10 @@ static INLINE void runner_doself_grav_pp_truncated( error("gpi not drifted to current time"); if (pjd < gcount && gparts[pjd].ti_drift != e->ti_current) error("gpj not drifted to current time"); + + /* Check that the particle was initialised */ + if (gparts[pid].initialised == 0) + error("Adding forces to an un-initialised gpart."); #endif /* Interact! */ -- GitLab