Commit 2d165a9d authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'removing_particles' into 'master'

Add functions to permanently remove particles from a simulation

Closes #470

See merge request !632
parents ac37e89a 4427cce3
......@@ -54,7 +54,7 @@ infile = args.input
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
"init_grav", "init_grav_out", "ghost_in", "ghost", "ghost_out", "extra_ghost", "drift_part", "drift_gpart",
"end_force", "kick1", "kick2", "timestep", "send", "recv", "grav_long_range", "grav_mm", "grav_down_in",
"grav_down", "grav_mesh", "cooling", "sourceterms",
"grav_down", "grav_mesh", "cooling", "star_formation", "sourceterms",
"stars_ghost_in", "stars_ghost", "stars_ghost_out",
"count"]
......
......@@ -923,6 +923,9 @@ int main(int argc, char *argv[]) {
if (with_structure_finding)
engine_policies |= engine_policy_structure_finding;
// MATTHIEU: Temporary star formation law
// engine_policies |= engine_policy_star_formation;
/* Initialize the engine with the space and policies. */
if (myrank == 0) clocks_gettime(&tic);
engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
......
......@@ -112,7 +112,7 @@ pl.rcParams.update(PLOT_PARAMS)
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
"init_grav", "init_grav_out", "ghost_in", "ghost", "ghost_out", "extra_ghost", "drift_part", "drift_gpart",
"end_force", "kick1", "kick2", "timestep", "send", "recv", "grav_long_range", "grav_mm", "grav_down_in",
"grav_down", "grav_mesh", "cooling", "sourceterms",
"grav_down", "grav_mesh", "cooling", "star_formation", "sourceterms",
"stars_ghost_in", "stars_ghost", "stars_ghost_out",
"count"]
......
......@@ -278,6 +278,42 @@ __attribute__((always_inline)) INLINE static int spart_is_active(
return (spart_bin <= max_active_bin);
}
/**
* @brief Has this particle been inhibited?
*
* @param p The #part.
* @param e The #engine containing information about the current time.
* @return 1 if the #part is inhibited, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int part_is_inhibited(
const struct part *p, const struct engine *e) {
return p->time_bin == time_bin_inhibited;
}
/**
* @brief Has this gravity particle been inhibited?
*
* @param gp The #gpart.
* @param e The #engine containing information about the current time.
* @return 1 if the #part is inhibited, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int gpart_is_inhibited(
const struct gpart *gp, const struct engine *e) {
return gp->time_bin == time_bin_inhibited;
}
/**
* @brief Has this star particle been inhibited?
*
* @param sp The #spart.
* @param e The #engine containing information about the current time.
* @return 1 if the #part is inhibited, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int spart_is_inhibited(
const struct spart *sp, const struct engine *e) {
return sp->time_bin == time_bin_inhibited;
}
/* Are cells / particles active for kick1 tasks ? */
/**
......
......@@ -1246,7 +1246,8 @@ void cell_check_part_drift_point(struct cell *c, void *data) {
c->hydro.ti_old_part, ti_drift);
for (int i = 0; i < c->hydro.count; ++i)
if (c->hydro.parts[i].ti_drift != ti_drift)
if (c->hydro.parts[i].ti_drift != ti_drift &&
c->hydro.parts[i].time_bin != time_bin_inhibited)
error("part in an incorrect time-zone! p->ti_drift=%lld ti_drift=%lld",
c->hydro.parts[i].ti_drift, ti_drift);
#else
......@@ -1279,12 +1280,14 @@ void cell_check_gpart_drift_point(struct cell *c, void *data) {
c->grav.ti_old_part, ti_drift);
for (int i = 0; i < c->grav.count; ++i)
if (c->grav.parts[i].ti_drift != ti_drift)
if (c->grav.parts[i].ti_drift != ti_drift &&
c->grav.parts[i].time_bin != time_bin_inhibited)
error("g-part in an incorrect time-zone! gp->ti_drift=%lld ti_drift=%lld",
c->grav.parts[i].ti_drift, ti_drift);
for (int i = 0; i < c->stars.count; ++i)
if (c->stars.parts[i].ti_drift != ti_drift)
if (c->stars.parts[i].ti_drift != ti_drift &&
c->stars.parts[i].time_bin != time_bin_inhibited)
error("s-part in an incorrect time-zone! sp->ti_drift=%lld ti_drift=%lld",
c->stars.parts[i].ti_drift, ti_drift);
#else
......@@ -2745,6 +2748,8 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) {
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
if (c->end_force != NULL) scheduler_activate(s, c->end_force);
if (c->hydro.cooling != NULL) scheduler_activate(s, c->hydro.cooling);
if (c->hydro.star_formation != NULL)
scheduler_activate(s, c->hydro.star_formation);
if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
}
......@@ -2892,6 +2897,9 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
if (c->end_force != NULL) scheduler_activate(s, c->end_force);
if ((e->policy & engine_policy_cooling) && c->hydro.cooling != NULL)
scheduler_activate(s, c->hydro.cooling);
if ((e->policy & engine_policy_star_formation) &&
c->hydro.star_formation != NULL)
scheduler_activate(s, c->hydro.star_formation);
if (c->grav.down != NULL) scheduler_activate(s, c->grav.down);
if (c->grav.down_in != NULL) scheduler_activate(s, c->grav.down_in);
if (c->grav.mesh != NULL) scheduler_activate(s, c->grav.mesh);
......@@ -3263,6 +3271,9 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
struct part *const p = &parts[k];
struct xpart *const xp = &xparts[k];
/* Ignore inhibited particles */
if (part_is_inhibited(p, e)) continue;
/* Drift... */
drift_part(p, xp, dt_drift, dt_kick_hydro, dt_kick_grav, dt_therm,
ti_old_part, ti_current);
......@@ -3417,6 +3428,9 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
/* Get a handle on the gpart. */
struct gpart *const gp = &gparts[k];
/* Ignore inhibited particles */
if (gpart_is_inhibited(gp, e)) continue;
/* Drift... */
drift_gpart(gp, dt_drift, ti_old_gpart, ti_current);
......@@ -3462,6 +3476,9 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
/* Get a handle on the spart. */
struct spart *const sp = &sparts[k];
/* Ignore inhibited particles */
if (spart_is_inhibited(sp, e)) continue;
/* Drift... */
drift_spart(sp, dt_drift, ti_old_gpart, ti_current);
......@@ -3600,6 +3617,165 @@ void cell_check_timesteps(struct cell *c) {
#endif
}
/**
* @brief "Remove" a gas particle from the calculation.
*
* The particle is inhibited and will officially be removed at the next rebuild.
*
* @param e The #engine running on this node.
* @param c The #cell from which to remove the particle.
* @param p The #part to remove.
* @param xp The extended data of the particle to remove.
*/
void cell_remove_part(const struct engine *e, struct cell *c, struct part *p,
struct xpart *xp) {
/* Quick cross-check */
if (c->nodeID != e->nodeID)
error("Can't remove a particle in a foreign cell.");
/* Mark the particle as inhibited */
p->time_bin = time_bin_inhibited;
/* Mark the gpart as inhibited and stand-alone */
if (p->gpart) {
p->gpart->time_bin = time_bin_inhibited;
p->gpart->id_or_neg_offset = p->id;
p->gpart->type = swift_type_dark_matter;
}
/* Un-link the part */
p->gpart = NULL;
}
/**
* @brief "Remove" a gravity particle from the calculation.
*
* The particle is inhibited and will officially be removed at the next rebuild.
*
* @param e The #engine running on this node.
* @param c The #cell from which to remove the particle.
* @param gp The #gpart to remove.
*/
void cell_remove_gpart(const struct engine *e, struct cell *c,
struct gpart *gp) {
/* Quick cross-check */
if (c->nodeID != e->nodeID)
error("Can't remove a particle in a foreign cell.");
if (gp->type != swift_type_dark_matter)
error("Trying to remove a non-dark matter gpart.");
/* Mark the particle as inhibited */
gp->time_bin = time_bin_inhibited;
}
/**
* @brief "Remove" a star particle from the calculation.
*
* The particle is inhibited and will officially be removed at the next rebuild.
*
* @param e The #engine running on this node.
* @param c The #cell from which to remove the particle.
* @param sp The #spart to remove.
*/
void cell_remove_spart(const struct engine *e, struct cell *c,
struct spart *sp) {
/* Quick cross-check */
if (c->nodeID != e->nodeID)
error("Can't remove a particle in a foreign cell.");
/* Mark the particle as inhibited and stand-alone */
sp->time_bin = time_bin_inhibited;
if (sp->gpart) {
sp->gpart->time_bin = time_bin_inhibited;
sp->gpart->id_or_neg_offset = sp->id;
sp->gpart->type = swift_type_dark_matter;
}
/* Un-link the spart */
sp->gpart = NULL;
}
/**
* @brief "Remove" a gas particle from the calculation and convert its gpart
* friend to a dark matter particle.
*
* The particle is inhibited and will officially be removed at the next rebuild.
*
* @param e The #engine running on this node.
* @param c The #cell from which to remove the particle.
* @param p The #part to remove.
* @param xp The extended data of the particle to remove.
*/
void cell_convert_part_to_gpart(const struct engine *e, struct cell *c,
struct part *p, struct xpart *xp) {
/* Quick cross-checks */
if (c->nodeID != e->nodeID)
error("Can't remove a particle in a foreign cell.");
if (p->gpart == NULL)
error("Trying to convert part without gpart friend to dark matter!");
/* Get a handle */
struct gpart *gp = p->gpart;
/* Mark the particle as inhibited */
p->time_bin = time_bin_inhibited;
/* Un-link the part */
p->gpart = NULL;
/* Mark the gpart as dark matter */
gp->type = swift_type_dark_matter;
gp->id_or_neg_offset = p->id;
#ifdef SWIFT_DEBUG_CHECKS
gp->ti_kick = p->ti_kick;
#endif
}
/**
* @brief "Remove" a spart particle from the calculation and convert its gpart
* friend to a dark matter particle.
*
* The particle is inhibited and will officially be removed at the next rebuild.
*
* @param e The #engine running on this node.
* @param c The #cell from which to remove the particle.
* @param sp The #spart to remove.
*/
void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c,
struct spart *sp) {
/* Quick cross-check */
if (c->nodeID != e->nodeID)
error("Can't remove a particle in a foreign cell.");
if (sp->gpart == NULL)
error("Trying to convert spart without gpart friend to dark matter!");
/* Get a handle */
struct gpart *gp = sp->gpart;
/* Mark the particle as inhibited */
sp->time_bin = time_bin_inhibited;
/* Un-link the spart */
sp->gpart = NULL;
/* Mark the gpart as dark matter */
gp->type = swift_type_dark_matter;
gp->id_or_neg_offset = sp->id;
#ifdef SWIFT_DEBUG_CHECKS
gp->ti_kick = sp->ti_kick;
#endif
}
/**
* @brief Can we use the MM interactions fo a given pair of cells?
*
......
......@@ -270,6 +270,9 @@ struct cell {
/*! Number of #part updated in this cell. */
int updated;
/*! Number of #part inhibited in this cell. */
int inhibited;
/*! Is the #part data of this cell being used in a sub-cell? */
int hold;
......@@ -330,6 +333,9 @@ struct cell {
/*! Task for cooling */
struct task *cooling;
/*! Task for star formation */
struct task *star_formation;
#ifdef SWIFT_DEBUG_CHECKS
/*! Last (integer) time the cell's sort arrays were updated. */
......@@ -380,6 +386,9 @@ struct cell {
/*! Number of #gpart updated in this cell. */
int updated;
/*! Number of #gpart inhibited in this cell. */
int inhibited;
/*! Is the #gpart data of this cell being used in a sub-cell? */
int phold;
......@@ -460,6 +469,9 @@ struct cell {
/*! Number of #spart updated in this cell. */
int updated;
/*! Number of #spart inhibited in this cell. */
int inhibited;
/*! Is the #spart data of this cell being used in a sub-cell? */
int hold;
......@@ -635,6 +647,16 @@ void cell_activate_sorts(struct cell *c, int sid, struct scheduler *s);
void cell_clear_drift_flags(struct cell *c, void *data);
void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data);
int cell_has_tasks(struct cell *c);
void cell_remove_part(const struct engine *e, struct cell *c, struct part *p,
struct xpart *xp);
void cell_remove_gpart(const struct engine *e, struct cell *c,
struct gpart *gp);
void cell_remove_spart(const struct engine *e, struct cell *c,
struct spart *sp);
void cell_convert_part_to_gpart(const struct engine *e, struct cell *c,
struct part *p, struct xpart *xp);
void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c,
struct spart *sp);
int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj,
const struct engine *e, const struct space *s);
int cell_can_use_pair_mm_rebuild(const struct cell *ci, const struct cell *cj,
......
......@@ -36,7 +36,8 @@
/* Local collections for MPI reduces. */
struct mpicollectgroup1 {
long long updates, g_updates, s_updates;
long long updated, g_updated, s_updated;
long long inhibited, g_inhibited, s_inhibited;
integertime_t ti_hydro_end_min;
integertime_t ti_gravity_end_min;
int forcerebuild;
......@@ -85,9 +86,12 @@ void collectgroup1_apply(struct collectgroup1 *grp1, struct engine *e) {
e->ti_end_min = min(e->ti_hydro_end_min, e->ti_gravity_end_min);
e->ti_end_max = max(e->ti_hydro_end_max, e->ti_gravity_end_max);
e->ti_beg_max = max(e->ti_hydro_beg_max, e->ti_gravity_beg_max);
e->updates = grp1->updates;
e->g_updates = grp1->g_updates;
e->s_updates = grp1->s_updates;
e->updates = grp1->updated;
e->g_updates = grp1->g_updated;
e->s_updates = grp1->s_updated;
e->nr_inhibited_parts = grp1->inhibited;
e->nr_inhibited_gparts = grp1->g_inhibited;
e->nr_inhibited_sparts = grp1->s_inhibited;
e->forcerebuild = grp1->forcerebuild;
}
......@@ -95,10 +99,16 @@ void collectgroup1_apply(struct collectgroup1 *grp1, struct engine *e) {
* @brief Initialises a collectgroup1 struct ready for processing.
*
* @param grp1 The #collectgroup1 to initialise
* @param updates the number of updated hydro particles on this node this step.
* @param g_updates the number of updated gravity particles on this node this
* @param updated the number of updated hydro particles on this node this step.
* @param g_updated the number of updated gravity particles on this node this
* step.
* @param s_updated the number of updated star particles on this node this step.
* @param inhibited the number of inhibited hydro particles on this node this
* step.
* @param g_inhibited the number of inhibited gravity particles on this node
* this step.
* @param s_inhibited the number of inhibited star particles on this node this
* step.
* @param s_updates the number of updated star particles on this node this step.
* @param ti_hydro_end_min the minimum end time for next hydro time step after
* this step.
* @param ti_hydro_end_max the maximum end time for next hydro time step after
......@@ -113,17 +123,22 @@ void collectgroup1_apply(struct collectgroup1 *grp1, struct engine *e) {
* after this step.
* @param forcerebuild whether a rebuild is required after this step.
*/
void collectgroup1_init(struct collectgroup1 *grp1, size_t updates,
size_t g_updates, size_t s_updates,
void collectgroup1_init(struct collectgroup1 *grp1, size_t updated,
size_t g_updated, size_t s_updated, size_t inhibited,
size_t g_inhibited, size_t s_inhibited,
integertime_t ti_hydro_end_min,
integertime_t ti_hydro_end_max,
integertime_t ti_hydro_beg_max,
integertime_t ti_gravity_end_min,
integertime_t ti_gravity_end_max,
integertime_t ti_gravity_beg_max, int forcerebuild) {
grp1->updates = updates;
grp1->g_updates = g_updates;
grp1->s_updates = s_updates;
grp1->updated = updated;
grp1->g_updated = g_updated;
grp1->s_updated = s_updated;
grp1->inhibited = inhibited;
grp1->g_inhibited = g_inhibited;
grp1->s_inhibited = s_inhibited;
grp1->ti_hydro_end_min = ti_hydro_end_min;
grp1->ti_hydro_end_max = ti_hydro_end_max;
grp1->ti_hydro_beg_max = ti_hydro_beg_max;
......@@ -147,9 +162,12 @@ void collectgroup1_reduce(struct collectgroup1 *grp1) {
/* Populate an MPI group struct and reduce this across all nodes. */
struct mpicollectgroup1 mpigrp11;
mpigrp11.updates = grp1->updates;
mpigrp11.g_updates = grp1->g_updates;
mpigrp11.s_updates = grp1->s_updates;
mpigrp11.updated = grp1->updated;
mpigrp11.g_updated = grp1->g_updated;
mpigrp11.s_updated = grp1->s_updated;
mpigrp11.inhibited = grp1->inhibited;
mpigrp11.g_inhibited = grp1->g_inhibited;
mpigrp11.s_inhibited = grp1->s_inhibited;
mpigrp11.ti_hydro_end_min = grp1->ti_hydro_end_min;
mpigrp11.ti_gravity_end_min = grp1->ti_gravity_end_min;
mpigrp11.forcerebuild = grp1->forcerebuild;
......@@ -160,9 +178,12 @@ void collectgroup1_reduce(struct collectgroup1 *grp1) {
error("Failed to reduce mpicollection1.");
/* And update. */
grp1->updates = mpigrp12.updates;
grp1->g_updates = mpigrp12.g_updates;
grp1->s_updates = mpigrp12.s_updates;
grp1->updated = mpigrp12.updated;
grp1->g_updated = mpigrp12.g_updated;
grp1->s_updated = mpigrp12.s_updated;
grp1->inhibited = mpigrp12.inhibited;
grp1->g_inhibited = mpigrp12.g_inhibited;
grp1->s_inhibited = mpigrp12.s_inhibited;
grp1->ti_hydro_end_min = mpigrp12.ti_hydro_end_min;
grp1->ti_gravity_end_min = mpigrp12.ti_gravity_end_min;
grp1->forcerebuild = mpigrp12.forcerebuild;
......@@ -182,9 +203,14 @@ static void doreduce1(struct mpicollectgroup1 *mpigrp11,
/* Do what is needed for each part of the collection. */
/* Sum of updates. */
mpigrp11->updates += mpigrp12->updates;
mpigrp11->g_updates += mpigrp12->g_updates;
mpigrp11->s_updates += mpigrp12->s_updates;
mpigrp11->updated += mpigrp12->updated;
mpigrp11->g_updated += mpigrp12->g_updated;
mpigrp11->s_updated += mpigrp12->s_updated;
/* Sum of inhibited */
mpigrp11->inhibited += mpigrp12->inhibited;
mpigrp11->g_inhibited += mpigrp12->g_inhibited;
mpigrp11->s_inhibited += mpigrp12->s_inhibited;
/* Minimum end time. */
mpigrp11->ti_hydro_end_min =
......
......@@ -35,7 +35,10 @@ struct engine;
struct collectgroup1 {
/* Number of particles updated */
long long updates, g_updates, s_updates;
long long updated, g_updated, s_updated;
/* Number of particles inhibited */
long long inhibited, g_inhibited, s_inhibited;
/* Times for the time-step */
integertime_t ti_hydro_end_min, ti_hydro_end_max, ti_hydro_beg_max;
......@@ -47,8 +50,9 @@ struct collectgroup1 {
void collectgroup_init(void);
void collectgroup1_apply(struct collectgroup1 *grp1, struct engine *e);
void collectgroup1_init(struct collectgroup1 *grp1, size_t updates,
size_t g_updates, size_t s_updates,
void collectgroup1_init(struct collectgroup1 *grp1, size_t updated,
size_t g_updated, size_t s_updated, size_t inhibited,
size_t g_inhibited, size_t s_inhibited,
integertime_t ti_hydro_end_min,
integertime_t ti_hydro_end_max,
integertime_t ti_hydro_beg_max,
......
......@@ -787,35 +787,109 @@ void io_duplicate_stars_gparts(struct threadpool* tp,
}
/**
* @brief Copy every DM #gpart into the dmparts array.
* @brief Copy every non-inhibited #part into the parts_written array.
*
* @param parts The array of #part containing all particles.
* @param xparts The array of #xpart containing all particles.
* @param parts_written The array of #part to fill with particles we want to
* write.
* @param xparts_written The array of #xpart to fill with particles we want to
* write.
* @param Nparts The total number of #part.
* @param Nparts_written The total number of #part to write.
*/
void io_collect_parts_to_write(const struct part* restrict parts,
const struct xpart* restrict xparts,
struct part* restrict parts_written,
struct xpart* restrict xparts_written,
const size_t Nparts,
const size_t Nparts_written) {
size_t count = 0;
/* Loop over all parts */
for (size_t i = 0; i < Nparts; ++i) {
/* And collect the ones that have not been removed */
if (parts[i].time_bin != time_bin_inhibited) {
parts_written[count] = parts[i];
xparts_written[count] = xparts[i];
count++;
}
}
/* Check that everything is fine */
if (count != Nparts_written)
error("Collected the wrong number of particles (%zu vs. %zu expected)",
count, Nparts_written);
}
/**
* @brief Copy every non-inhibited #spart into the sparts_written array.
*
* @param sparts The array of #spart containing all particles.
* @param sparts_written The array of #spart to fill with particles we want to
* write.
* @param Nsparts The total number of #part.
* @param Nsparts_written The total number of #part to write.
*/
void io_collect_sparts_to_write(const struct spart* restrict sparts,
struct spart* restrict sparts_written,
const size_t Nsparts,
const size_t Nsparts_written) {
size_t count = 0;
/* Loop over all parts */
for (size_t i = 0; i < Nsparts; ++i) {
/* And collect the ones that have not been removed */
if (sparts[i].time_bin != time_bin_inhibited) {
sparts_written[count] = sparts[i];
count++;
}
}
/* Check that everything is fine */
if (count != Nsparts_written)
error("Collected the wrong number of s-particles (%zu vs. %zu expected)",
count, Nsparts_written);
}
/**
* @brief Copy every non-inhibited DM #gpart into the gparts_written array.
*
* @param gparts The array of #gpart containing all particles.
* @param Ntot The number of #gpart.
* @param dmparts The array of #gpart containg DM particles to be filled.
* @param Ndm The number of DM particles.
* @param gparts_written The array of #gpart to fill with particles we want to
* write.
* @param Ngparts The total number of #part.
* @param Ngparts_written The total number of #part to write.
*/
void io_collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
struct gpart* const dmparts, size_t Ndm) {
void io_collect_gparts_to_write(const struct gpart* restrict gparts,
struct gpart* restrict gparts_written,
const size_t Ngparts,
const size_t Ngparts_written) {
size_t count = 0;
/* Loop over all gparts */
for (size_t i = 0; i < Ntot; ++i) {
/* Loop over all parts */
for (size_t i = 0; i < Ngparts; ++i) {
/* message("i=%zd count=%zd id=%lld part=%p", i, count, gparts[i].id,
* gparts[i].part); */
/* And collect the ones that have not been removed */
if ((gparts[i].time_bin != time_bin_inhibited) &&