From ba24398d38470261ff8bcb4a80b199cf2f4da467 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Thu, 1 Dec 2016 10:56:49 +0000
Subject: [PATCH] Replace ti_begin and ti_end by a time bin. Still need to fix
cooling(), do_recv() statistics() and propagate types everywhere.
---
examples/SedovBlast_3D/makeIC.py | 1 +
src/active.h | 27 +++++++++++------
src/engine.c | 4 +++
src/gravity/Default/gravity.h | 2 +-
src/gravity/Default/gravity_debug.h | 5 ++-
src/hydro/Gadget2/hydro.h | 6 ++--
src/hydro/Gadget2/hydro_debug.h | 4 +--
src/kick.h | 46 ++++++++++++++++++----------
src/runner.c | 47 +++++++++++++++++++++++------
src/statistics.c | 10 ++++--
src/timeline.h | 28 +++++++++++++++--
src/timestep.h | 31 ++++++++++---------
12 files changed, 147 insertions(+), 64 deletions(-)
diff --git a/examples/SedovBlast_3D/makeIC.py b/examples/SedovBlast_3D/makeIC.py
index e1b743c6cd..7d1a78188c 100644
--- a/examples/SedovBlast_3D/makeIC.py
+++ b/examples/SedovBlast_3D/makeIC.py
@@ -54,6 +54,7 @@ u[:] = P0 / (rho0 * (gamma - 1))
# Make the central particles detonate
index = argsort(r)
u[index[0:N_inject]] = E0 / (N_inject * m[0])
+print ids[index[0:N_inject]]
#--------------------------------------------------
diff --git a/src/active.h b/src/active.h
index cd0fb6cb21..2f23581342 100644
--- a/src/active.h
+++ b/src/active.h
@@ -101,14 +101,18 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active(
__attribute__((always_inline)) INLINE static int part_is_active(
const struct part *p, const struct engine *e) {
+ const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
#ifdef SWIFT_DEBUG_CHECKS
- if (p->ti_end < e->ti_current)
- error("particle in an impossible time-zone! p->ti_end=%d e->ti_current=%d",
- p->ti_end, e->ti_current);
+ if (ti_end < ti_current)
+ error(
+ "particle in an impossible time-zone! p->ti_end=%lld "
+ "e->ti_current=%lld",
+ ti_end, ti_current);
#endif
- // return (p->ti_end == e->ti_current);
- return (p->time_bin < get_max_active_bin(e->ti_current));
+ return (ti_end == ti_current);
}
/**
@@ -120,15 +124,18 @@ __attribute__((always_inline)) INLINE static int part_is_active(
__attribute__((always_inline)) INLINE static int gpart_is_active(
const struct gpart *gp, const struct engine *e) {
+ const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
+
#ifdef SWIFT_DEBUG_CHECKS
- if (gp->ti_end < e->ti_current)
+ if (ti_end < ti_current)
error(
- "g-particle in an impossible time-zone! gp->ti_end=%d e->ti_current=%d",
- gp->ti_end, e->ti_current);
+ "g-particle in an impossible time-zone! gp->ti_end=%lld "
+ "e->ti_current=%lld",
+ ti_end, ti_current);
#endif
- // return (gp->ti_end == e->ti_current);
- return (gp->time_bin < get_max_active_bin(e->ti_current));
+ return (ti_end == ti_current);
}
#endif /* SWIFT_ACTIVE_H */
diff --git a/src/engine.c b/src/engine.c
index e989aefd53..af6e222d11 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2643,6 +2643,10 @@ void engine_step(struct engine *e) {
engine_launch(e, e->nr_threads);
TIMER_TOC(timer_runners);
+ for (size_t i = 0; i < e->s->nr_parts; ++i) {
+ if (e->s->parts[i].time_bin == 0) error("Particle in bin 0");
+ }
+
/* Save some statistics */
if (e->time - e->timeLastStatistics >= e->deltaTimeStatistics) {
engine_print_stats(e);
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 4175040dc1..89c7ac66d4 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -53,7 +53,7 @@ gravity_compute_timestep_self(const struct gpart* const gp) {
__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
struct gpart* gp) {
- gp->time_bin = -1;
+ gp->time_bin = 0;
gp->epsilon = 0.; // MATTHIEU
}
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index c284f543b3..4c3c41c943 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -24,10 +24,9 @@ __attribute__((always_inline)) INLINE static void gravity_debug_particle(
printf(
"x=[%.3e,%.3e,%.3e], "
"v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
- "mass=%.3e t_begin=%d, t_end=%d\n",
+ "mass=%.3e time_bin=%d\n",
p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
- p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->ti_begin,
- p->ti_end);
+ p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->time_bin);
}
#endif /* SWIFT_DEFAULT_GRAVITY_DEBUG_H */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 7a7d196ca3..e80ea9f0c0 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -212,7 +212,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part *restrict p, struct xpart *restrict xp) {
- p->time_bin = -1;
+ p->time_bin = 0;
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
@@ -385,8 +385,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->rho *= expf(w2);
/* Drift the pressure */
- const integertime_t ti_begin = get_integer_time_begin(t1, p->time_bin);
- const integertime_t ti_end = get_integer_time_end(t1, p->time_bin);
+ const integertime_t ti_begin = get_integer_time_begin(t0, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(t0, p->time_bin);
const float dt_entr = (t1 - (ti_begin + ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, dt_entr);
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 656299b383..a9630a128d 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -27,14 +27,14 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"h=%.3e, wcount=%.3f, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, "
"P=%.3e, P_over_rho2=%.3e, S=%.3e, dS/dt=%.3e, c=%.3e\n"
"divV=%.3e, rotV=[%.3e,%.3e,%.3e], balsara=%.3e \n "
- "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
+ "v_sig=%e dh/dt=%.3e time_bin=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
p->rho, hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy,
p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
- p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end);
+ p->force.v_sig, p->force.h_dt, p->time_bin);
}
#endif /* SWIFT_GADGET2_HYDRO_DEBUG_H */
diff --git a/src/kick.h b/src/kick.h
index 91177f91b7..c2515e4f4d 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -32,20 +32,27 @@
*
* @param gp The #gpart to kick.
* @param new_dti The (integer) time-step for this kick.
+ * @param ti_current The current (integer) time.
* @param timeBase The minimal allowed time-step size.
*/
__attribute__((always_inline)) INLINE static void kick_gpart(
- struct gpart *restrict gp, int new_dti, double timeBase) {
+ struct gpart *restrict gp, integertime_t new_dti, integertime_t ti_current,
+ double timeBase) {
/* Compute the time step for this kick */
- const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
- const int ti_end = gp->ti_end + new_dti / 2;
+ const integertime_t old_ti_begin =
+ get_integer_time_begin(ti_current, gp->time_bin);
+ const integertime_t old_ti_end =
+ get_integer_time_end(ti_current, gp->time_bin);
+ const int ti_start = (old_ti_begin + old_ti_end) / 2;
+ const int ti_end = old_ti_end + new_dti / 2;
const float dt = (ti_end - ti_start) * timeBase;
- const float half_dt = (ti_end - gp->ti_end) * timeBase;
+ const float half_dt = (ti_end - old_ti_end) * timeBase;
/* Move particle forward in time */
- gp->ti_begin = gp->ti_end;
- gp->ti_end = gp->ti_begin + new_dti;
+ // gp->ti_begin = gp->ti_end;
+ // gp->ti_end = gp->ti_begin + new_dti;
+ gp->time_bin = get_time_bin(new_dti);
/* Kick particles in momentum space */
gp->v_full[0] += gp->a_grav[0] * dt;
@@ -62,26 +69,33 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
* @param p The #part to kick.
* @param xp The #xpart of the particle.
* @param new_dti The (integer) time-step for this kick.
+ * @param ti_current The current (integer) time.
* @param timeBase The minimal allowed time-step size.
*/
__attribute__((always_inline)) INLINE static void kick_part(
struct part *restrict p, struct xpart *restrict xp, int new_dti,
- double timeBase) {
+ integertime_t ti_current, double timeBase) {
/* Compute the time step for this kick */
- const integertime_t ti_begin = get_integer_time_begin(t1, p->time_bin);
- const integertime_t ti_end = get_integer_time_end(t1, p->time_bin);
- const int ti_start = (ti_begin + ti_end) / 2;
- const int ti_end = ti_end + new_dti / 2;
+ const integertime_t old_ti_begin =
+ get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t old_ti_end =
+ get_integer_time_end(ti_current, p->time_bin);
+ const int ti_start = (old_ti_begin + old_ti_end) / 2;
+ const int ti_end = old_ti_end + new_dti / 2;
const float dt = (ti_end - ti_start) * timeBase;
- const float half_dt = (ti_end - p->ti_end) * timeBase;
+ const float half_dt = (ti_end - old_ti_end) * timeBase;
/* Move particle forward in time */
- p->ti_begin = p->ti_end;
- p->ti_end = p->ti_begin + new_dti;
+ // p->ti_begin = p->ti_end;
+ // p->ti_end = p->ti_begin + new_dti;
+ p->time_bin = get_time_bin(new_dti);
+ // if (p->id == 116650) message("Time bin=%d new_dti=%d", p->time_bin,
+ // new_dti);
if (p->gpart != NULL) {
- p->gpart->ti_begin = p->ti_begin;
- p->gpart->ti_end = p->ti_end;
+ // p->gpart->ti_begin = p->ti_begin;
+ // p->gpart->ti_end = p->ti_end;
+ p->gpart->time_bin = get_time_bin(new_dti);
}
/* Get the acceleration */
diff --git a/src/runner.c b/src/runner.c
index 23e8aa0f69..a4da82a234 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -113,6 +113,9 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
#include "runner_doiact_fft.h"
#include "runner_doiact_grav.h"
+#undef ICHECK
+#define ICHECK -116650
+
/**
* @brief Perform source terms
*
@@ -205,7 +208,7 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
const int count = c->count;
- const int ti_current = r->e->ti_current;
+ // const int ti_current = r->e->ti_current;
const struct cooling_function_data *cooling_func = r->e->cooling_func;
const struct phys_const *constants = r->e->physical_constants;
const struct UnitSystem *us = r->e->internalUnits;
@@ -227,9 +230,11 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xp = &xparts[i];
/* Kick has already updated ti_end, so need to check ti_begin */
- if (p->ti_begin == ti_current) {
+ // if (p->ti_begin == ti_current) {
+ {
- const double dt = (p->ti_end - p->ti_begin) * timeBase;
+ // const double dt = (p->ti_end - p->ti_begin) * timeBase;
+ const double dt = 1. * timeBase; // MATTHIEU
cooling_cool_part(constants, us, cooling_func, p, xp, dt);
}
@@ -906,6 +911,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const double timeBase = e->timeBase;
+ const int ti_current = e->ti_current;
const int count = c->count;
const int gcount = c->gcount;
struct part *restrict parts = c->parts;
@@ -946,15 +952,17 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
const int new_dti = get_gpart_timestep(gp, e);
/* Now we have a time step, proceed with the kick */
- kick_gpart(gp, new_dti, timeBase);
+ kick_gpart(gp, new_dti, e->ti_current, timeBase);
/* Number of updated g-particles */
g_updated++;
}
/* Minimal time for next end of time-step */
- ti_end_min = min(gp->ti_end, ti_end_min);
- ti_end_max = max(gp->ti_end, ti_end_max);
+ const integertime_t ti_end =
+ get_integer_time_end(ti_current, gp->time_bin);
+ ti_end_min = min((int)ti_end, ti_end_min);
+ ti_end_max = max((int)ti_end, ti_end_max);
}
}
@@ -967,9 +975,18 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
struct part *restrict p = &parts[k];
struct xpart *restrict xp = &xparts[k];
+ /* Minimal time for next end of time-step */
+ integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ if (p->id == ICHECK)
+ message("Particle in kick ti_end=%lld ti_current=%d", ti_end,
+ e->ti_current);
+
/* If particle needs to be kicked */
if (part_is_active(p, e)) {
+ if (p->id == ICHECK) message("Particle active in kick");
+
/* First, finish the force loop */
hydro_end_force(p);
if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
@@ -977,17 +994,27 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
/* Compute the next timestep (hydro condition) */
const int new_dti = get_part_timestep(p, xp, e);
+ if (p->id == ICHECK)
+ message("time_step=%d (%e)", new_dti, new_dti * e->timeBase);
+
/* Now we have a time step, proceed with the kick */
- kick_part(p, xp, new_dti, timeBase);
+ kick_part(p, xp, new_dti, e->ti_current, timeBase);
+
+ // if (p->id == ICHECK) printParticle_single(p, xp);
/* Number of updated particles */
updated++;
if (p->gpart != NULL) g_updated++;
+
+ ti_end += get_integer_timestep(p->time_bin);
}
- /* Minimal time for next end of time-step */
- ti_end_min = min(p->ti_end, ti_end_min);
- ti_end_max = max(p->ti_end, ti_end_max);
+ if (p->id == ICHECK)
+ message("ti_current = %d dti=%lld ti_end=%lld (%f)", ti_current,
+ get_integer_timestep(p->time_bin), ti_end,
+ ti_end * e->timeBase);
+ ti_end_min = min((int)ti_end, ti_end_min);
+ ti_end_max = max((int)ti_end, ti_end_max);
}
}
diff --git a/src/statistics.c b/src/statistics.c
index 7a567a447a..1130a790f4 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -125,7 +125,10 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
const struct gpart *gp = (p->gpart != NULL) ? gp = p->gpart : NULL;
/* Get useful variables */
- const float dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+ const float dt = (ti_current - (ti_begin + ti_end) / 2) * timeBase;
const double x[3] = {p->x[0], p->x[1], p->x[2]};
float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
if (gp != NULL) {
@@ -206,7 +209,10 @@ void stats_collect_gpart_mapper(void *map_data, int nr_gparts,
if (gp->id_or_neg_offset < 0) continue;
/* Get useful variables */
- const float dt = (ti_current - (gp->ti_begin + gp->ti_end) / 2) * timeBase;
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current, gp->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
+ const float dt = (ti_current - (ti_begin + ti_end) / 2) * timeBase;
const double x[3] = {gp->x[0], gp->x[1], gp->x[2]};
const float v[3] = {gp->v_full[0] + gp->a_grav[0] * dt,
gp->v_full[1] + gp->a_grav[1] * dt,
diff --git a/src/timeline.h b/src/timeline.h
index 2d217a3822..23d2e55e3c 100644
--- a/src/timeline.h
+++ b/src/timeline.h
@@ -24,8 +24,11 @@
/* Local headers. */
#include "inline.h"
+#include "intrinsics.h"
-typedef long long integertime_t;
+#include <math.h>
+
+typedef unsigned long long integertime_t;
typedef char timebin_t;
/*! The number of time bins */
@@ -34,6 +37,8 @@ typedef char timebin_t;
/*! The maximal number of timesteps in a simulation */
#define max_nr_timesteps (1LL << (num_time_bins + 1))
+//#define ceil_div(x ,y) (x + y - 1) / y
+
/**
* @brief Returns the integer time interval corresponding to a time bin
*
@@ -45,6 +50,17 @@ static INLINE integertime_t get_integer_timestep(timebin_t bin) {
return 1LL << (bin + 1);
}
+/**
+ * @brief Returns the time bin corresponding to a given time_step size.
+ *
+ * Assumes that integertime_t maps to an unsigned long long.
+ */
+static INLINE timebin_t get_time_bin(integertime_t time_step) {
+
+ /* ((int) log_2(time_step)) - 1 */
+ return 62 - intrinsics_clzll(time_step);
+}
+
/**
* @brief Returns the physical time interval corresponding to a time bin.
*
@@ -67,7 +83,10 @@ static INLINE integertime_t get_integer_time_begin(integertime_t ti_current,
timebin_t bin) {
const integertime_t dti = get_integer_timestep(bin);
- return ti_current % dti;
+ if (dti == 0)
+ return 0;
+ else
+ return dti * ((ti_current - 1) / dti);
}
/**
@@ -81,7 +100,10 @@ static INLINE integertime_t get_integer_time_end(integertime_t ti_current,
timebin_t bin) {
const integertime_t dti = get_integer_timestep(bin);
- return ti_current % dti + dti;
+ if (dti == 0)
+ return 0;
+ else
+ return dti * ceil((double)ti_current / (double)dti);
}
/**
diff --git a/src/timestep.h b/src/timestep.h
index db52911ec1..1ac934d3a3 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -23,31 +23,34 @@
#include "../config.h"
/* Local headers. */
-#include "const.h"
#include "cooling.h"
#include "debug.h"
+#include "timeline.h"
+
/**
* @brief Compute a valid integer time-step form a given time-step
*
* @param new_dt The time-step to convert.
- * @param ti_begin The (integer) start of the previous time-step.
- * @param ti_end The (integer) end of the previous time-step.
+ * @param old_bin The old time bin.
+ * @param ti_current The current time on the integer time-line.
* @param timeBase_inv The inverse of the system's minimal time-step.
*/
-__attribute__((always_inline)) INLINE static int get_integer_timestep(
- float new_dt, int ti_begin, int ti_end, double timeBase_inv) {
+__attribute__((always_inline)) INLINE static int make_integer_timestep(
+ float new_dt, timebin_t old_bin, integertime_t ti_current,
+ double timeBase_inv) {
/* Convert to integer time */
- int new_dti = (int)(new_dt * timeBase_inv);
+ integertime_t new_dti = (integertime_t)(new_dt * timeBase_inv);
- /* Recover the current timestep */
- const int current_dti = ti_end - ti_begin;
+ /* Current time-step */
+ integertime_t current_dti = get_integer_timestep(old_bin);
+ integertime_t ti_end = get_integer_time_end(ti_current, old_bin);
/* Limit timestep increase */
- if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
+ if (old_bin > 0) new_dti = min(new_dti, 2 * current_dti);
/* Put this timestep on the time line */
- int dti_timeline = max_nr_timesteps;
+ integertime_t dti_timeline = max_nr_timesteps;
while (new_dti < dti_timeline) dti_timeline /= 2;
new_dti = dti_timeline;
@@ -82,8 +85,8 @@ __attribute__((always_inline)) INLINE static int get_gpart_timestep(
new_dt = max(new_dt, e->dt_min);
/* Convert to integer time */
- const int new_dti =
- get_integer_timestep(new_dt, gp->ti_begin, gp->ti_end, e->timeBase_inv);
+ const int new_dti = make_integer_timestep(new_dt, gp->time_bin, e->ti_current,
+ e->timeBase_inv);
return new_dti;
}
@@ -138,8 +141,8 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
new_dt = max(new_dt, e->dt_min);
/* Convert to integer time */
- const int new_dti =
- get_integer_timestep(new_dt, p->ti_begin, p->ti_end, e->timeBase_inv);
+ const int new_dti = make_integer_timestep(new_dt, p->time_bin, e->ti_current,
+ e->timeBase_inv);
return new_dti;
}
--
GitLab