Skip to content
Snippets Groups Projects
Commit 8c4afc4a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Common part of time-step calculation in a separate function.

parent 47aba411
No related branches found
No related tags found
1 merge request!167Kick task for fixdt and proper treatment of conserved quantities
Hide whitespace changes
Inline Side-by-side
src/engine.c
+ 2
− 0
View file @ 8c4afc4a
......@@ -2455,6 +2455,7 @@ void engine_init(struct engine *e, struct space *s,
e->ti_current = 0;
e->timeStep = 0.;
e->timeBase = 0.;
e->timeBase_inv = 0.;
e->timeFirstSnapshot =
parser_get_param_double(params, "Snapshots:time_first");
e->deltaTimeSnapshot =
......@@ -2596,6 +2597,7 @@ void engine_init(struct engine *e, struct space *s,
/* Deal with timestep */
e->timeBase = (e->timeEnd - e->timeBegin) / max_nr_timesteps;
e->timeBase_inv = 1.0 / e->timeBase;
e->ti_current = 0;
/* Fixed time-step case */
......
src/engine.h
+ 1
− 0
View file @ 8c4afc4a
......@@ -132,6 +132,7 @@ struct engine {
/* Time base */
double timeBase;
double timeBase_inv;
/* Snapshot information */
double timeFirstSnapshot;
......
src/runner.c
+ 4
− 14
View file @ 8c4afc4a
......@@ -1079,20 +1079,14 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
*/
void runner_do_kick(struct runner *r, struct cell *c, int timer) {
const double global_dt_min = r->e->dt_min;
const double global_dt_max = r->e->dt_max;
const struct engine *e = r->e;
const double timeBase = e->timeBase;
const int 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 gcount = c->gcount;
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
const struct external_potential *potential = r->e->external_potential;
const struct hydro_props *hydro_properties = r->e->hydro_properties;
const struct phys_const *constants = r->e->physical_constants;
const float ln_max_h_change = hydro_properties->log_max_h_change;
int updated = 0, g_updated = 0;
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
......@@ -1124,9 +1118,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
gravity_end_force(gp);
/* Compute the next timestep */
const int new_dti =
get_gpart_timestep(gp, potential, constants, global_dt_min,
global_dt_max, timeBase_inv);
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);
......@@ -1161,9 +1153,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
if (p->gpart != NULL) gravity_end_force(p->gpart);
/* Compute the next timestep (hydro condition) */
const int new_dti = get_part_timestep(
p, xp, hydro_properties, potential, constants, global_dt_min,
global_dt_max, timeBase_inv, ln_max_h_change);
const int new_dti = get_part_timestep(p, xp, e);
/* Now we have a time step, proceed with the kick */
kick_part(p, xp, new_dti, timeBase);
......
src/timestep.h
+ 38
− 47
View file @ 8c4afc4a
......@@ -26,26 +26,14 @@
#include "const.h"
#include "debug.h"
__attribute__((always_inline)) INLINE static int get_gpart_timestep(
const struct gpart *gp, const struct external_potential *potential,
const struct phys_const *constants, double global_dt_min,
double global_dt_max, double timeBase_inv) {
const float new_dt_external =
gravity_compute_timestep_external(potential, constants, gp);
const float new_dt_self = gravity_compute_timestep_self(constants, gp);
float new_dt = fminf(new_dt_external, new_dt_self);
/* Limit timestep within the allowed range */
new_dt = fminf(new_dt, global_dt_max);
new_dt = fmaxf(new_dt, global_dt_min);
__attribute__((always_inline)) INLINE static int get_integer_timestep(
float new_dt, int ti_begin, int ti_end, double timeBase_inv) {
/* Convert to integer time */
int new_dti = new_dt * timeBase_inv;
int new_dti = (int)(new_dt * timeBase_inv);
/* Recover the current timestep */
const int current_dti = gp->ti_end - gp->ti_begin;
const int current_dti = ti_end - ti_begin;
/* Limit timestep increase */
if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
......@@ -58,30 +46,47 @@ __attribute__((always_inline)) INLINE static int get_gpart_timestep(
/* Make sure we are allowed to increase the timestep size */
if (new_dti > current_dti) {
if ((max_nr_timesteps - gp->ti_end) % new_dti > 0) new_dti = current_dti;
if ((max_nr_timesteps - ti_end) % new_dti > 0) new_dti = current_dti;
}
return new_dti;
}
__attribute__((always_inline)) INLINE static int get_gpart_timestep(
const struct gpart *gp, const struct engine *e) {
const float new_dt_external = gravity_compute_timestep_external(
e->external_potential, e->physical_constants, gp);
const float new_dt_self =
gravity_compute_timestep_self(e->physical_constants, gp);
float new_dt = fminf(new_dt_external, new_dt_self);
/* Limit timestep within the allowed range */
new_dt = fminf(new_dt, e->dt_max);
new_dt = fmaxf(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);
return new_dti;
}
__attribute__((always_inline)) INLINE static int get_part_timestep(
const struct part *p, const struct xpart *xp,
const struct hydro_props *hydro_properties,
const struct external_potential *potential,
const struct phys_const *constants, double global_dt_min,
double global_dt_max, double timeBase_inv, double ln_max_h_change) {
const struct part *p, const struct xpart *xp, const struct engine *e) {
/* Compute the next timestep (hydro condition) */
const float new_dt_hydro = hydro_compute_timestep(p, xp, hydro_properties);
const float new_dt_hydro = hydro_compute_timestep(p, xp, e->hydro_properties);
/* Compute the next timestep (gravity condition) */
float new_dt_grav = FLT_MAX;
if (p->gpart != NULL) {
const float new_dt_external =
gravity_compute_timestep_external(potential, constants, p->gpart);
const float new_dt_external = gravity_compute_timestep_external(
e->external_potential, e->physical_constants, p->gpart);
const float new_dt_self =
gravity_compute_timestep_self(constants, p->gpart);
gravity_compute_timestep_self(e->physical_constants, p->gpart);
new_dt_grav = fminf(new_dt_external, new_dt_self);
}
......@@ -91,33 +96,19 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
/* Limit change in h */
const float dt_h_change =
(p->h_dt != 0.0f) ? fabsf(ln_max_h_change * p->h / p->h_dt) : FLT_MAX;
(p->h_dt != 0.0f)
? fabsf(e->hydro_properties->log_max_h_change * p->h / p->h_dt)
: FLT_MAX;
new_dt = fminf(new_dt, dt_h_change);
/* Limit timestep within the allowed range */
new_dt = fminf(new_dt, global_dt_max);
new_dt = fmaxf(new_dt, global_dt_min);
new_dt = fminf(new_dt, e->dt_max);
new_dt = fmaxf(new_dt, e->dt_min);
/* Convert to integer time */
int new_dti = new_dt * timeBase_inv;
/* Recover the current timestep */
const int current_dti = p->ti_end - p->ti_begin;
/* Limit timestep increase */
if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
/* Put this timestep on the time line */
int dti_timeline = max_nr_timesteps;
while (new_dti < dti_timeline) dti_timeline /= 2;
new_dti = dti_timeline;
/* Make sure we are allowed to increase the timestep size */
if (new_dti > current_dti) {
if ((max_nr_timesteps - p->ti_end) % new_dti > 0) new_dti = current_dti;
}
const int new_dti =
get_integer_timestep(new_dt, p->ti_begin, p->ti_end, e->timeBase_inv);
return new_dti;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment