/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
* Matthieu Schaller (schaller@strw.leidenuniv.nl)
* 2015 Peter W. Draper (p.w.draper@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/>.
*
******************************************************************************/
/* Config parameters. */
#include <config.h>
/* This object's header. */
#include "runner.h"
/* Local headers. */
#include "active.h"
#include "black_holes.h"
#include "cell.h"
#include "engine.h"
#include "feedback.h"
#include "kick.h"
#include "multipole.h"
#include "timers.h"
#include "timestep.h"
#include "timestep_limiter.h"
#include "timestep_sync.h"
#include "tracers.h"
/**
* @brief Initialize the multipoles before the gravity calculation.
*
* @param r The runner thread.
* @param c The cell.
* @param timer 1 if the time is to be recorded.
*/
void runner_do_init_grav(struct runner *r, struct cell *c, const int timer) {
const struct engine *e = r->e;
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
if (!(e->policy & engine_policy_self_gravity))
error("Grav-init task called outside of self-gravity calculation");
#endif
/* Anything to do here? */
if (!cell_is_active_gravity(c, e)) return;
/* Does the multipole need drifting? */
if (c->grav.ti_old_multipole < e->ti_current) cell_drift_multipole(c, e);
/* Reset the gravity acceleration tensors */
gravity_field_tensors_init(&c->grav.multipole->pot, e->ti_current);
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL)
runner_do_init_grav(r, c->progeny[k], /*timer=*/0);
}
}
if (timer) TIMER_TOC(timer_init_grav);
}
/**
* @brief Perform the first half-kick on all the active particles in a cell.
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
void runner_do_kick1(struct runner *r, struct cell *c, const int timer) {
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const struct hydro_props *hydro_props = e->hydro_properties;
const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
const int periodic = e->s->periodic;
const int with_cosmology = (e->policy & engine_policy_cosmology);
struct part *restrict parts = c->hydro.parts;
struct xpart *restrict xparts = c->hydro.xparts;
struct gpart *restrict gparts = c->grav.parts;
struct spart *restrict sparts = c->stars.parts;
struct sink *restrict sinks = c->sinks.parts;
struct bpart *restrict bparts = c->black_holes.parts;
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
const int sink_count = c->sinks.count;
const int bcount = c->black_holes.count;
const integertime_t ti_current = e->ti_current;
const double time_base = e->time_base;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_starting_hydro(c, e) && !cell_is_starting_gravity(c, e) &&
!cell_is_starting_stars(c, e) && !cell_is_starting_sinks(c, e) &&
!cell_is_starting_black_holes(c, e))
return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_kick1(r, c->progeny[k], /*timer=*/0);
} else {
integertime_t ti_begin_mesh = -1;
integertime_t ti_end_mesh = 0;
double dt_kick_mesh_grav = 0.;
/* Are we at a step where we do mesh-gravity time-integration? */
if (periodic && e->mesh->ti_beg_mesh_next == e->ti_current) {
const integertime_t ti_step =
e->mesh->ti_end_mesh_next - e->mesh->ti_beg_mesh_next;
ti_begin_mesh = e->mesh->ti_beg_mesh_next;
ti_end_mesh = e->mesh->ti_beg_mesh_next + ti_step / 2;
/* Time interval for this mesh gravity half-kick */
dt_kick_mesh_grav = kick_get_grav_kick_dt(
ti_begin_mesh, ti_end_mesh, time_base, with_cosmology, cosmo);
}
/* Loop over the parts in this cell. */
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
struct part *restrict p = &parts[k];
struct xpart *restrict xp = &xparts[k];
/* If particle needs to be kicked */
if (part_is_starting(p, e)) {
#ifdef SWIFT_DEBUG_CHECKS
if (p->limiter_data.wakeup != time_bin_not_awake)
error("Woken-up particle that has not been processed in kick1");
#endif
const integertime_t ti_step = get_integer_timestep(p->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current + 1, p->time_bin);
const integertime_t ti_end = ti_begin + ti_step / 2;
#ifdef SWIFT_DEBUG_CHECKS
if (ti_begin != ti_current)
error(
"Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
"ti_step=%lld time_bin=%d wakeup=%d ti_current=%lld",
ti_end, ti_begin, ti_step, p->time_bin, p->limiter_data.wakeup,
ti_current);
#endif
/* Time intervals for this half-kick */
const double dt_kick_grav = kick_get_grav_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
const double dt_kick_hydro = kick_get_hydro_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
const double dt_kick_therm = kick_get_therm_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
const double dt_kick_corr = kick_get_corr_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
/* Do the kick */
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_mesh_grav,
dt_kick_therm, dt_kick_corr, cosmo, hydro_props,
entropy_floor, ti_begin, ti_end, ti_begin_mesh, ti_end_mesh);
/* Update the accelerations to be used in the drift for hydro */
if (p->gpart != NULL) {
xp->a_grav[0] = p->gpart->a_grav[0] + p->gpart->a_grav_mesh[0];
xp->a_grav[1] = p->gpart->a_grav[1] + p->gpart->a_grav_mesh[1];
xp->a_grav[2] = p->gpart->a_grav[2] + p->gpart->a_grav_mesh[2];
}
}
}
/* Loop over the gparts in this cell. */
for (int k = 0; k < gcount; k++) {
/* Get a handle on the part. */
struct gpart *restrict gp = &gparts[k];
#ifdef SWIFT_DEBUG_CHECKS
if (ti_begin_mesh != -1 && !gpart_is_starting(gp, e) &&
!gpart_is_inhibited(gp, e)) {
error(
"Particle on a time-step longer than the mesh synchronization "
"step!");
}
#endif
/* If the g-particle has no counterpart and needs to be kicked */
if ((gp->type == swift_type_dark_matter ||
gp->type == swift_type_dark_matter_background ||
gp->type == swift_type_neutrino) &&
gpart_is_starting(gp, e)) {
const integertime_t ti_step = get_integer_timestep(gp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current + 1, gp->time_bin);
const integertime_t ti_end = ti_begin + ti_step / 2;
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_end_check =
get_integer_time_end(ti_current + 1, gp->time_bin);
if (ti_begin != ti_current)
error(
"G-particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
"ti_step=%lld time_bin=%d ti_current=%lld",
ti_end_check, ti_begin, ti_step, gp->time_bin, ti_current);
#endif
/* Time interval for this gravity half-kick */
const double dt_kick_grav = kick_get_grav_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
/* Do the kick */
kick_gpart(gp, dt_kick_grav, ti_begin, ti_end, dt_kick_mesh_grav,
ti_begin_mesh, ti_end_mesh);
}
}
/* Loop over the stars particles in this cell. */
for (int k = 0; k < scount; k++) {
/* Get a handle on the s-part. */
struct spart *restrict sp = &sparts[k];
/* If particle needs to be kicked */
if (spart_is_starting(sp, e)) {
const integertime_t ti_step = get_integer_timestep(sp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current + 1, sp->time_bin);
const integertime_t ti_end = ti_begin + ti_step / 2;
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_end_check =
get_integer_time_end(ti_current + 1, sp->time_bin);
if (ti_begin != ti_current)
error(
"S-particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
"ti_step=%lld time_bin=%d ti_current=%lld",
ti_end_check, ti_begin, ti_step, sp->time_bin, ti_current);
#endif
/* Time interval for this gravity half-kick */
const double dt_kick_grav = kick_get_grav_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
/* Do the kick */
kick_spart(sp, dt_kick_grav, ti_begin, ti_end, dt_kick_mesh_grav,
ti_begin_mesh, ti_end_mesh);
}
}
/* Loop over the sink particles in this cell. */
for (int k = 0; k < sink_count; k++) {
/* Get a handle on the s-part. */
struct sink *restrict sink = &sinks[k];
/* If particle needs to be kicked */
if (sink_is_starting(sink, e)) {
const integertime_t ti_step = get_integer_timestep(sink->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current + 1, sink->time_bin);
const integertime_t ti_end = ti_begin + ti_step / 2;
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_end_check =
get_integer_time_end(ti_current + 1, sink->time_bin);
if (ti_begin != ti_current)
error(
"Sink-particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
"ti_step=%lld time_bin=%d ti_current=%lld",
ti_end_check, ti_begin, ti_step, sink->time_bin, ti_current);
#endif
/* Time interval for this gravity half-kick */
const double dt_kick_grav = kick_get_grav_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
/* Do the kick */
kick_sink(sink, dt_kick_grav, ti_begin, ti_end, dt_kick_mesh_grav,
ti_begin_mesh, ti_end_mesh);
}
}
/* Loop over the black hole particles in this cell. */
for (int k = 0; k < bcount; k++) {
/* Get a handle on the s-part. */
struct bpart *restrict bp = &bparts[k];
/* If particle needs to be kicked */
if (bpart_is_starting(bp, e)) {
const integertime_t ti_step = get_integer_timestep(bp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current + 1, bp->time_bin);
const integertime_t ti_end = ti_begin + ti_step / 2;
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_end_check =
get_integer_time_end(ti_current + 1, bp->time_bin);
if (ti_begin != ti_current)
error(
"B-particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
"ti_step=%lld time_bin=%d ti_current=%lld",
ti_end_check, ti_begin, ti_step, bp->time_bin, ti_current);
#endif
/* Time interval for this gravity half-kick */
const double dt_kick_grav = kick_get_grav_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
/* Do the kick */
kick_bpart(bp, dt_kick_grav, ti_begin, ti_end, dt_kick_mesh_grav,
ti_begin_mesh, ti_end_mesh);
}
}
}
if (timer) TIMER_TOC(timer_kick1);
}
/**
* @brief Perform the second half-kick on all the active particles in a cell.
*
* Also prepares particles to be drifted.
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
void runner_do_kick2(struct runner *r, struct cell *c, const int timer) {
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const struct hydro_props *hydro_props = e->hydro_properties;
const struct pressure_floor_props *pressure_floor = e->pressure_floor_props;
const struct entropy_floor_properties *entropy_floor = e->entropy_floor;
const int with_cosmology = (e->policy & engine_policy_cosmology);
const int periodic = e->s->periodic;
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
const int sink_count = c->sinks.count;
const int bcount = c->black_holes.count;
struct part *restrict parts = c->hydro.parts;
struct xpart *restrict xparts = c->hydro.xparts;
struct gpart *restrict gparts = c->grav.parts;
struct spart *restrict sparts = c->stars.parts;
struct sink *restrict sinks = c->sinks.parts;
struct bpart *restrict bparts = c->black_holes.parts;
const integertime_t ti_current = e->ti_current;
const double time_base = e->time_base;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e) &&
!cell_is_active_stars(c, e) && !cell_is_active_sinks(c, e) &&
!cell_is_active_black_holes(c, e))
return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_kick2(r, c->progeny[k], /*timer=*/0);
} else {
integertime_t ti_begin_mesh = -1;
integertime_t ti_end_mesh = 0;
double dt_kick_mesh_grav = 0.;
/* Are we at a step where we do mesh-gravity time-integration? */
if (periodic && e->mesh->ti_end_mesh_last == e->ti_current) {
const integertime_t ti_step =
e->mesh->ti_end_mesh_last - e->mesh->ti_beg_mesh_last;
ti_begin_mesh = e->mesh->ti_beg_mesh_last + ti_step / 2;
ti_end_mesh = e->mesh->ti_end_mesh_last;
/* Time interval for this mesh gravity half-kick */
dt_kick_mesh_grav = kick_get_grav_kick_dt(
ti_begin_mesh, ti_end_mesh, time_base, with_cosmology, cosmo);
}
/* Loop over the particles in this cell. */
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
struct part *restrict p = &parts[k];
struct xpart *restrict xp = &xparts[k];
/* If particle needs to be kicked */
if (part_is_active(p, e)) {
#ifdef SWIFT_DEBUG_CHECKS
if (p->limiter_data.wakeup != time_bin_not_awake)
error("Woken-up particle that has not been processed in kick1");
#endif
/* Time-step length on the integer timeline */
const integertime_t ti_step = get_integer_timestep(p->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current, p->time_bin) + ti_step / 2;
const integertime_t ti_end = ti_begin + ti_step / 2;
#ifdef SWIFT_DEBUG_CHECKS
if (ti_end != ti_current)
error(
"Particle in wrong time-bin, ti_begin=%lld, ti_step=%lld "
"time_bin=%d wakeup=%d ti_current=%lld",
ti_begin, ti_step, p->time_bin, p->limiter_data.wakeup,
ti_current);
#endif
/* Time intervals for this half-kick */
const double dt_kick_grav = kick_get_grav_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
const double dt_kick_hydro = kick_get_hydro_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
const double dt_kick_therm = kick_get_therm_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
const double dt_kick_corr = kick_get_corr_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
/* Finish the time-step with a second half-kick */
kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_mesh_grav,
dt_kick_therm, dt_kick_corr, cosmo, hydro_props,
entropy_floor, ti_begin, ti_end, ti_begin_mesh, ti_end_mesh);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
if (p->ti_drift != p->ti_kick) error("Error integrating part in time.");
#endif
/* Prepare the values to be drifted */
hydro_reset_predicted_values(p, xp, cosmo, pressure_floor);
}
}
/* Loop over the g-particles in this cell. */
for (int k = 0; k < gcount; k++) {
/* Get a handle on the part. */
struct gpart *restrict gp = &gparts[k];
#ifdef SWIFT_DEBUG_CHECKS
if (ti_begin_mesh != -1 && !gpart_is_active(gp, e)) {
error(
"Particle on a time-step longer than the mesh synchronization "
"step!");
}
#endif
/* If the g-particle has no counterpart and needs to be kicked */
if ((gp->type == swift_type_dark_matter ||
gp->type == swift_type_dark_matter_background ||
gp->type == swift_type_neutrino) &&
gpart_is_active(gp, e)) {
const integertime_t ti_step = get_integer_timestep(gp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current, gp->time_bin) + ti_step / 2;
const integertime_t ti_end = ti_begin + ti_step / 2;
#ifdef SWIFT_DEBUG_CHECKS
if (ti_end != ti_current) error("Particle in wrong time-bin");
#endif
/* Time interval for this gravity half-kick */
const double dt_kick_grav = kick_get_grav_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
/* Finish the time-step with a second half-kick */
kick_gpart(gp, dt_kick_grav, ti_begin, ti_end, dt_kick_mesh_grav,
ti_begin_mesh, ti_end_mesh);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
if (gp->ti_drift != gp->ti_kick)
error("Error integrating g-part in time.");
#endif
/* Prepare the values to be drifted */
gravity_reset_predicted_values(gp);
}
}
/* Loop over the particles in this cell. */
for (int k = 0; k < scount; k++) {
/* Get a handle on the part. */
struct spart *restrict sp = &sparts[k];
/* If particle needs to be kicked */
if (spart_is_active(sp, e)) {
const integertime_t ti_step = get_integer_timestep(sp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current, sp->time_bin) + ti_step / 2;
const integertime_t ti_end = ti_begin + ti_step / 2;
#ifdef SWIFT_DEBUG_CHECKS
if (ti_end != ti_current) error("Particle in wrong time-bin");
#endif
/* Time interval for this gravity half-kick */
const double dt_kick_grav = kick_get_grav_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
/* Finish the time-step with a second half-kick */
kick_spart(sp, dt_kick_grav, ti_begin, ti_end, dt_kick_mesh_grav,
ti_begin_mesh, ti_end_mesh);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
if (sp->ti_drift != sp->ti_kick)
error("Error integrating s-part in time.");
#endif
/* Prepare the values to be drifted */
stars_reset_predicted_values(sp);
}
}
/* Loop over the sink particles in this cell. */
for (int k = 0; k < sink_count; k++) {
/* Get a handle on the part. */
struct sink *restrict sink = &sinks[k];
/* If particle needs to be kicked */
if (sink_is_active(sink, e)) {
const integertime_t ti_step = get_integer_timestep(sink->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current, sink->time_bin) + ti_step / 2;
const integertime_t ti_end = ti_begin + ti_step / 2;
#ifdef SWIFT_DEBUG_CHECKS
if (ti_end != ti_current) error("Particle in wrong time-bin");
#endif
/* Time interval for this gravity half-kick */
const double dt_kick_grav = kick_get_grav_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
/* Finish the time-step with a second half-kick */
kick_sink(sink, dt_kick_grav, ti_begin, ti_end, dt_kick_mesh_grav,
ti_begin_mesh, ti_end_mesh);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
if (sink->ti_drift != sink->ti_kick)
error("Error integrating sink-part in time.");
#endif
/* Prepare the values to be drifted */
sink_reset_predicted_values(sink);
}
}
/* Loop over the particles in this cell. */
for (int k = 0; k < bcount; k++) {
/* Get a handle on the part. */
struct bpart *restrict bp = &bparts[k];
/* If particle needs to be kicked */
if (bpart_is_active(bp, e)) {
const integertime_t ti_step = get_integer_timestep(bp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current, bp->time_bin) + ti_step / 2;
const integertime_t ti_end = ti_begin + ti_step / 2;
#ifdef SWIFT_DEBUG_CHECKS
if (ti_end != ti_current) error("Particle in wrong time-bin");
#endif
/* Time interval for this gravity half-kick */
const double dt_kick_grav = kick_get_grav_kick_dt(
ti_begin, ti_end, time_base, with_cosmology, cosmo);
/* Finish the time-step with a second half-kick */
kick_bpart(bp, dt_kick_grav, ti_begin, ti_end, dt_kick_mesh_grav,
ti_begin_mesh, ti_end_mesh);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
if (bp->ti_drift != bp->ti_kick)
error("Error integrating b-part in time.");
#endif
/* Prepare the values to be drifted */
black_holes_reset_predicted_values(bp);
}
}
}
if (timer) TIMER_TOC(timer_kick2);
}
/**
* @brief Computes the next time-step of all active particles in this cell
* and update the cell's statistics.
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
void runner_do_timestep(struct runner *r, struct cell *c, const int timer) {
const struct engine *e = r->e;
const integertime_t ti_current = e->ti_current;
const integertime_t ti_current_subcycle = e->ti_current_subcycle;
const struct cosmology *cosmo = e->cosmology;
const struct feedback_props *feedback_props = e->feedback_props;
const struct unit_system *us = e->internal_units;
const struct phys_const *phys_const = e->physical_constants;
const int with_cosmology = (e->policy & engine_policy_cosmology);
const int with_feedback = (e->policy & engine_policy_feedback);
const int with_rt = (e->policy & engine_policy_rt);
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
const int sink_count = c->sinks.count;
const int bcount = c->black_holes.count;
struct part *restrict parts = c->hydro.parts;
struct xpart *restrict xparts = c->hydro.xparts;
struct gpart *restrict gparts = c->grav.parts;
struct spart *restrict sparts = c->stars.parts;
struct sink *restrict sinks = c->sinks.parts;
struct bpart *restrict bparts = c->black_holes.parts;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e) &&
!cell_is_active_stars(c, e) && !cell_is_active_sinks(c, e) &&
!cell_is_active_black_holes(c, e)) {
/* Note: cell_is_rt_active is deliberately skipped. We only change
* the RT subcycling time steps when particles are hydro active. */
c->hydro.updated = 0;
c->grav.updated = 0;
c->stars.updated = 0;
c->sinks.updated = 0;
c->black_holes.updated = 0;
c->rt.updated = 0;
return;
}
int updated = 0, g_updated = 0, s_updated = 0, sink_updated = 0,
b_updated = 0;
integertime_t ti_hydro_end_min = max_nr_timesteps, ti_hydro_beg_max = 0;
integertime_t ti_rt_end_min = max_nr_timesteps, ti_rt_beg_max = 0;
integertime_t ti_rt_min_step_size = max_nr_timesteps;
integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_beg_max = 0;
integertime_t ti_stars_end_min = max_nr_timesteps, ti_stars_beg_max = 0;
integertime_t ti_sinks_end_min = max_nr_timesteps, ti_sinks_beg_max = 0;
integertime_t ti_black_holes_end_min = max_nr_timesteps,
ti_black_holes_beg_max = 0;
/* No children? */
if (!c->split) {
/* Loop over the particles in this cell. */
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
struct part *restrict p = &parts[k];
struct xpart *restrict xp = &xparts[k];
/* If particle needs updating */
if (part_is_active(p, e)) {
#ifdef SWIFT_DEBUG_CHECKS
/* Current end of time-step */
const integertime_t ti_end =
get_integer_time_end(ti_current, p->time_bin);
if (ti_end != ti_current)
error("Computing time-step of rogue particle.");
if (with_rt) {
const integertime_t ti_rt_end = get_integer_time_end(
ti_current_subcycle, p->rt_time_data.time_bin);
if (ti_rt_end != ti_current_subcycle)
error("Computing RT time-step of rogue particle");
}
#endif
/* Old time-step length in physical units */
const integertime_t ti_old_step = get_integer_timestep(p->time_bin);
double old_time_step_length;
if (with_cosmology) {
old_time_step_length = cosmology_get_delta_time(
e->cosmology, e->ti_current - ti_old_step, e->ti_current);
} else {
old_time_step_length = get_timestep(p->time_bin, e->time_base);
}
/* Get new time-step */
integertime_t ti_rt_new_step = get_part_rt_timestep(p, xp, e);
const integertime_t ti_new_step =
get_part_timestep(p, xp, e, ti_rt_new_step);
/* Enforce RT time-step size <= hydro step size. */
ti_rt_new_step = min(ti_new_step, ti_rt_new_step);
#ifdef SWIFT_RT_DEBUG_CHECKS
/* For the DEBUG RT scheme, this sets the RT time step to be
* (dt_hydro / max_nr_sub_cycles). For others, this does a proper
* debugging/consistency check. */
rt_debugging_check_timestep(p, &ti_rt_new_step, &ti_new_step,
e->max_nr_rt_subcycles, e->time_base);
if (ti_rt_new_step <= 0LL)
error("Got integer time step <= 0? %lld %lld",
get_part_rt_timestep(p, xp, e), ti_new_step);
#endif
/* Update particle */
p->time_bin = get_time_bin(ti_new_step);
if (p->gpart != NULL) p->gpart->time_bin = p->time_bin;
/* Update the tracers properties */
tracers_after_timestep_part(
p, xp, e->internal_units, e->physical_constants, with_cosmology,
e->cosmology, e->hydro_properties, e->cooling_func, e->time,
old_time_step_length, e->snapshot_recording_triggers_started_part);
/* Number of updated particles */
updated++;
if (p->gpart != NULL) g_updated++;
/* What is the next sync-point ? */
ti_hydro_end_min = min(ti_current + ti_new_step, ti_hydro_end_min);
/* What is the next starting point for this cell ? */
ti_hydro_beg_max = max(ti_current, ti_hydro_beg_max);
if (p->gpart != NULL) {
/* What is the next sync-point ? */
ti_gravity_end_min =
min(ti_current + ti_new_step, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
}
/* Same for RT */
if (with_rt) {
p->rt_time_data.time_bin = get_time_bin(ti_rt_new_step);
ti_rt_end_min =
min(ti_current_subcycle + ti_rt_new_step, ti_rt_end_min);
ti_rt_beg_max =
max(ti_current_subcycle + ti_rt_new_step, ti_rt_beg_max);
ti_rt_min_step_size = min(ti_rt_min_step_size, ti_rt_new_step);
}
}
else { /* part is inactive */
if (!part_is_inhibited(p, e)) {
const integertime_t ti_end =
get_integer_time_end(ti_current, p->time_bin);
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, p->time_bin);
/* What is the next sync-point ? */
ti_hydro_end_min = min(ti_end, ti_hydro_end_min);
/* What is the next starting point for this cell ? */
ti_hydro_beg_max = max(ti_beg, ti_hydro_beg_max);
/* Same for RT. */
if (with_rt) {
/* Here we assume that the particle is inactive, which is true for
* hydro, but not necessarily for a RT subcycle. RT time steps are
* only changed while the particle is hydro active. This allows to
* end up with results ti_rt_end == ti_current_subcyle, so we need
* to pretend we're past ti_current_subcycle already. */
integertime_t ti_rt_end = get_integer_time_end(
ti_current_subcycle + 1, p->rt_time_data.time_bin);
const integertime_t ti_rt_beg = get_integer_time_begin(
ti_current_subcycle + 1, p->rt_time_data.time_bin);
ti_rt_end_min = min(ti_rt_end, ti_rt_end_min);
ti_rt_beg_max = max(ti_rt_beg, ti_rt_beg_max);
integertime_t ti_rt_step =
get_integer_timestep(p->rt_time_data.time_bin);
ti_rt_min_step_size = min(ti_rt_min_step_size, ti_rt_step);
}
if (p->gpart != NULL) {
/* What is the next sync-point ? */
ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
}
}
}
}
/* Loop over the g-particles in this cell. */
for (int k = 0; k < gcount; k++) {
/* Get a handle on the part. */
struct gpart *restrict gp = &gparts[k];
/* If the g-particle has no counterpart */
if (gp->type == swift_type_dark_matter ||
gp->type == swift_type_dark_matter_background ||
gp->type == swift_type_neutrino) {
/* need to be updated ? */
if (gpart_is_active(gp, e)) {
#ifdef SWIFT_DEBUG_CHECKS
/* Current end of time-step */
const integertime_t ti_end =
get_integer_time_end(ti_current, gp->time_bin);
if (ti_end != ti_current)
error("Computing time-step of rogue particle.");
#endif
/* Get new time-step */
const integertime_t ti_new_step = get_gpart_timestep(gp, e);
/* Update particle */
gp->time_bin = get_time_bin(ti_new_step);
/* Number of updated g-particles */
g_updated++;
/* What is the next sync-point ? */
ti_gravity_end_min =
min(ti_current + ti_new_step, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
} else { /* gpart is inactive */
if (!gpart_is_inhibited(gp, e)) {
const integertime_t ti_end =
get_integer_time_end(ti_current, gp->time_bin);
/* What is the next sync-point ? */
ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, gp->time_bin);
/* What is the next starting point for this cell ? */
ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
}
}
}
}
/* Loop over the star particles in this cell. */
for (int k = 0; k < scount; k++) {
/* Get a handle on the part. */
struct spart *restrict sp = &sparts[k];
/* need to be updated ? */
if (spart_is_active(sp, e)) {
#ifdef SWIFT_DEBUG_CHECKS
/* Current end of time-step */
const integertime_t ti_end =
get_integer_time_end(ti_current, sp->time_bin);
if (ti_end != ti_current)
error("Computing time-step of rogue particle.");
#endif
/* Old time-step length in physical units */
const integertime_t ti_old_step = get_integer_timestep(sp->time_bin);
double old_time_step_length;
if (with_cosmology) {
old_time_step_length = cosmology_get_delta_time(
e->cosmology, e->ti_current - ti_old_step, e->ti_current);
} else {
old_time_step_length = get_timestep(sp->time_bin, e->time_base);
}
/* Get new time-step */
const integertime_t ti_new_step = get_spart_timestep(sp, e);
/* Update particle */
sp->time_bin = get_time_bin(ti_new_step);
sp->gpart->time_bin = get_time_bin(ti_new_step);
/* Update the tracers properties */
tracers_after_timestep_spart(
sp, e->internal_units, e->physical_constants, with_cosmology,
e->cosmology, old_time_step_length,
e->snapshot_recording_triggers_started_spart);
/* Update feedback related counters */
if (with_feedback) {
feedback_will_do_feedback(sp, feedback_props, with_cosmology, cosmo,
e->time, us, phys_const, e->ti_current,
e->time_base);
}
/* Number of updated s-particles */
s_updated++;
g_updated++;
ti_stars_end_min = min(ti_current + ti_new_step, ti_stars_end_min);
ti_gravity_end_min = min(ti_current + ti_new_step, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_stars_beg_max = max(ti_current, ti_stars_beg_max);
ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
/* star particle is inactive but not inhibited */
} else {
if (!spart_is_inhibited(sp, e)) {
const integertime_t ti_end =
get_integer_time_end(ti_current, sp->time_bin);
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, sp->time_bin);
ti_stars_end_min = min(ti_end, ti_stars_end_min);
ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_stars_beg_max = max(ti_beg, ti_stars_beg_max);
ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
}
}
}
/* Loop over the sink particles in this cell. */
for (int k = 0; k < sink_count; k++) {
/* Get a handle on the part. */
struct sink *restrict sink = &sinks[k];
/* need to be updated ? */
if (sink_is_active(sink, e)) {
#ifdef SWIFT_DEBUG_CHECKS
/* Current end of time-step */
const integertime_t ti_end =
get_integer_time_end(ti_current, sink->time_bin);
if (ti_end != ti_current)
error("Computing time-step of rogue particle.");
#endif
/* Get new time-step */
const integertime_t ti_new_step = get_sink_timestep(sink, e);
/* Update particle */
sink->time_bin = get_time_bin(ti_new_step);
sink->gpart->time_bin = get_time_bin(ti_new_step);
/* Number of updated sink-particles */
sink_updated++;
g_updated++;
ti_sinks_end_min = min(ti_current + ti_new_step, ti_sinks_end_min);
ti_gravity_end_min = min(ti_current + ti_new_step, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_sinks_beg_max = max(ti_current, ti_sinks_beg_max);
ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
/* sink particle is inactive but not inhibited */
} else {
if (!sink_is_inhibited(sink, e)) {
const integertime_t ti_end =
get_integer_time_end(ti_current, sink->time_bin);
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, sink->time_bin);
ti_sinks_end_min = min(ti_end, ti_sinks_end_min);
ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_sinks_beg_max = max(ti_beg, ti_sinks_beg_max);
ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
}
}
}
/* Loop over the black hole particles in this cell. */
for (int k = 0; k < bcount; k++) {
/* Get a handle on the part. */
struct bpart *restrict bp = &bparts[k];
/* need to be updated ? */
if (bpart_is_active(bp, e)) {
#ifdef SWIFT_DEBUG_CHECKS
/* Current end of time-step */
const integertime_t ti_end =
get_integer_time_end(ti_current, bp->time_bin);
if (ti_end != ti_current)
error("Computing time-step of rogue particle.");
#endif
/* Old time-step length in physical units */
const integertime_t ti_old_step = get_integer_timestep(bp->time_bin);
double old_time_step_length;
if (with_cosmology) {
old_time_step_length = cosmology_get_delta_time(
e->cosmology, e->ti_current - ti_old_step, e->ti_current);
} else {
old_time_step_length = get_timestep(bp->time_bin, e->time_base);
}
/* Get new time-step */
const integertime_t ti_new_step = get_bpart_timestep(bp, e);
/* Update particle */
bp->time_bin = get_time_bin(ti_new_step);
bp->gpart->time_bin = get_time_bin(ti_new_step);
/* Update the tracers properties */
tracers_after_timestep_bpart(
bp, e->internal_units, e->physical_constants, with_cosmology,
e->cosmology, old_time_step_length,
e->snapshot_recording_triggers_started_bpart);
/* Number of updated s-particles */
b_updated++;
g_updated++;
ti_black_holes_end_min =
min(ti_current + ti_new_step, ti_black_holes_end_min);
ti_gravity_end_min = min(ti_current + ti_new_step, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_black_holes_beg_max = max(ti_current, ti_black_holes_beg_max);
ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
/* star particle is inactive but not inhibited */
} else {
if (!bpart_is_inhibited(bp, e)) {
const integertime_t ti_end =
get_integer_time_end(ti_current, bp->time_bin);
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, bp->time_bin);
ti_black_holes_end_min = min(ti_end, ti_black_holes_end_min);
ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_black_holes_beg_max = max(ti_beg, ti_black_holes_beg_max);
ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
}
}
}
} else {
/* Loop over the progeny. */
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
struct cell *restrict cp = c->progeny[k];
/* Recurse */
runner_do_timestep(r, cp, /*timer=*/0);
/* And aggregate */
updated += cp->hydro.updated;
g_updated += cp->grav.updated;
sink_updated += cp->sinks.updated;
s_updated += cp->stars.updated;
b_updated += cp->black_holes.updated;
ti_hydro_end_min = min(cp->hydro.ti_end_min, ti_hydro_end_min);
ti_hydro_beg_max = max(cp->hydro.ti_beg_max, ti_hydro_beg_max);
ti_rt_end_min = min(cp->rt.ti_rt_end_min, ti_rt_end_min);
ti_rt_beg_max = max(cp->rt.ti_rt_beg_max, ti_rt_beg_max);
ti_rt_min_step_size =
min(cp->rt.ti_rt_min_step_size, ti_rt_min_step_size);
ti_gravity_end_min = min(cp->grav.ti_end_min, ti_gravity_end_min);
ti_gravity_beg_max = max(cp->grav.ti_beg_max, ti_gravity_beg_max);
ti_stars_end_min = min(cp->stars.ti_end_min, ti_stars_end_min);
ti_stars_beg_max = max(cp->stars.ti_beg_max, ti_stars_beg_max);
ti_sinks_end_min = min(cp->sinks.ti_end_min, ti_sinks_end_min);
ti_sinks_beg_max = max(cp->sinks.ti_beg_max, ti_sinks_beg_max);
ti_black_holes_end_min =
min(cp->black_holes.ti_end_min, ti_black_holes_end_min);
ti_black_holes_beg_max =
max(cp->grav.ti_beg_max, ti_black_holes_beg_max);
}
}
}
/* Flag something may have changed */
if (c->top == c) space_mark_cell_as_updated(r->e->s, c);
/* Store the values. */
c->hydro.updated = updated;
c->grav.updated = g_updated;
c->stars.updated = s_updated;
c->sinks.updated = sink_updated;
c->black_holes.updated = b_updated;
/* We don't count the RT updates here because the
* timestep tasks aren't active during sub-cycles.
* We do that in rt_advanced_cell_time instead. */
c->hydro.ti_end_min = ti_hydro_end_min;
c->hydro.ti_beg_max = ti_hydro_beg_max;
c->rt.ti_rt_end_min = ti_rt_end_min;
c->rt.ti_rt_beg_max = ti_rt_beg_max;
if (cell_is_starting_hydro(c, e)) {
/* We only change the RT time steps when the cell is also hydro active.
* Without this check here, ti_rt_min_step_size = max_nr_steps... */
c->rt.ti_rt_min_step_size = ti_rt_min_step_size;
}
c->grav.ti_end_min = ti_gravity_end_min;
c->grav.ti_beg_max = ti_gravity_beg_max;
c->stars.ti_end_min = ti_stars_end_min;
c->stars.ti_beg_max = ti_stars_beg_max;
c->sinks.ti_end_min = ti_sinks_end_min;
c->sinks.ti_beg_max = ti_sinks_beg_max;
c->black_holes.ti_end_min = ti_black_holes_end_min;
c->black_holes.ti_beg_max = ti_black_holes_beg_max;
#ifdef SWIFT_DEBUG_CHECKS
if (c->hydro.ti_end_min == e->ti_current &&
c->hydro.ti_end_min < max_nr_timesteps)
error("End of next hydro step is current time!");
if (c->grav.ti_end_min == e->ti_current &&
c->grav.ti_end_min < max_nr_timesteps)
error("End of next gravity step is current time!");
if (c->stars.ti_end_min == e->ti_current &&
c->stars.ti_end_min < max_nr_timesteps)
error("End of next stars step is current time!");
if (c->sinks.ti_end_min == e->ti_current &&
c->sinks.ti_end_min < max_nr_timesteps)
error("End of next sinks step is current time!");
if (c->black_holes.ti_end_min == e->ti_current &&
c->black_holes.ti_end_min < max_nr_timesteps)
error("End of next black holes step is current time!");
/* Contrary to sinks, stars, bhs etc, we may have "rt particles"
* without running with RT. So additional if (with_rt) check is
* needed here. */
if (with_rt && (c->rt.ti_rt_end_min == e->ti_current &&
c->rt.ti_rt_end_min < max_nr_timesteps))
error("Cell %lld End of next RT step is current time!", c->cellID);
#endif
if (timer) TIMER_TOC(timer_timestep);
}
/**
* @brief Recursively collect the end-of-timestep information from the top-level
* to the super level.
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
void runner_do_timestep_collect(struct runner *r, struct cell *c,
const int timer) {
/* Early stop if we are at the super level.
* The time-step task would have set things at this level already */
if (c->super == c) return;
/* Counters for the different quantities. */
size_t h_updated = 0;
size_t g_updated = 0;
size_t s_updated = 0;
size_t b_updated = 0;
size_t si_updated = 0;
size_t rt_updated = 0;
integertime_t ti_hydro_end_min = max_nr_timesteps, ti_hydro_beg_max = 0;
integertime_t ti_rt_end_min = max_nr_timesteps, ti_rt_beg_max = 0;
integertime_t ti_grav_end_min = max_nr_timesteps, ti_grav_beg_max = 0;
integertime_t ti_stars_end_min = max_nr_timesteps, ti_stars_beg_max = 0;
integertime_t ti_black_holes_end_min = max_nr_timesteps,
ti_black_holes_beg_max = 0;
integertime_t ti_sinks_end_min = max_nr_timesteps, ti_sinks_beg_max = 0;
/* Collect the values from the progeny. */
for (int k = 0; k < 8; k++) {
struct cell *cp = c->progeny[k];
if (cp != NULL) {
/* Recurse */
runner_do_timestep_collect(r, cp, 0);
/* And update */
ti_hydro_end_min = min(ti_hydro_end_min, cp->hydro.ti_end_min);
ti_hydro_beg_max = max(ti_hydro_beg_max, cp->hydro.ti_beg_max);
ti_rt_end_min = min(cp->rt.ti_rt_end_min, ti_rt_end_min);
ti_rt_beg_max = max(cp->rt.ti_rt_beg_max, ti_rt_beg_max);
ti_grav_end_min = min(ti_grav_end_min, cp->grav.ti_end_min);
ti_grav_beg_max = max(ti_grav_beg_max, cp->grav.ti_beg_max);
ti_stars_end_min = min(ti_stars_end_min, cp->stars.ti_end_min);
ti_stars_beg_max = max(ti_stars_beg_max, cp->stars.ti_beg_max);
ti_black_holes_end_min =
min(ti_black_holes_end_min, cp->black_holes.ti_end_min);
ti_black_holes_beg_max =
max(ti_black_holes_beg_max, cp->black_holes.ti_beg_max);
ti_sinks_end_min = min(ti_sinks_end_min, cp->sinks.ti_end_min);
ti_sinks_beg_max = max(ti_sinks_beg_max, cp->sinks.ti_beg_max);
h_updated += cp->hydro.updated;
g_updated += cp->grav.updated;
s_updated += cp->stars.updated;
b_updated += cp->black_holes.updated;
si_updated += cp->sinks.updated;
rt_updated += cp->rt.updated;
/* Collected, so clear for next time. */
cp->hydro.updated = 0;
cp->grav.updated = 0;
cp->stars.updated = 0;
cp->black_holes.updated = 0;
cp->sinks.updated = 0;
cp->rt.updated = 0;
}
}
/* Flag something may have changed */
if (c->top == c) space_mark_cell_as_updated(r->e->s, c);
/* Store the collected values in the cell. */
c->hydro.ti_end_min = ti_hydro_end_min;
c->hydro.ti_beg_max = ti_hydro_beg_max;
c->rt.ti_rt_end_min = ti_rt_end_min;
c->rt.ti_rt_beg_max = ti_rt_beg_max;
c->grav.ti_end_min = ti_grav_end_min;
c->grav.ti_beg_max = ti_grav_beg_max;
c->stars.ti_end_min = ti_stars_end_min;
c->stars.ti_beg_max = ti_stars_beg_max;
c->black_holes.ti_end_min = ti_black_holes_end_min;
c->black_holes.ti_beg_max = ti_black_holes_beg_max;
c->sinks.ti_end_min = ti_sinks_end_min;
c->sinks.ti_beg_max = ti_sinks_beg_max;
c->hydro.updated = h_updated;
c->grav.updated = g_updated;
c->stars.updated = s_updated;
c->black_holes.updated = b_updated;
c->sinks.updated = si_updated;
c->rt.updated = rt_updated;
}
/**
* @brief Apply the time-step limiter to all awaken particles in a cell
* hierarchy.
*
* @param r The task #runner.
* @param c The #cell.
* @param force Limit the particles irrespective of the #cell flags.
* @param timer Are we timing this ?
*/
void runner_do_limiter(struct runner *r, struct cell *c, int force,
const int timer) {
const struct engine *e = r->e;
const int count = c->hydro.count;
struct part *restrict parts = c->hydro.parts;
struct xpart *restrict xparts = c->hydro.xparts;
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
/* Check that we only limit local cells. */
if (c->nodeID != engine_rank) error("Limiting dt of a foreign cell is nope.");
#endif
integertime_t ti_hydro_end_min = max_nr_timesteps, ti_hydro_beg_max = 0;
integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_beg_max = 0;
/* Limit irrespective of cell flags? */
force = (force || cell_get_flag(c, cell_flag_do_hydro_limiter));
/* Early abort? */
if (c->hydro.count == 0) {
/* Clear the limiter flags. */
cell_clear_flag(
c, cell_flag_do_hydro_limiter | cell_flag_do_hydro_sub_limiter);
return;
}
/* Loop over the progeny ? */
if (c->split && (force || cell_get_flag(c, cell_flag_do_hydro_sub_limiter))) {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
struct cell *restrict cp = c->progeny[k];
/* Recurse */
runner_do_limiter(r, cp, force, /*timer=*/0);
/* And aggregate */
ti_hydro_end_min = min(cp->hydro.ti_end_min, ti_hydro_end_min);
ti_hydro_beg_max = max(cp->hydro.ti_beg_max, ti_hydro_beg_max);
ti_gravity_end_min = min(cp->grav.ti_end_min, ti_gravity_end_min);
ti_gravity_beg_max = max(cp->grav.ti_beg_max, ti_gravity_beg_max);
}
}
/* Flag something may have changed */
if (c->top == c) space_mark_cell_as_updated(r->e->s, c);
/* Store the updated values */
c->hydro.ti_end_min = min(c->hydro.ti_end_min, ti_hydro_end_min);
c->hydro.ti_beg_max = max(c->hydro.ti_beg_max, ti_hydro_beg_max);
c->grav.ti_end_min = min(c->grav.ti_end_min, ti_gravity_end_min);
c->grav.ti_beg_max = max(c->grav.ti_beg_max, ti_gravity_beg_max);
} else if (!c->split && force) {
ti_hydro_end_min = c->hydro.ti_end_min;
ti_hydro_beg_max = c->hydro.ti_beg_max;
ti_gravity_end_min = c->grav.ti_end_min;
ti_gravity_beg_max = c->grav.ti_beg_max;
/* Loop over the gas particles in this cell. */
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
struct part *restrict p = &parts[k];
struct xpart *restrict xp = &xparts[k];
#ifdef SWIFT_HYDRO_DENSITY_CHECKS
/* Finish the limiter loop by adding a (fake) self-contribution */
p->limiter_data.N_limiter++;
const float h_inv_dim = pow_dimension(1. / p->h); /* 1/h^d */
p->limiter_data.n_limiter += kernel_root;
p->limiter_data.n_limiter *= h_inv_dim;
#endif
/* Avoid inhibited particles */
if (part_is_inhibited(p, e)) continue;
/* Bip, bip, bip... wake-up time */
if (p->limiter_data.wakeup != time_bin_not_awake) {
if (!part_is_active(p, e) && p->limiter_data.to_be_synchronized) {
warning(
"Not limiting particle with id %lld because it needs to be "
"synced.",
p->id);
continue;
}
// message("Limiting particle %lld in cell %lld", p->id, c->cellID);
/* Apply the limiter and get the new end of time-step */
const integertime_t ti_end_new = timestep_limit_part(p, xp, e);
const timebin_t new_bin = p->time_bin;
const integertime_t ti_beg_new =
ti_end_new - get_integer_timestep(new_bin);
/* Mark this particle has not needing synchronization */
p->limiter_data.to_be_synchronized = 0;
#ifdef SWIFT_HYDRO_DENSITY_CHECKS
p->limited_part = 1;
#endif
/* What is the next sync-point ? */
ti_hydro_end_min = min(ti_end_new, ti_hydro_end_min);
/* What is the next starting point for this cell ? */
ti_hydro_beg_max = max(ti_beg_new, ti_hydro_beg_max);
/* Also limit the gpart counter-part */
if (p->gpart != NULL) {
/* Register the time-bin */
p->gpart->time_bin = p->time_bin;
/* What is the next sync-point ? */
ti_gravity_end_min = min(ti_end_new, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_gravity_beg_max = max(ti_beg_new, ti_gravity_beg_max);
}
}
}
/* Flag something may have changed */
if (c->top == c) space_mark_cell_as_updated(r->e->s, c);
/* Store the updated values */
c->hydro.ti_end_min = min(c->hydro.ti_end_min, ti_hydro_end_min);
c->hydro.ti_beg_max = max(c->hydro.ti_beg_max, ti_hydro_beg_max);
c->grav.ti_end_min = min(c->grav.ti_end_min, ti_gravity_end_min);
c->grav.ti_beg_max = max(c->grav.ti_beg_max, ti_gravity_beg_max);
}
/* Clear the limiter flags. */
cell_clear_flag(c,
cell_flag_do_hydro_limiter | cell_flag_do_hydro_sub_limiter);
if (timer) TIMER_TOC(timer_do_limiter);
}
/**
* @brief Apply the time-step synchronization proceduere to all flagged
* particles in a cell hierarchy.
*
* @param r The task #runner.
* @param c The #cell.
* @param force Limit the particles irrespective of the #cell flags.
* @param timer Are we timing this ?
*/
void runner_do_sync(struct runner *r, struct cell *c, int force,
const int timer) {
const struct engine *e = r->e;
const integertime_t ti_current = e->ti_current;
const struct cosmology *cosmo = e->cosmology;
const int with_cosmology = (e->policy & engine_policy_cosmology);
const int count = c->hydro.count;
struct part *restrict parts = c->hydro.parts;
struct xpart *restrict xparts = c->hydro.xparts;
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
/* Check that we only sync local cells. */
if (c->nodeID != engine_rank) error("Syncing of a foreign cell is nope.");
#endif
integertime_t ti_hydro_end_min = max_nr_timesteps, ti_hydro_beg_max = 0;
integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_beg_max = 0;
/* Limit irrespective of cell flags? */
force = (force || cell_get_flag(c, cell_flag_do_hydro_sync));
/* Early abort? */
if (c->hydro.count == 0) {
/* Clear the sync flags. */
cell_clear_flag(c, cell_flag_do_hydro_sync | cell_flag_do_hydro_sub_sync);
return;
}
/* Loop over the progeny ? */
if (c->split && (force || cell_get_flag(c, cell_flag_do_hydro_sub_sync))) {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
struct cell *restrict cp = c->progeny[k];
/* Recurse */
runner_do_sync(r, cp, force, /*timer=*/0);
/* And aggregate */
ti_hydro_end_min = min(cp->hydro.ti_end_min, ti_hydro_end_min);
ti_hydro_beg_max = max(cp->hydro.ti_beg_max, ti_hydro_beg_max);
ti_gravity_end_min = min(cp->grav.ti_end_min, ti_gravity_end_min);
ti_gravity_beg_max = max(cp->grav.ti_beg_max, ti_gravity_beg_max);
}
}
/* Flag something may have changed */
if (c->top == c) space_mark_cell_as_updated(r->e->s, c);
/* Store the updated values */
c->hydro.ti_end_min = min(c->hydro.ti_end_min, ti_hydro_end_min);
c->hydro.ti_beg_max = max(c->hydro.ti_beg_max, ti_hydro_beg_max);
c->grav.ti_end_min = min(c->grav.ti_end_min, ti_gravity_end_min);
c->grav.ti_beg_max = max(c->grav.ti_beg_max, ti_gravity_beg_max);
} else if (!c->split && force) {
ti_hydro_end_min = c->hydro.ti_end_min;
ti_hydro_beg_max = c->hydro.ti_beg_max;
ti_gravity_end_min = c->grav.ti_end_min;
ti_gravity_beg_max = c->grav.ti_beg_max;
/* Loop over the gas particles in this cell. */
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
struct part *restrict p = &parts[k];
struct xpart *restrict xp = &xparts[k];
/* Avoid inhibited particles */
if (part_is_inhibited(p, e)) continue;
/* If the particle is active no need to sync it */
if (part_is_active(p, e) && p->limiter_data.to_be_synchronized) {
p->limiter_data.to_be_synchronized = 0;
}
if (p->limiter_data.to_be_synchronized) {
/* Finish this particle's time-step */
timestep_process_sync_part(p, xp, e, cosmo);
/* Note that at this moment the new RT time step is only used to
* limit the hydro time step here. */
integertime_t ti_rt_new_step = get_part_rt_timestep(p, xp, e);
/* Get new time-step */
integertime_t ti_new_step = get_part_timestep(p, xp, e, ti_rt_new_step);
timebin_t new_time_bin = get_time_bin(ti_new_step);
/* Enforce RT time-step size <= hydro step size. */
/* On the commented out line below: We should be doing this once we
* correctly add RT to this part of the code. */
/* ti_rt_new_step = min(ti_new_step, ti_rt_new_step); */
/* Apply the limiter if necessary */
if (p->limiter_data.wakeup != time_bin_not_awake) {
new_time_bin = min(new_time_bin, -p->limiter_data.wakeup + 2);
p->limiter_data.wakeup = time_bin_not_awake;
}
/* Limit the time-bin to what is allowed in this step */
new_time_bin = min(new_time_bin, e->max_active_bin);
ti_new_step = get_integer_timestep(new_time_bin);
/* Time-step length in physical units */
// MATTHIEU: TODO: think about this one!
double time_step_length;
if (with_cosmology) {
time_step_length = cosmology_get_delta_time(
e->cosmology, e->ti_current, e->ti_current + ti_new_step);
} else {
time_step_length = get_timestep(new_time_bin, e->time_base);
}
/* Update particle */
p->time_bin = new_time_bin;
if (p->gpart != NULL) p->gpart->time_bin = new_time_bin;
/* Update the tracers properties */
tracers_after_timestep_part(
p, xp, e->internal_units, e->physical_constants, with_cosmology,
e->cosmology, e->hydro_properties, e->cooling_func, e->time,
0 * time_step_length, e->snapshot_recording_triggers_started_part);
#ifdef SWIFT_HYDRO_DENSITY_CHECKS
p->limited_part = 1;
#endif
/* What is the next sync-point ? */
ti_hydro_end_min = min(ti_current + ti_new_step, ti_hydro_end_min);
/* What is the next starting point for this cell ? */
ti_hydro_beg_max = max(ti_current, ti_hydro_beg_max);
/* Also limit the gpart counter-part */
if (p->gpart != NULL) {
/* Register the time-bin */
p->gpart->time_bin = p->time_bin;
/* What is the next sync-point ? */
ti_gravity_end_min =
min(ti_current + ti_new_step, ti_gravity_end_min);
/* What is the next starting point for this cell ? */
ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
}
}
}
/* Flag something may have changed */
if (c->top == c) space_mark_cell_as_updated(r->e->s, c);
/* Store the updated values */
c->hydro.ti_end_min = min(c->hydro.ti_end_min, ti_hydro_end_min);
c->hydro.ti_beg_max = max(c->hydro.ti_beg_max, ti_hydro_beg_max);
c->grav.ti_end_min = min(c->grav.ti_end_min, ti_gravity_end_min);
c->grav.ti_beg_max = max(c->grav.ti_beg_max, ti_gravity_beg_max);
}
/* Clear the sync flags. */
cell_clear_flag(c, cell_flag_do_hydro_sync | cell_flag_do_hydro_sub_sync);
if (timer) TIMER_TOC(timer_do_sync);
}
/**
* @brief Update the cell's t_rt_end_min so that the sub-cycling can proceed
* with correct cell times. (During sub-cycles, the regular timestep and
* timestep_collect tasks do not run. This replaces the collection of cell
* times of timestep tasks during sub-cycles. )
*
* @param r The #runner thread.
* @param c The #cell.
* @param timer Are we timing this ?
*/
void runner_do_rt_advance_cell_time(struct runner *r, struct cell *c,
int timer) {
struct engine *e = r->e;
const int count = c->hydro.count;
/* Reset update count regardless whether cell is active or not */
c->rt.updated = 0;
#ifdef SWIFT_RT_DEBUG_CHECKS
if (c->super == c) c->rt.advanced_time = 1;
#endif
/* Anything to do here? */
if (count == 0) return;
if (!cell_is_rt_active(c, e)) return;
TIMER_TIC;
int rt_updated = 0;
if (c->split) {
for (int k = 0; k < 8; k++) {
struct cell *cp = c->progeny[k];
if (cp != NULL) {
runner_do_rt_advance_cell_time(r, cp, 0);
rt_updated += cp->rt.updated;
}
}
} else {
/* Do some debugging stuff on active particles before setting the cell time.
* This test is not reliable on foreign cells. After a rebuild, we may end
* up with an active foreign cell which was not updated in this step because
* it had no active neighbouring cells, and its particle data may be random
* junk. */
if (c->nodeID == engine_rank) {
struct part *restrict parts = c->hydro.parts;
/* Loop over the gas particles in this cell. */
for (int i = 0; i < count; i++) {
/* Get a handle on the part. */
struct part *restrict p = &parts[i];
/* Skip inhibited parts */
if (part_is_inhibited(p, e)) continue;
/* Skip inactive parts */
if (!part_is_rt_active(p, e)) continue;
#ifdef SWIFT_RT_DEBUG_CHECKS
/* Run checks. */
rt_debug_sequence_check(p, 5, __func__);
/* Mark that the subcycling has happened */
rt_debugging_count_subcycle(p);
#endif
rt_updated++;
}
}
}
c->rt.updated = rt_updated;
/* Note: c->rt.ti_rt_min_step_size may be greater than
* c->super->rt.ti_rt_min_step_size. This is expected behaviour.
* We only update the cell's own time after it's been active. */
c->rt.ti_rt_end_min += c->rt.ti_rt_min_step_size;
if (timer) TIMER_TOC(timer_do_rt_advance_cell_time);
}
/**
* @brief Recursively collect the end-of-timestep information from the top-level
* to the super level for the RT sub-cycling. (During sub-cycles, the regular
* timestep and timestep_collect tasks do not run. This replaces the
* timestep_collect task.)
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
void runner_do_collect_rt_times(struct runner *r, struct cell *c,
const int timer) {
TIMER_TIC;
const struct engine *e = r->e;
size_t rt_updated = 0;
if (e->ti_current == e->ti_current_subcycle)
error("called collect_rt_times during a main step");
/* Early stop if we are at the super level.
* The time-step/rt_advance_cell_time tasks would have set things at
* this level already. */
if (c->super == c) {
#ifdef SWIFT_RT_DEBUG_CHECKS
/* Do a check before the early exit.
* rt_advanced_cell_time should be called exactly once before
* collect times. Except on the first subcycle, because the
* collect_rt_times task shouldn't be called in the main steps.
* In that case, it should be exactly 2.
* This is only valid if the cell has been active this step.
* Otherwise, rt_advance_cell_time will not have run, yet the
* rt_collect_cell_times call may still be executed if the top
* level cell is above the super level cell. */
if (!cell_is_rt_active(c, e)) return;
if (e->ti_current_subcycle - c->rt.ti_rt_end_min == e->ti_current) {
/* This is the first subcycle */
if (c->rt.advanced_time != 2)
error("Called cell with wrong advanced_time counter. Expect=2, got=%d",
c->rt.advanced_time);
} else {
if (c->rt.advanced_time != 1)
error("Called cell with wrong advanced_time counter. Expect=1, got=%d",
c->rt.advanced_time);
}
c->rt.advanced_time = 0;
#endif
return;
}
integertime_t ti_rt_end_min = max_nr_timesteps, ti_rt_beg_max = 0;
/* Collect the values from the progeny. */
for (int k = 0; k < 8; k++) {
struct cell *cp = c->progeny[k];
if (cp != NULL) {
/* Recurse */
runner_do_collect_rt_times(r, cp, 0);
/* And update */
ti_rt_end_min = min(cp->rt.ti_rt_end_min, ti_rt_end_min);
ti_rt_beg_max = max(cp->rt.ti_rt_beg_max, ti_rt_beg_max);
/* Collected, so clear for next time. */
rt_updated += cp->rt.updated;
cp->rt.updated = 0;
}
}
/* Store the collected values in the cell. */
c->rt.ti_rt_end_min = ti_rt_end_min;
c->rt.ti_rt_beg_max = ti_rt_beg_max;
c->rt.updated = rt_updated;
if (timer) TIMER_TOC(timer_do_rt_collect_times);
}