From f4fef206aa91070a0a3a185b65f0b79696f71115 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <matthieu.schaller@durham.ac.uk> Date: Sat, 20 Feb 2016 18:09:28 +0000 Subject: [PATCH] Propagated the changes everywhere --- src/cell.c | 15 ++--- src/cell.h | 5 +- src/debug.c | 8 +-- src/engine.c | 74 ++++++++++++----------- src/engine.h | 26 +++++--- src/gravity/Default/gravity_part.h | 4 +- src/hydro/Gadget2/hydro.h | 12 ++-- src/hydro/Gadget2/hydro_debug.h | 4 +- src/hydro/Gadget2/hydro_part.h | 4 +- src/runner.c | 97 +++++++++++++++++------------- src/runner_doiact.h | 88 +++++++++++++-------------- src/runner_doiact_grav.h | 14 ++--- src/space.c | 22 ++++--- 13 files changed, 202 insertions(+), 171 deletions(-) diff --git a/src/cell.c b/src/cell.c index 6d84dafdaa..032efcac4d 100644 --- a/src/cell.c +++ b/src/cell.c @@ -88,8 +88,8 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) { /* Unpack the current pcell. */ c->h_max = pc->h_max; - c->t_end_min = pc->t_end_min; - c->t_end_max = pc->t_end_max; + c->ti_end_min = pc->ti_end_min; + c->ti_end_max = pc->ti_end_max; c->count = pc->count; c->tag = pc->tag; @@ -162,8 +162,8 @@ int cell_pack(struct cell *c, struct pcell *pc) { /* Start by packing the data of the current cell. */ pc->h_max = c->h_max; - pc->t_end_min = c->t_end_min; - pc->t_end_max = c->t_end_max; + pc->ti_end_min = c->ti_end_min; + pc->ti_end_max = c->ti_end_max; pc->count = c->count; c->tag = pc->tag = atomic_inc(&cell_next_tag) % cell_max_tag; @@ -559,8 +559,8 @@ void cell_init_parts(struct cell *c, void *data) { struct xpart *xp = c->xparts; for (int i = 0; i < c->count; ++i) { - p[i].t_begin = 0.; - p[i].t_end = 0.; + p[i].ti_begin = 0; + p[i].ti_end = 0; xp[i].v_full[0] = p[i].v[0]; xp[i].v_full[1] = p[i].v[1]; xp[i].v_full[2] = p[i].v[2]; @@ -568,7 +568,8 @@ void cell_init_parts(struct cell *c, void *data) { hydro_init_part(&p[i]); hydro_reset_acceleration(&p[i]); } - c->t_end_min = 0.; + c->ti_end_min = 0; + c->ti_end_max = 0; } /** diff --git a/src/cell.h b/src/cell.h index 0ba63d1f03..b0451b311f 100644 --- a/src/cell.h +++ b/src/cell.h @@ -40,7 +40,8 @@ extern int cell_next_tag; struct pcell { /* Stats on this cell's particles. */ - double h_max, t_end_min, t_end_max; + double h_max; + int ti_end_min, ti_end_max; /* Number of particles in this cell. */ int count; @@ -65,7 +66,7 @@ struct cell { double h_max; /* Minimum and maximum end of time step in this cell. */ - double t_end_min, t_end_max; + int ti_end_min, ti_end_max; /* Minimum dimension, i.e. smallest edge of this cell. */ float dmin; diff --git a/src/debug.c b/src/debug.c index 2be981175e..a9d702b4aa 100644 --- a/src/debug.c +++ b/src/debug.c @@ -75,12 +75,12 @@ void printgParticle(struct gpart *parts, long long int id, int N) { if (parts[i].id == -id || (parts[i].id > 0 && parts[i].part->id == id)) { printf( "## gParticle[%d]: id=%lld, x=[%.16e,%.16e,%.16e], " - "v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], m=%.3e, t_begin=%.3e, " - "t_end=%.3e\n", + "v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], m=%.3e, t_begin=%d, " + "t_end=%d\n", i, parts[i].part->id, parts[i].x[0], parts[i].x[1], parts[i].x[2], parts[i].v[0], parts[i].v[1], parts[i].v[2], parts[i].a[0], - parts[i].a[1], parts[i].a[2], parts[i].mass, parts[i].t_begin, - parts[i].t_end); + parts[i].a[1], parts[i].a[2], parts[i].mass, parts[i].ti_begin, + parts[i].ti_end); found = 1; } diff --git a/src/engine.c b/src/engine.c index 2cca53610c..f72a4fd7df 100644 --- a/src/engine.c +++ b/src/engine.c @@ -52,6 +52,7 @@ #include "cycle.h" #include "debug.h" #include "error.h" +#include "minmax.h" #include "part.h" #include "timers.h" @@ -1322,7 +1323,7 @@ int engine_marktasks(struct engine *e) { struct scheduler *s = &e->sched; int k, nr_tasks = s->nr_tasks, *ind = s->tasks_ind; struct task *t, *tasks = s->tasks; - float t_end = e->time; + float ti_end = e->ti_current; struct cell *ci, *cj; // ticks tic = getticks(); @@ -1384,7 +1385,7 @@ int engine_marktasks(struct engine *e) { (t->type == task_type_sub && t->cj == NULL)) { /* Set this task's skip. */ - t->skip = (t->ci->t_end_min > t_end); + t->skip = (t->ci->ti_end_min > ti_end); } /* Pair? */ @@ -1396,7 +1397,7 @@ int engine_marktasks(struct engine *e) { cj = t->cj; /* Set this task's skip. */ - t->skip = (ci->t_end_min > t_end && cj->t_end_min > t_end); + t->skip = (ci->ti_end_min > ti_end && cj->ti_end_min > ti_end); /* Too much particle movement? */ if (t->tight && @@ -1421,7 +1422,7 @@ int engine_marktasks(struct engine *e) { /* Kick? */ else if (t->type == task_type_kick) { - t->skip = (t->ci->t_end_min > t_end); + t->skip = (t->ci->ti_end_min > ti_end); } /* Drift? */ @@ -1431,7 +1432,7 @@ int engine_marktasks(struct engine *e) { /* Init? */ else if (t->type == task_type_init) { /* Set this task's skip. */ - t->skip = (t->ci->t_end_min > t_end); + t->skip = (t->ci->ti_end_min > ti_end); } /* None? */ @@ -1617,7 +1618,7 @@ void engine_barrier(struct engine *e, int tid) { void engine_collect_kick(struct cell *c) { int updated = 0; - float t_end_min = FLT_MAX, t_end_max = 0.0f; + int ti_end_min = max_nr_timesteps, ti_end_max = 0; double e_kin = 0.0, e_int = 0.0, e_pot = 0.0; float mom[3] = {0.0f, 0.0f, 0.0f}, ang[3] = {0.0f, 0.0f, 0.0f}; struct cell *cp; @@ -1639,8 +1640,8 @@ void engine_collect_kick(struct cell *c) { engine_collect_kick(cp); /* And update */ - t_end_min = fminf(t_end_min, cp->t_end_min); - t_end_max = fmaxf(t_end_max, cp->t_end_max); + ti_end_min = min(ti_end_min, cp->ti_end_min); + ti_end_max = max(ti_end_max, cp->ti_end_max); updated += cp->updated; e_kin += cp->e_kin; e_int += cp->e_int; @@ -1655,8 +1656,8 @@ void engine_collect_kick(struct cell *c) { } /* Store the collected values in the cell. */ - c->t_end_min = t_end_min; - c->t_end_max = t_end_max; + c->ti_end_min = ti_end_min; + c->ti_end_max = ti_end_max; c->updated = updated; c->e_kin = e_kin; c->e_int = e_int; @@ -1797,7 +1798,7 @@ void engine_step(struct engine *e) { int k; int updates = 0; - float t_end_min = FLT_MAX, t_end_max = 0.f; + int ti_end_min = max_nr_timesteps, ti_end_max = 0; double e_pot = 0.0, e_int = 0.0, e_kin = 0.0; float mom[3] = {0.0, 0.0, 0.0}; float ang[3] = {0.0, 0.0, 0.0}; @@ -1815,8 +1816,8 @@ void engine_step(struct engine *e) { engine_collect_kick(c); /* And aggregate */ - t_end_min = fminf(t_end_min, c->t_end_min); - t_end_max = fmaxf(t_end_max, c->t_end_max); + ti_end_min = min(ti_end_min, c->ti_end_min); + ti_end_max = max(ti_end_max, c->ti_end_max); e_kin += c->e_kin; e_int += c->e_int; e_pot += c->e_pot; @@ -1831,37 +1832,40 @@ void engine_step(struct engine *e) { /* Aggregate the data from the different nodes. */ #ifdef WITH_MPI - double in[4], out[4]; - out[0] = t_end_min; - if (MPI_Allreduce(out, in, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD) != + int in_i[4], out_i[4]; + out_t[0] = ti_end_min; + if (MPI_Allreduce(out_i, in_i, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD) != MPI_SUCCESS) error("Failed to aggregate t_end_min."); - t_end_min = in[0]; - out[0] = t_end_max; - if (MPI_Allreduce(out, in, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD) != + ti_end_min = in_i[0]; + out_i[0] = ti_end_max; + if (MPI_Allreduce(out_i, in_i, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD) != MPI_SUCCESS) error("Failed to aggregate t_end_max."); - t_end_max = in[0]; - out[0] = updates; - out[1] = e_kin; - out[2] = e_int; - out[3] = e_pot; - if (MPI_Allreduce(out, in, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) != + ti_end_max = in_i[0]; + doubles in_d[4], out_d[4]; + out_d[0] = updates; + out_d[1] = e_kin; + out_d[2] = e_int; + out_d[3] = e_pot; + if (MPI_Allreduce(out_d, in_d, 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS) error("Failed to aggregate energies."); - updates = in[0]; - e_kin = in[1]; - e_int = in[2]; - e_pot = in[3]; + updates = in_d[0]; + e_kin = in_d[1]; + e_int = in_d[2]; + e_pot = in_d[3]; #endif // message("\nDRIFT\n"); /* Move forward in time */ - e->timeOld = e->time; - e->time = t_end_min; + e->ti_old = e->ti_current; + e->ti_current = ti_end_min; e->step += 1; - e->timeStep = e->time - e->timeOld; + e->time = e->ti_current * e->timeBase; + e->timeOld = e->ti_old * e->timeBase; + e->timeStep = (e->ti_current - e->ti_old) * e->timeBase; /* Drift everybody */ engine_launch(e, e->nr_threads, 1 << task_type_drift, 0); @@ -2216,6 +2220,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads, e->step = 0; e->nullstep = 0; e->time = 0.0; + e->ti_current = 0; e->nr_nodes = nr_nodes; e->nodeID = nodeID; e->proxy_ind = NULL; @@ -2224,8 +2229,6 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads, e->forcerepart = 0; e->links = NULL; e->nr_links = 0; - e->timeBegin = timeBegin; - e->timeEnd = timeEnd; e->timeOld = timeBegin; e->time = timeBegin; e->timeStep = 0.; @@ -2261,11 +2264,12 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads, /* Print policy */ engine_print_policy(e); - + /* Deal with timestep */ if ((e->policy & engine_policy_fixdt) == engine_policy_fixdt) { e->dt_min = e->dt_max; } + e->timeBase = (timeEnd - timeBegin) / max_nr_timesteps; /* Construct types for MPI communications */ #ifdef WITH_MPI diff --git a/src/engine.h b/src/engine.h index 46c89a331f..b5562e3d34 100644 --- a/src/engine.h +++ b/src/engine.h @@ -67,6 +67,9 @@ extern const char *engine_policy_names[]; /* The rank of the engine as a global variable (for messages). */ extern int engine_rank; +/* The maximal number of timesteps in a simulation */ +#define max_nr_timesteps (1 << 28) + /* Mini struct to link cells to density/force tasks. */ struct link { @@ -99,20 +102,27 @@ struct engine { float dt_min, dt_max; /* Time of the simulation beginning */ - float timeBegin; - + //float timeBegin; + //int ti_begin; + /* Time of the simulation end */ - float timeEnd; - + //float timeEnd; + //int ti_end; + /* The previous system time. */ - float timeOld; - + double timeOld; + int ti_old; + /* The current system time. */ - float time; + double time; + int ti_current; /* Time step */ - float timeStep; + double timeStep; + /* Time base */ + double timeBase; + /* File for statistics */ FILE *file_stats; diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h index bad5c7318e..17262f9d41 100644 --- a/src/gravity/Default/gravity_part.h +++ b/src/gravity/Default/gravity_part.h @@ -35,10 +35,10 @@ struct gpart { float mass; /* Particle time of beginning of time-step. */ - float t_begin; + int ti_begin; /* Particle time of end of time-step. */ - float t_end; + int ti_end; /* Anonymous union for id/part. */ union { diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index 5194bccc6c..b135bed20c 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -124,10 +124,11 @@ __attribute__((always_inline)) * * @param p The particle to act upon * @param xp The extended particle data to act upon - * @param time The current time + * @param ti_current The current time (on the timeline) + * @param timeBase The minimal time-step size */ __attribute__((always_inline)) INLINE static void hydro_prepare_force( - struct part* p, struct xpart* xp, float time) { + struct part* p, struct xpart* xp, int ti_current, double timeBase) { /* Compute the norm of the curl */ p->force.curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] + @@ -135,7 +136,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( p->density.rot_v[2] * p->density.rot_v[2]); /* Compute the pressure */ - const float dt = time - 0.5f * (p->t_begin + p->t_end); + const float dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase; p->force.pressure = (p->entropy + p->entropy_dt * dt) * powf(p->rho, const_hydro_gamma); @@ -175,12 +176,13 @@ __attribute__((always_inline)) * @param xp The extended data of the particle * @param t0 The time at the start of the drift * @param t1 The time at the end of the drift + * @param timeBase The minimal time-step size */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( - struct part* p, struct xpart* xp, float t0, float t1) { + struct part* p, struct xpart* xp, int t0, int t1, double timeBase) { /* Drift the pressure */ - const float dt_entr = t1 - 0.5f * (p->t_begin + p->t_end); + const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase; p->force.pressure = (p->entropy + p->entropy_dt * dt_entr) * powf(p->rho, const_hydro_gamma); diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h index 5ac83e227c..46e156bb99 100644 --- a/src/hydro/Gadget2/hydro_debug.h +++ b/src/hydro/Gadget2/hydro_debug.h @@ -26,11 +26,11 @@ __attribute__((always_inline)) "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, " "dS/dt=%.3e, c=%.3e\n" "divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e] \n " - "v_sig=%e dh/dt=%.3e t_begin=%.3e, t_end=%.3e\n", + "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%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, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->force.pressure, p->entropy, p->entropy_dt, p->force.soundspeed, p->div_v, p->force.curl_v, p->density.rot_v[0], p->density.rot_v[1], - p->density.rot_v[2], p->force.v_sig, p->h_dt, p->t_begin, p->t_end); + p->density.rot_v[2], p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end); } diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h index 41898513aa..298cfc5ac1 100644 --- a/src/hydro/Gadget2/hydro_part.h +++ b/src/hydro/Gadget2/hydro_part.h @@ -47,10 +47,10 @@ struct part { float h_dt; /* Particle time of beginning of time-step. */ - float t_begin; + int ti_begin; /* Particle time of end of time-step. */ - float t_end; + int ti_end; /* Particle density. */ float rho; diff --git a/src/runner.c b/src/runner.c index 157616ad29..a50770dce7 100644 --- a/src/runner.c +++ b/src/runner.c @@ -40,12 +40,13 @@ #include "debug.h" #include "engine.h" #include "error.h" +#include "gravity.h" +#include "hydro.h" +#include "minmax.h" #include "scheduler.h" #include "space.h" #include "task.h" #include "timers.h" -#include "hydro.h" -#include "gravity.h" /* Orientation of the cell pairs */ const float runner_shift[13 * 3] = { @@ -496,7 +497,7 @@ void runner_doinit(struct runner *r, struct cell *c, int timer) { struct part *p, *parts = c->parts; const int count = c->count; - const float t_end = r->e->time; + const int ti_current = r->e->ti_current; TIMER_TIC; @@ -513,7 +514,7 @@ void runner_doinit(struct runner *r, struct cell *c, int timer) { /* Get a direct pointer on the part. */ p = &parts[i]; - if (p->t_end <= t_end) { + if (p->ti_end <= ti_current) { /* Get ready for a density calculation */ hydro_init_part(p); @@ -547,8 +548,9 @@ void runner_doghost(struct runner *r, struct cell *c) { int redo, count = c->count; int *pid; float h_corr; - float t_end = r->e->time; - + const int ti_current = r->e->ti_current; + const double timeBase = r->e->timeBase; + TIMER_TIC; /* Recurse? */ @@ -578,10 +580,10 @@ void runner_doghost(struct runner *r, struct cell *c) { xp = &xparts[pid[i]]; /* Is this part within the timestep? */ - if (p->t_end <= t_end) { + if (p->ti_end <= ti_current) { /* Finish the density calculation */ - hydro_end_density(p, t_end); + hydro_end_density(p, ti_current); //MATTHIEU /* If no derivative, double the smoothing length. */ if (p->density.wcount_dh == 0.0f) h_corr = p->h; @@ -618,7 +620,7 @@ void runner_doghost(struct runner *r, struct cell *c) { /* As of here, particle force variables will be set. */ /* Compute variables required for the force loop */ - hydro_prepare_force(p, xp, t_end); + hydro_prepare_force(p, xp, ti_current, timeBase); /* The particle force values are now set. Do _NOT_ try to read any particle density variables! */ @@ -696,7 +698,10 @@ void runner_doghost(struct runner *r, struct cell *c) { void runner_dodrift(struct runner *r, struct cell *c, int timer) { const int nr_parts = c->count; - const float dt = r->e->time - r->e->timeOld; + const double timeBase = r->e->timeBase; + const double dt = (r->e->ti_current - r->e->ti_old) * timeBase; + const float ti_old = r->e->ti_old; + const float ti_current = r->e->ti_current; struct part *restrict p, *restrict parts = c->parts; struct xpart *restrict xp, *restrict xparts = c->xparts; float dx_max = 0.f, h_max = 0.f; @@ -742,7 +747,7 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) { p->rho *= expf(w); /* Predict the values of the extra fields */ - hydro_predict_extra(p, xp, r->e->timeOld, r->e->time); + hydro_predict_extra(p, xp, ti_old, ti_current, timeBase); /* Compute motion since last cell construction */ const float dx = @@ -795,18 +800,21 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) { void runner_dokick(struct runner *r, struct cell *c, int timer) { - const float dt_max_timeline = r->e->timeEnd - r->e->timeBegin; - const float global_dt_min = r->e->dt_min, global_dt_max = r->e->dt_max; - const float t_current = r->e->time; + //const float dt_max_timeline = r->e->timeEnd - r->e->timeBegin; + const float global_dt_min = r->e->dt_min; + const float global_dt_max = r->e->dt_max; + const float ti_current = r->e->ti_current; + const double timeBase = r->e->timeBase; + const double timeBase_inv = 1.0 / r->e->timeBase; const int count = c->count; const int is_fixdt = (r->e->policy & engine_policy_fixdt) == engine_policy_fixdt; - float new_dt; - float dt_timeline; + int new_dti; + int dti_timeline; int updated = 0; - float t_end_min = FLT_MAX, t_end_max = 0.f; + int ti_end_min = max_nr_timesteps, ti_end_max = 0; double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0; float mom[3] = {0.0f, 0.0f, 0.0f}; float ang[3] = {0.0f, 0.0f, 0.0f}; @@ -832,7 +840,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { x[2] = p->x[2]; /* If particle needs to be kicked */ - if (is_fixdt || p->t_end <= t_current) { + if (is_fixdt || p->ti_end <= ti_current) { /* First, finish the force loop */ p->h_dt *= p->h * 0.333333333f; @@ -845,7 +853,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { if (is_fixdt) { /* Now we have a time step, proceed with the kick */ - new_dt = global_dt_max; + new_dti = global_dt_max * timeBase_inv; } else { @@ -853,7 +861,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { const float new_dt_hydro = hydro_compute_timestep(p, xp); const float new_dt_grav = gravity_compute_timestep(p, xp); - new_dt = fminf(new_dt_hydro, new_dt_grav); + float new_dt = fminf(new_dt_hydro, new_dt_grav); /* Limit change in h */ const float dt_h_change = @@ -862,33 +870,36 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { new_dt = fminf(new_dt, dt_h_change); - /* Recover the current timestep */ - const float current_dt = p->t_end - p->t_begin; - - /* Limit timestep increase */ - if (current_dt > 0.0f) new_dt = fminf(new_dt, 2.0f * current_dt); - /* Limit timestep within the allowed range */ new_dt = fminf(new_dt, global_dt_max); new_dt = fmaxf(new_dt, global_dt_min); + /* Convert to integer time */ + new_dti = new_dt * timeBase_inv; + + /* Recover the current timestep */ + const float current_dti = p->ti_end - p->ti_begin; + + /* Limit timestep increase */ + if (current_dti > 0) new_dti = min(new_dt, 2 * current_dti); + /* Put this timestep on the time line */ - dt_timeline = dt_max_timeline; - while (new_dt < dt_timeline) dt_timeline /= 2.; + dti_timeline = max_nr_timesteps; + while (new_dti < dti_timeline) dti_timeline /= 2; /* Now we have a time step, proceed with the kick */ - new_dt = dt_timeline; + new_dti = dti_timeline; } /* Compute the time step for this kick */ - const float t_start = 0.5f * (p->t_begin + p->t_end); - const float t_end = p->t_end + 0.5f * new_dt; - const float dt = t_end - t_start; - const float half_dt = t_end - p->t_end; + const int ti_start = (p->ti_begin + p->ti_end) / 2; + const int ti_end = p->ti_end + new_dti / 2; + const float dt = (ti_end - ti_start) * timeBase; + const float half_dt = (ti_end - p->ti_end) * timeBase; /* Move particle forward in time */ - p->t_begin = p->t_end; - p->t_end = p->t_begin + new_dt; + p->ti_begin = p->ti_end; + p->ti_end = p->ti_begin + new_dti; /* Kick particles in momentum space */ xp->v_full[0] += p->a_hydro[0] * dt; @@ -929,8 +940,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { e_int += hydro_get_internal_energy(p); /* Minimal time for next end of time-step */ - t_end_min = fminf(p->t_end, t_end_min); - t_end_max = fmaxf(p->t_end, t_end_max); + ti_end_min = min(p->ti_end, ti_end_min); + ti_end_max = max(p->ti_end, ti_end_max); /* Number of updated particles */ updated++; @@ -961,8 +972,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { ang[0] += cp->ang[0]; ang[1] += cp->ang[1]; ang[2] += cp->ang[2]; - t_end_min = fminf(cp->t_end_min, t_end_min); - t_end_max = fmaxf(cp->t_end_max, t_end_max); + ti_end_min = min(cp->ti_end_min, ti_end_min); + ti_end_max = max(cp->ti_end_max, ti_end_max); } } @@ -978,8 +989,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { c->ang[0] = ang[0]; c->ang[1] = ang[1]; c->ang[2] = ang[2]; - c->t_end_min = t_end_min; - c->t_end_max = t_end_max; + c->ti_end_min = ti_end_min; + c->ti_end_max = ti_end_max; if (timer) { #ifdef TIMER_VERBOSE @@ -1085,8 +1096,8 @@ void *runner_main(void *data) { case task_type_recv: parts = ci->parts; nr_parts = ci->count; - ci->t_end_min = ci->t_end_max = FLT_MAX; - for (k = 0; k < nr_parts; k++) parts[k].t_end = FLT_MAX; + ci->ti_end_min = ci->ti_end_max = max_nr_timesteps; + for (k = 0; k < nr_parts; k++) parts[k].ti_end = max_nr_timesteps; break; case task_type_grav_pp: if (t->cj == NULL) diff --git a/src/runner_doiact.h b/src/runner_doiact.h index 1b987eb140..5b29576c5a 100644 --- a/src/runner_doiact.h +++ b/src/runner_doiact.h @@ -111,7 +111,7 @@ void DOPAIR_NAIVE(struct runner *r, struct cell *restrict ci, struct part *restrict pi, *restrict pj; double pix[3]; float dx[3], hi, hig2, r2; - float t_current = e->time; + const int ti_current = e->ti_current; #ifdef VECTORIZE int icount = 0; float r2q[VEC_SIZE] __attribute__((aligned(16))); @@ -123,7 +123,7 @@ void DOPAIR_NAIVE(struct runner *r, struct cell *restrict ci, TIMER_TIC /* Anything to do here? */ - if (ci->t_end_min > t_current && cj->t_end_min > t_current) return; + if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; /* Get the relative distance between the pairs, wrapping. */ for (k = 0; k < 3; k++) { @@ -220,7 +220,7 @@ void DOSELF_NAIVE(struct runner *r, struct cell *restrict c) { struct part *restrict pi, *restrict pj; double pix[3] = {0.0, 0.0, 0.0}; float dx[3], hi, hig2, r2; - float t_current = r->e->time; + const int ti_current = r->e->ti_current; #ifdef VECTORIZE int icount = 0; float r2q[VEC_SIZE] __attribute__((aligned(16))); @@ -232,7 +232,7 @@ void DOSELF_NAIVE(struct runner *r, struct cell *restrict c) { TIMER_TIC /* Anything to do here? */ - if (c->t_end_min > t_current) return; + if (c->ti_end_min > ti_current) return; /* printf( "runner_dopair_naive: doing pair [ %g %g %g ]/[ %g %g %g ] with %i/%i parts and shift = [ %g %g %g ].\n" , @@ -770,7 +770,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) { double hi_max, hj_max; double di_max, dj_min; int count_i, count_j; - float t_current = e->time; + const int ti_current = e->ti_current; #ifdef VECTORIZE int icount = 0; float r2q[VEC_SIZE] __attribute__((aligned(16))); @@ -782,7 +782,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) { TIMER_TIC /* Anything to do here? */ - if (ci->t_end_min > t_current && cj->t_end_min > t_current) return; + if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; /* Get the sort ID. */ sid = space_getsid(e->s, &ci, &cj, shift); @@ -816,7 +816,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) { /* Get a hold of the ith part in ci. */ pi = &parts_i[sort_i[pid].i]; - if (pi->t_end > t_current) continue; + if (pi->ti_end > ti_current) continue; hi = pi->h; di = sort_i[pid].d + hi * kernel_gamma + dx_max - rshift; if (di < dj_min) continue; @@ -880,7 +880,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) { /* Get a hold of the jth part in cj. */ pj = &parts_j[sort_j[pjd].i]; - if (pj->t_end > t_current) continue; + if (pj->ti_end > ti_current) continue; hj = pj->h; dj = sort_j[pjd].d - hj * kernel_gamma - dx_max - rshift; if (dj > di_max) continue; @@ -967,7 +967,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { double hi_max, hj_max; double di_max, dj_min; int count_i, count_j; - float t_current = e->time; + const int ti_current = e->ti_current; #ifdef VECTORIZE int icount1 = 0; float r2q1[VEC_SIZE] __attribute__((aligned(16))); @@ -985,7 +985,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { TIMER_TIC /* Anything to do here? */ - if (ci->t_end_min > t_current && cj->t_end_min > t_current) return; + if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; /* Get the shift ID. */ sid = space_getsid(e->s, &ci, &cj, shift); @@ -1014,28 +1014,28 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { dx_max = (ci->dx_max + cj->dx_max); /* Collect the number of parts left and right below dt. */ - if (ci->t_end_max <= t_current) { + if (ci->ti_end_max <= ti_current) { sortdt_i = sort_i; countdt_i = count_i; - } else if (ci->t_end_min <= t_current) { + } else if (ci->ti_end_min <= ti_current) { if ((sortdt_i = (struct entry *)alloca(sizeof(struct entry) * count_i)) == NULL) error("Failed to allocate dt sortlists."); for (k = 0; k < count_i; k++) - if (parts_i[sort_i[k].i].t_end <= t_current) { + if (parts_i[sort_i[k].i].ti_end <= ti_current) { sortdt_i[countdt_i] = sort_i[k]; countdt_i += 1; } } - if (cj->t_end_max <= t_current) { + if (cj->ti_end_max <= ti_current) { sortdt_j = sort_j; countdt_j = count_j; - } else if (cj->t_end_min <= t_current) { + } else if (cj->ti_end_min <= ti_current) { if ((sortdt_j = (struct entry *)alloca(sizeof(struct entry) * count_j)) == NULL) error("Failed to allocate dt sortlists."); for (k = 0; k < count_j; k++) - if (parts_j[sort_j[k].i].t_end <= t_current) { + if (parts_j[sort_j[k].i].ti_end <= ti_current) { sortdt_j[countdt_j] = sort_j[k]; countdt_j += 1; } @@ -1055,7 +1055,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { for (k = 0; k < 3; k++) pix[k] = pi->x[k] - shift[k]; /* Look at valid dt parts only? */ - if (pi->t_end > t_current) { + if (pi->ti_end > ti_current) { /* Loop over the parts in cj within dt. */ for (pjd = 0; pjd < countdt_j && sortdt_j[pjd].d < di; pjd++) { @@ -1127,7 +1127,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { #ifndef VECTORIZE /* Does pj need to be updated too? */ - if (pj->t_end <= t_current) + if (pj->ti_end <= ti_current) IACT(r2, dx, hi, hj, pi, pj); else IACT_NONSYM(r2, dx, hi, hj, pi, pj); @@ -1135,7 +1135,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { #else /* Does pj need to be updated too? */ - if (pj->t_end <= t_current) { + if (pj->ti_end <= ti_current) { /* Add this interaction to the symmetric queue. */ r2q2[icount2] = r2; @@ -1200,7 +1200,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { hjg2 = hj * hj * kernel_gamma2; /* Is this particle outside the dt? */ - if (pj->t_end > t_current) { + if (pj->ti_end > ti_current) { /* Loop over the parts in ci. */ for (pid = countdt_i - 1; pid >= 0 && sortdt_i[pid].d > dj; pid--) { @@ -1271,7 +1271,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) { #ifndef VECTORIZE /* Does pi need to be updated too? */ - if (pi->t_end <= t_current) + if (pi->ti_end <= ti_current) IACT(r2, dx, hj, hi, pj, pi); else IACT_NONSYM(r2, dx, hj, hi, pj, pi); @@ -1361,7 +1361,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { double pix[3]; float dx[3], hi, hj, hig2, r2; struct part *restrict parts = c->parts, *restrict pi, *restrict pj; - float t_current = r->e->time; + const int ti_current = r->e->ti_current; int firstdt = 0, countdt = 0, *indt = NULL, doj; #ifdef VECTORIZE int icount1 = 0; @@ -1380,13 +1380,13 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { TIMER_TIC /* Set up indt if needed. */ - if (c->t_end_min > t_current) + if (c->ti_end_min > ti_current) return; - else if (c->t_end_max > t_current) { + else if (c->ti_end_max > ti_current) { if ((indt = (int *)alloca(sizeof(int) * count)) == NULL) error("Failed to allocate indt."); for (k = 0; k < count; k++) - if (parts[k].t_end <= t_current) { + if (parts[k].ti_end <= ti_current) { indt[countdt] = k; countdt += 1; } @@ -1404,7 +1404,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { hig2 = hi * hi * kernel_gamma2; /* Is the ith particle inactive? */ - if (pi->t_end > t_current) { + if (pi->ti_end > ti_current) { /* Loop over the other particles .*/ for (pjd = firstdt; pjd < countdt; pjd++) { @@ -1472,7 +1472,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) { dx[k] = pix[k] - pj->x[k]; r2 += dx[k] * dx[k]; } - doj = (pj->t_end <= t_current) && (r2 < hj * hj * kernel_gamma2); + doj = (pj->ti_end <= ti_current) && (r2 < hj * hj * kernel_gamma2); /* Hit or miss? */ if (r2 < hig2 || doj) { @@ -1584,7 +1584,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { double pix[3]; float dx[3], hi, hj, hig2, r2; struct part *restrict parts = c->parts, *restrict pi, *restrict pj; - float t_current = r->e->time; + const int ti_current = r->e->ti_current; int firstdt = 0, countdt = 0, *indt = NULL; #ifdef VECTORIZE int icount1 = 0; @@ -1603,13 +1603,13 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { TIMER_TIC /* Set up indt if needed. */ - if (c->t_end_min > t_current) + if (c->ti_end_min > ti_current) return; - else if (c->t_end_max > t_current) { + else if (c->ti_end_max > ti_current) { if ((indt = (int *)alloca(sizeof(int) * count)) == NULL) error("Failed to allocate indt."); for (k = 0; k < count; k++) - if (parts[k].t_end <= t_current) { + if (parts[k].ti_end <= ti_current) { indt[countdt] = k; countdt += 1; } @@ -1627,7 +1627,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { hig2 = hi * hi * kernel_gamma2; /* Is the ith particle not active? */ - if (pi->t_end > t_current) { + if (pi->ti_end > ti_current) { /* Loop over the other particles .*/ for (pjd = firstdt; pjd < countdt; pjd++) { @@ -1702,7 +1702,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { #ifndef VECTORIZE /* Does pj need to be updated too? */ - if (pj->t_end <= t_current) + if (pj->ti_end <= ti_current) IACT(r2, dx, hi, hj, pi, pj); else IACT_NONSYM(r2, dx, hi, hj, pi, pj); @@ -1710,7 +1710,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) { #else /* Does pj need to be updated too? */ - if (pj->t_end <= t_current) { + if (pj->ti_end <= ti_current) { /* Add this interaction to the symmetric queue. */ r2q2[icount2] = r2; @@ -1795,7 +1795,7 @@ void DOSUB1(struct runner *r, struct cell *ci, struct cell *cj, int sid, double shift[3]; float h; struct space *s = r->e->s; - float t_current = r->e->time; + const int ti_current = r->e->ti_current; TIMER_TIC @@ -1803,7 +1803,7 @@ void DOSUB1(struct runner *r, struct cell *ci, struct cell *cj, int sid, if (cj == NULL) { /* Should we even bother? */ - if (ci->t_end_min > t_current) return; + if (ci->ti_end_min > ti_current) return; /* Recurse? */ if (ci->split) { @@ -1829,7 +1829,7 @@ void DOSUB1(struct runner *r, struct cell *ci, struct cell *cj, int sid, else { /* Should we even bother? */ - if (ci->t_end_min > t_current && cj->t_end_min > t_current) return; + if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; /* Get the cell dimensions. */ h = fmin(ci->h[0], fmin(ci->h[1], ci->h[2])); @@ -2042,7 +2042,7 @@ void DOSUB1(struct runner *r, struct cell *ci, struct cell *cj, int sid, } /* Otherwise, compute the pair directly. */ - else if (ci->t_end_min <= t_current || cj->t_end_min <= t_current) { + else if (ci->ti_end_min <= ti_current || cj->ti_end_min <= ti_current) { /* Do any of the cells need to be sorted first? */ if (!(ci->sorted & (1 << sid))) runner_dosort(r, ci, (1 << sid), 1); @@ -2070,7 +2070,7 @@ void DOSUB2(struct runner *r, struct cell *ci, struct cell *cj, int sid, double shift[3]; float h; struct space *s = r->e->s; - float t_current = r->e->time; + const int ti_current = r->e->ti_current; TIMER_TIC @@ -2078,7 +2078,7 @@ void DOSUB2(struct runner *r, struct cell *ci, struct cell *cj, int sid, if (cj == NULL) { /* Should we even bother? */ - if (ci->t_end_min > t_current) return; + if (ci->ti_end_min > ti_current) return; /* Recurse? */ if (ci->split) { @@ -2104,7 +2104,7 @@ void DOSUB2(struct runner *r, struct cell *ci, struct cell *cj, int sid, else { /* Should we even bother? */ - if (ci->t_end_min > t_current && cj->t_end_min > t_current) return; + if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; /* Get the cell dimensions. */ h = fmin(ci->h[0], fmin(ci->h[1], ci->h[2])); @@ -2317,7 +2317,7 @@ void DOSUB2(struct runner *r, struct cell *ci, struct cell *cj, int sid, } /* Otherwise, compute the pair directly. */ - else if (ci->t_end_min <= t_current || cj->t_end_min <= t_current) { + else if (ci->ti_end_min <= ti_current || cj->ti_end_min <= ti_current) { /* Do any of the cells need to be sorted first? */ if (!(ci->sorted & (1 << sid))) runner_dosort(r, ci, (1 << sid), 1); @@ -2346,7 +2346,7 @@ void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts, float h; struct space *s = r->e->s; struct cell *sub = NULL; - float t_current = r->e->time; + const int ti_current = r->e->ti_current; TIMER_TIC @@ -2906,7 +2906,7 @@ void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts, } /* Otherwise, compute the pair directly. */ - else if (ci->t_end_min <= t_current || cj->t_end_min <= t_current) { + else if (ci->ti_end_min <= ti_current || cj->ti_end_min <= ti_current) { /* Get the relative distance between the pairs, wrapping. */ for (k = 0; k < 3; k++) { diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index dba2321365..4f4e855789 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -42,7 +42,7 @@ void runner_dopair_grav_new(struct runner *r, struct cell *ci, double pix[3]; float dx[3], r2, h_max, di, dj; int count_i, count_j, cnj, cnj_new; - float t_end = e->time; + const int ti_current = e->ti_current; struct multipole m; #ifdef VECTORIZE int icount = 0; @@ -53,7 +53,7 @@ void runner_dopair_grav_new(struct runner *r, struct cell *ci, TIMER_TIC /* Anything to do here? */ - if (ci->t_end_min > t_end && cj->t_end_min > t_end) return; + if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; /* Get the sort ID. */ sid = space_getsid(e->s, &ci, &cj, shift); @@ -93,7 +93,7 @@ void runner_dopair_grav_new(struct runner *r, struct cell *ci, /* Get a hold of the ith part in ci. */ pi = &parts_i[sort_i[pid].i]; - if (pi->t_end > t_end) continue; + if (pi->ti_end > ti_current) continue; di = sort_i[pid].d + h_max - rshift; for (k = 0; k < 3; k++) pix[k] = pi->x[k] - shift[k]; @@ -347,7 +347,7 @@ void runner_dopair_grav(struct runner *r, struct cell *restrict ci, struct gpart *restrict pi, *restrict pj; double pix[3]; float dx[3], r2; - float t_end = r->e->time; + const int ti_current = r->e->ti_current; #ifdef VECTORIZE int icount = 0; float r2q[VEC_SIZE] __attribute__((aligned(16))); @@ -357,7 +357,7 @@ void runner_dopair_grav(struct runner *r, struct cell *restrict ci, TIMER_TIC /* Anything to do here? */ - if (ci->t_end_min > t_end && cj->t_end_min > t_end) return; + if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return; /* Get the relative distance between the pairs, wrapping. */ if (e->s->periodic) @@ -454,7 +454,7 @@ void runner_doself_grav(struct runner *r, struct cell *restrict c) { struct gpart *restrict pi, *restrict pj; double pix[3] = {0.0, 0.0, 0.0}; float dx[3], r2; - float t_end = r->e->time; + const int ti_current = r->e->ti_current; #ifdef VECTORIZE int icount = 0; float r2q[VEC_SIZE] __attribute__((aligned(16))); @@ -464,7 +464,7 @@ void runner_doself_grav(struct runner *r, struct cell *restrict c) { TIMER_TIC /* Anything to do here? */ - if (c->t_end_min > t_end) return; + if (c->ti_end_min > ti_current) return; /* Loop over every part in c. */ for (pid = 0; pid < count; pid++) { diff --git a/src/space.c b/src/space.c index 60f86bcde8..a4f905ff57 100644 --- a/src/space.c +++ b/src/space.c @@ -41,6 +41,7 @@ #include "error.h" #include "kernel.h" #include "lock.h" +#include "minmax.h" #include "runner.h" /* Shared sort structure. */ @@ -976,7 +977,8 @@ void space_map_cells_pre(struct space *s, int full, void space_split(struct space *s, struct cell *c) { int k, count = c->count, gcount = c->gcount, maxdepth = 0; - float h, h_max = 0.0f, t_end_min = FLT_MAX, t_end_max = 0., t_end; + float h, h_max = 0.0f; + int ti_end_min = max_nr_timesteps, ti_end_max = 0, ti_end; struct cell *temp; struct part *p, *parts = c->parts; struct xpart *xp, *xparts = c->xparts; @@ -1025,16 +1027,16 @@ void space_split(struct space *s, struct cell *c) { } else { space_split(s, c->progeny[k]); h_max = fmaxf(h_max, c->progeny[k]->h_max); - t_end_min = fminf(t_end_min, c->progeny[k]->t_end_min); - t_end_max = fmaxf(t_end_max, c->progeny[k]->t_end_max); + ti_end_min = min(ti_end_min, c->progeny[k]->ti_end_min); + ti_end_max = max(ti_end_max, c->progeny[k]->ti_end_max); if (c->progeny[k]->maxdepth > maxdepth) maxdepth = c->progeny[k]->maxdepth; } /* Set the values for this cell. */ c->h_max = h_max; - c->t_end_min = t_end_min; - c->t_end_max = t_end_max; + c->ti_end_min = ti_end_min; + c->ti_end_max = ti_end_max; c->maxdepth = maxdepth; } @@ -1056,14 +1058,14 @@ void space_split(struct space *s, struct cell *c) { xp->x_old[1] = p->x[1]; xp->x_old[2] = p->x[2]; h = p->h; - t_end = p->t_end; + ti_end = p->ti_end; if (h > h_max) h_max = h; - if (t_end < t_end_min) t_end_min = t_end; - if (t_end > t_end_max) t_end_max = t_end; + if (ti_end < ti_end_min) ti_end_min = ti_end; + if (ti_end > ti_end_max) ti_end_max = ti_end; } c->h_max = h_max; - c->t_end_min = t_end_min; - c->t_end_max = t_end_max; + c->ti_end_min = ti_end_min; + c->ti_end_max = ti_end_max; } /* Set ownership according to the start of the parts array. */ -- GitLab