Commit 47aba411 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Timestep computation moved to a separate file

parent 5ee2948e
......@@ -48,6 +48,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h kick.h \
timestep.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
......
......@@ -25,7 +25,7 @@
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
struct part* p, struct xpart* xp,
const struct part* p, const struct xpart* xp,
const struct hydro_props* hydro_properties) {
const float CFL_condition = hydro_properties->CFL_condition;
......
......@@ -53,6 +53,7 @@
#include "space.h"
#include "task.h"
#include "timers.h"
#include "timestep.h"
/* Orientation of the cell pairs */
const float runner_shift[13 * 3] = {
......@@ -1122,38 +1123,10 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
/* First, finish the force calculation */
gravity_end_force(gp);
/* Compute the next timestep (gravity condition) */
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);
/* Convert to integer time */
int new_dti = new_dt * timeBase_inv;
/* Recover the current timestep */
const int current_dti = gp->ti_end - gp->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 - gp->ti_end) % new_dti > 0)
new_dti = current_dti;
}
/* Compute the next timestep */
const int new_dti =
get_gpart_timestep(gp, potential, constants, global_dt_min,
global_dt_max, timeBase_inv);
/* Now we have a time step, proceed with the kick */
kick_gpart(gp, new_dti, timeBase);
......@@ -1188,55 +1161,9 @@ 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 float new_dt_hydro =
hydro_compute_timestep(p, xp, 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_self =
gravity_compute_timestep_self(constants, p->gpart);
new_dt_grav = fminf(new_dt_external, new_dt_self);
}
/* Final time-step is minimum of hydro and gravity */
float new_dt = fminf(new_dt_hydro, new_dt_grav);
/* 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;
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);
/* 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_part_timestep(
p, xp, hydro_properties, potential, constants, global_dt_min,
global_dt_max, timeBase_inv, ln_max_h_change);
/* Now we have a time step, proceed with the kick */
kick_part(p, xp, new_dti, timeBase);
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_TIMESTEP_H
#define SWIFT_TIMESTEP_H
/* Config parameters. */
#include "../config.h"
/* Local headers. */
#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);
/* Convert to integer time */
int new_dti = new_dt * timeBase_inv;
/* Recover the current timestep */
const int current_dti = gp->ti_end - gp->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 - gp->ti_end) % new_dti > 0) new_dti = current_dti;
}
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) {
/* Compute the next timestep (hydro condition) */
const float new_dt_hydro = hydro_compute_timestep(p, xp, 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_self =
gravity_compute_timestep_self(constants, p->gpart);
new_dt_grav = fminf(new_dt_external, new_dt_self);
}
/* Final time-step is minimum of hydro and gravity */
float new_dt = fminf(new_dt_hydro, new_dt_grav);
/* 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;
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);
/* 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;
}
return new_dti;
}
#endif /* SWIFT_TIMESTEP_H */
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