Commit 14c0fe90 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Updated the kick1 and kick2 task to use the cosmological terms.

parent 839047c7
......@@ -180,7 +180,7 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
const struct cooling_function_data *cooling_func = e->cooling_func;
const struct phys_const *constants = e->physical_constants;
const struct unit_system *us = e->internal_units;
const double timeBase = e->timeBase;
const double time_base = e->time_base;
TIMER_TIC;
......@@ -202,8 +202,10 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
if (part_is_active(p, e)) {
// MATTHIEU -- cosmological dt for cooling??? */
/* Let's cool ! */
const double dt = get_timestep(p->time_bin, timeBase);
const double dt = get_timestep(p->time_bin, time_base);
cooling_cool_part(constants, us, cooling_func, p, xp, dt);
}
}
......@@ -948,6 +950,8 @@ void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer) {
void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const int with_cosmology = (e->policy & engine_policy_cosmology);
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
......@@ -956,7 +960,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
const int gcount = c->gcount;
const int scount = c->scount;
const integertime_t ti_current = e->ti_current;
const double timeBase = e->timeBase;
const double time_base = e->time_base;
TIMER_TIC;
......@@ -994,22 +998,31 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
ti_end, ti_begin, ti_step, p->time_bin, ti_current);
#endif
/* Time interval for this half-kick */
double dt_kick_grav, dt_kick_hydro, dt_kick_therm;
if (with_cosmology) {
dt_kick_hydro = cosmology_get_hydro_kick_factor(
cosmo, ti_begin, ti_begin + ti_step / 2);
dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_begin,
ti_begin + ti_step / 2);
dt_kick_therm = cosmology_get_therm_kick_factor(
cosmo, ti_begin, ti_begin + ti_step / 2);
} else {
dt_kick_hydro = (ti_step / 2) * time_base;
dt_kick_grav = (ti_step / 2) * time_base;
dt_kick_therm = (ti_step / 2) * time_base;
}
/* do the kick */
kick_part(p, xp, ti_begin, ti_begin + ti_step / 2, timeBase);
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, ti_begin,
ti_begin + ti_step / 2);
/* Update the accelerations to be used in the drift for hydro */
if (p->gpart != NULL) {
const float a_tot[3] = {p->a_hydro[0] + p->gpart->a_grav[0],
p->a_hydro[1] + p->gpart->a_grav[1],
p->a_hydro[2] + p->gpart->a_grav[2]};
p->a_hydro[0] = a_tot[0];
p->a_hydro[1] = a_tot[1];
p->a_hydro[2] = a_tot[2];
p->gpart->a_grav[0] = a_tot[0];
p->gpart->a_grav[1] = a_tot[1];
p->gpart->a_grav[2] = a_tot[2];
xp->a_grav[0] = p->gpart->a_grav[0];
xp->a_grav[1] = p->gpart->a_grav[1];
xp->a_grav[2] = p->gpart->a_grav[2];
}
}
}
......@@ -1038,8 +1051,17 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
ti_end, ti_begin, ti_step, gp->time_bin, ti_current);
#endif
/* Time interval for this half-kick */
double dt_kick_grav;
if (with_cosmology) {
dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_begin,
ti_begin + ti_step / 2);
} else {
dt_kick_grav = (ti_step / 2) * time_base;
}
/* do the kick */
kick_gpart(gp, ti_begin, ti_begin + ti_step / 2, timeBase);
kick_gpart(gp, dt_kick_grav, ti_begin, ti_begin + ti_step / 2);
}
}
......@@ -1067,8 +1089,17 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
ti_end, ti_begin, ti_step, sp->time_bin, ti_current);
#endif
/* Time interval for this half-kick */
double dt_kick_grav;
if (with_cosmology) {
dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_begin,
ti_begin + ti_step / 2);
} else {
dt_kick_grav = (ti_step / 2) * time_base;
}
/* do the kick */
kick_spart(sp, ti_begin, ti_begin + ti_step / 2, timeBase);
kick_spart(sp, dt_kick_grav, ti_begin, ti_begin + ti_step / 2);
}
}
}
......@@ -1088,8 +1119,8 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const integertime_t ti_current = e->ti_current;
const double timeBase = e->timeBase;
const struct cosmology *cosmo = e->cosmology;
const int with_cosmology = (e->policy & engine_policy_cosmology);
const int count = c->count;
const int gcount = c->gcount;
const int scount = c->scount;
......@@ -1097,6 +1128,8 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
struct spart *restrict sparts = c->sparts;
const integertime_t ti_current = e->ti_current;
const double time_base = e->time_base;
TIMER_TIC;
......@@ -1130,9 +1163,24 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
"time_bin=%d ti_current=%lld",
ti_begin, ti_step, p->time_bin, ti_current);
#endif
/* Time interval for this half-kick */
double dt_kick_grav, dt_kick_hydro, dt_kick_therm;
if (with_cosmology) {
dt_kick_hydro = cosmology_get_hydro_kick_factor(
cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
dt_kick_grav = cosmology_get_grav_kick_factor(
cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
dt_kick_therm = cosmology_get_therm_kick_factor(
cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
} else {
dt_kick_hydro = (ti_step / 2) * time_base;
dt_kick_grav = (ti_step / 2) * time_base;
dt_kick_therm = (ti_step / 2) * time_base;
}
/* Finish the time-step with a second half-kick */
kick_part(p, xp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm,
ti_begin + ti_step / 2, ti_begin + ti_step);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
......@@ -1162,8 +1210,18 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
error("Particle in wrong time-bin");
#endif
/* Time interval for this half-kick */
double dt_kick_grav;
if (with_cosmology) {
dt_kick_grav = cosmology_get_grav_kick_factor(
cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
} else {
dt_kick_grav = (ti_step / 2) * time_base;
}
/* Finish the time-step with a second half-kick */
kick_gpart(gp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
kick_gpart(gp, dt_kick_grav, ti_begin + ti_step / 2,
ti_begin + ti_step);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
......@@ -1194,8 +1252,18 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
error("Particle in wrong time-bin");
#endif
/* Time interval for this half-kick */
double dt_kick_grav;
if (with_cosmology) {
dt_kick_grav = cosmology_get_grav_kick_factor(
cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
} else {
dt_kick_grav = (ti_step / 2) * time_base;
}
/* Finish the time-step with a second half-kick */
kick_spart(sp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
kick_spart(sp, dt_kick_grav, ti_begin + ti_step / 2,
ti_begin + ti_step);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
......@@ -1230,7 +1298,6 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
struct spart *restrict sparts = c->sparts;
const double timeBase = e->timeBase;
TIMER_TIC;
......@@ -1277,10 +1344,6 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
p->time_bin = get_time_bin(ti_new_step);
if (p->gpart != NULL) p->gpart->time_bin = p->time_bin;
/* Tell the particle what the new physical time step is */
float dt = get_timestep(p->time_bin, timeBase);
hydro_timestep_extra(p, dt);
/* Number of updated particles */
updated++;
if (p->gpart != NULL) g_updated++;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment