diff --git a/src/const.h b/src/const.h index d164f71b3e0e044276bb6bf253eedc8d435e736e..c290a3e73a524e9003cadb63f8bd45e8b3c51dac 100644 --- a/src/const.h +++ b/src/const.h @@ -59,6 +59,4 @@ #define SOURCETERMS_NONE //#define SOURCETERMS_SN_FEEDBACK -#define ICHECK 0 - #endif /* SWIFT_CONST_H */ diff --git a/src/drift.h b/src/drift.h index 9eb86b86973719eb78baf8209bbeed658ae72e7d..7e6efceecdbe468d362aad20ebecd2289be35703 100644 --- a/src/drift.h +++ b/src/drift.h @@ -39,7 +39,7 @@ * @param ti_current Integer end of time-step */ __attribute__((always_inline)) INLINE static void drift_gpart( - struct gpart *restrict gp, double dt, double timeBase, integertime_t ti_old, + struct gpart *restrict gp, float dt, double timeBase, integertime_t ti_old, integertime_t ti_current) { #ifdef SWIFT_DEBUG_CHECKS @@ -53,9 +53,6 @@ __attribute__((always_inline)) INLINE static void drift_gpart( gp->ti_drift = ti_current; #endif - // message("dt= %e", dt); - // fprintf(files_timestep[gp->id_or_neg_offset], "drift: dt=%e\n", dt); - /* Drift... */ gp->x[0] += gp->v_full[0] * dt; gp->x[1] += gp->v_full[1] * dt; @@ -78,7 +75,7 @@ __attribute__((always_inline)) INLINE static void drift_gpart( * @param ti_current Integer end of time-step */ __attribute__((always_inline)) INLINE static void drift_part( - struct part *restrict p, struct xpart *restrict xp, double dt, + struct part *restrict p, struct xpart *restrict xp, float dt, double timeBase, integertime_t ti_old, integertime_t ti_current) { #ifdef SWIFT_DEBUG_CHECKS diff --git a/src/kick.h b/src/kick.h index 61bdd194a2fbee47c959626abc97362fe5972250..a3759b024c1305a38b8814ccec54c44e99a52e2d 100644 --- a/src/kick.h +++ b/src/kick.h @@ -40,7 +40,7 @@ __attribute__((always_inline)) INLINE static void kick_gpart( double timeBase) { /* Time interval for this half-kick */ - const double dt = (ti_end - ti_start) * timeBase; + const float dt = (ti_end - ti_start) * timeBase; #ifdef SWIFT_DEBUG_CHECKS if (gp->ti_kick != ti_start) @@ -75,7 +75,7 @@ __attribute__((always_inline)) INLINE static void kick_part( integertime_t ti_end, double timeBase) { /* Time interval for this half-kick */ - const double dt = (ti_end - ti_start) * timeBase; + const float dt = (ti_end - ti_start) * timeBase; #ifdef SWIFT_DEBUG_CHECKS if (p->ti_kick != ti_start) diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h index 6ccee4d8a7bfcd3212144a4a659e1dccc3e26ab3..81b51ab2009ad599d0201708d78c8c64cac991dc 100644 --- a/src/potential/point_mass/potential.h +++ b/src/potential/point_mass/potential.h @@ -64,33 +64,27 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( const struct phys_const* restrict phys_const, const struct gpart* restrict g) { - const double G_newton = phys_const->const_newton_G; - const double dx = g->x[0] - potential->x; - const double dy = g->x[1] - potential->y; - const double dz = g->x[2] - potential->z; - /* const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz); */ - /* const float rinv2 = rinv * rinv; */ - /* const float rinv3 = rinv2 * rinv; */ - /* const float drdv = (g->x[0] - potential->x) * (g->v_full[0]) + */ - /* (g->x[1] - potential->y) * (g->v_full[1]) + */ - /* (g->x[2] - potential->z) * (g->v_full[2]); */ - /* const float dota_x = G_newton * potential->mass * rinv3 * */ - /* (-g->v_full[0] + 3.f * rinv2 * drdv * dx); */ - /* const float dota_y = G_newton * potential->mass * rinv3 * */ - /* (-g->v_full[1] + 3.f * rinv2 * drdv * dy); */ - /* const float dota_z = G_newton * potential->mass * rinv3 * */ - /* (-g->v_full[2] + 3.f * rinv2 * drdv * dz); */ - /* const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z; - */ - /* const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] - */ - /* + g->a_grav[2] * g->a_grav[2]; */ - - /* return potential->timestep_mult * sqrtf(a_2 / dota_2); */ - - const double r = sqrt(dx * dx + dy * dy + dz * dz); - return potential->timestep_mult * - sqrt(r * r * r / (G_newton * potential->mass)); + const float G_newton = phys_const->const_newton_G; + const float dx = g->x[0] - potential->x; + const float dy = g->x[1] - potential->y; + const float dz = g->x[2] - potential->z; + const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz); + const float rinv2 = rinv * rinv; + const float rinv3 = rinv2 * rinv; + const float drdv = (g->x[0] - potential->x) * (g->v_full[0]) + + (g->x[1] - potential->y) * (g->v_full[1]) + + (g->x[2] - potential->z) * (g->v_full[2]); + const float dota_x = G_newton * potential->mass * rinv3 * + (-g->v_full[0] + 3.f * rinv2 * drdv * dx); + const float dota_y = G_newton * potential->mass * rinv3 * + (-g->v_full[1] + 3.f * rinv2 * drdv * dy); + const float dota_z = G_newton * potential->mass * rinv3 * + (-g->v_full[2] + 3.f * rinv2 * drdv * dz); + const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z; + const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] + + g->a_grav[2] * g->a_grav[2]; + + return potential->timestep_mult * sqrtf(a_2 / dota_2); } /** @@ -111,11 +105,11 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration( double time, const struct external_potential* restrict potential, const struct phys_const* restrict phys_const, struct gpart* restrict g) { - const double dx = g->x[0] - potential->x; - const double dy = g->x[1] - potential->y; - const double dz = g->x[2] - potential->z; - const double rinv = 1. / sqrt(dx * dx + dy * dy + dz * dz); - const double rinv3 = rinv * rinv * rinv; + const float dx = g->x[0] - potential->x; + const float dy = g->x[1] - potential->y; + const float dz = g->x[2] - potential->z; + const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz); + const float rinv3 = rinv * rinv * rinv; g->a_grav[0] += -potential->mass * dx * rinv3; g->a_grav[1] += -potential->mass * dy * rinv3;