Skip to content
GitLab
Commit f4fef206 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Propagated the changes everywhere

parent c8def96b
Hide whitespace changes
Inline Side-by-side
src/cell.c
View file @ f4fef206
......@@ -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;
}
/**
......
src/cell.h
View file @ f4fef206
......@@ -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;
......
src/debug.c
View file @ f4fef206
......@@ -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;
}
......
src/engine.c
View file @ f4fef206
......@@ -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
......
src/engine.h
View file @ f4fef206
......@@ -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;
......
src/gravity/Default/gravity_part.h
View file @ f4fef206
......@@ -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 {
......
src/hydro/Gadget2/hydro.h
View file @ f4fef206
......@@ -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);
......
src/hydro/Gadget2/hydro_debug.h
View file @ f4fef206
......@@ -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);
}
src/hydro/Gadget2/hydro_part.h
View file @ f4fef206
......@@ -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;
......
src/runner.c
View file @ f4fef206
......@@ -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;