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

Merge branch 'add_particles' into 'master'

Star formation (non-MPI)

See merge request !688
parents 5bcc4d36 7ba384f4
......@@ -6,7 +6,7 @@
/____/ |__/|__/___/_/ /_/
SPH With Inter-dependent Fine-grained Tasking
Website: www.swiftsim.com
Website: www.swiftsim.com
Twitter: @SwiftSimulation
See INSTALL.swift for install instructions.
......@@ -27,7 +27,7 @@ Parameters:
-D, --drift-all Always drift all particles even the ones
far from active particles. This emulates
Gadget-[23] and GIZMO's default behaviours.
-F, --sourceterms
-F, --star-formation Run with star formation
-g, --external-gravity Run with an external gravitational potential.
-G, --self-gravity Run with self-gravity.
-M, --multipole-reconstruction Reconstruct the multipoles every time-step.
......
......@@ -75,7 +75,7 @@ Parameters:
-D, --drift-all Always drift all particles even the ones
far from active particles. This emulates
Gadget-[23] and GIZMO's default behaviours.
-F, --sourceterms
-F, --star-formation Run with star formation
-g, --external-gravity Run with an external gravitational potential.
-G, --self-gravity Run with self-gravity.
-M, --multipole-reconstruction Reconstruct the multipoles every time-step.
......
......@@ -22,7 +22,7 @@ can be found by typing ``./swift -h``::
-D, --drift-all Always drift all particles even the ones
far from active particles. This emulates
Gadget-[23] and GIZMO's default behaviours.
-F, --sourceterms
-F, --star-formation Run with star formation
-g, --external-gravity Run with an external gravitational potential.
-G, --self-gravity Run with self-gravity.
-M, --multipole-reconstruction Reconstruct the multipoles every time-step.
......
......@@ -52,7 +52,7 @@
/* Global profiler. */
struct profiler prof;
// Usage string.
/* Usage string. */
static const char *const swift_usage[] = {
"swift [options] [[--] param-file]",
"swift [options] param-file",
......@@ -61,7 +61,7 @@ static const char *const swift_usage[] = {
NULL,
};
// Function to handle multiple -P arguments.
/* Function to handle multiple -P arguments. */
struct cmdparams {
const char *param[PARSER_MAX_NO_OF_PARAMS];
int nparam;
......@@ -97,7 +97,6 @@ int main(int argc, char *argv[]) {
struct stars_props stars_properties;
struct part *parts = NULL;
struct phys_const prog_const;
struct sourceterms sourceterms;
struct space s;
struct spart *sparts = NULL;
struct unit_system us;
......@@ -147,11 +146,11 @@ int main(int argc, char *argv[]) {
int restart = 0;
int with_cosmology = 0;
int with_external_gravity = 0;
int with_sourceterms = 0;
int with_cooling = 0;
int with_self_gravity = 0;
int with_hydro = 0;
int with_stars = 0;
int with_star_formation = 0;
int with_feedback = 0;
int with_fp_exceptions = 0;
int with_drift_all = 0;
......@@ -186,7 +185,8 @@ int main(int argc, char *argv[]) {
"particles. This emulates Gadget-[23] and GIZMO's default "
"behaviours.",
NULL, 0, 0),
OPT_BOOLEAN('F', "sourceterms", &with_sourceterms, "", NULL, 0, 0),
OPT_BOOLEAN('F', "star-formation", &with_star_formation,
"Run with star formation", NULL, 0, 0),
OPT_BOOLEAN('g', "external-gravity", &with_external_gravity,
"Run with an external gravitational potential.", NULL, 0, 0),
OPT_BOOLEAN('G', "self-gravity", &with_self_gravity,
......@@ -449,10 +449,9 @@ int main(int argc, char *argv[]) {
#ifdef WITH_MPI
if (with_mpole_reconstruction && nr_nodes > 1)
error("Cannot reconstruct m-poles every step over MPI (yet).");
#endif
#ifdef WITH_MPI
if (with_feedback) error("Can't run with feedback over MPI (yet).");
if (with_star_formation)
error("Can't run with star formation over MPI (yet)");
#endif
#if defined(WITH_MPI) && defined(HAVE_VELOCIRAPTOR)
......@@ -779,7 +778,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
space_init(&s, params, &cosmo, dim, parts, gparts, sparts, Ngas, Ngpart,
Nspart, periodic, replicate, generate_gas_in_ics, with_hydro,
with_self_gravity, talking, dry_run);
with_self_gravity, with_star_formation, talking, dry_run);
if (myrank == 0) {
clocks_gettime(&toc);
......@@ -860,10 +859,6 @@ int main(int argc, char *argv[]) {
chemistry_init(params, &us, &prog_const, &chemistry);
if (myrank == 0) chemistry_print(&chemistry);
/* Initialise the feedback properties */
if (with_sourceterms) sourceterms_init(params, &us, &sourceterms);
if (with_sourceterms && myrank == 0) sourceterms_print(&sourceterms);
/* Construct the engine policy */
int engine_policies = ENGINE_POLICY | engine_policy_steal;
if (with_drift_all) engine_policies |= engine_policy_drift_all;
......@@ -875,21 +870,18 @@ int main(int argc, char *argv[]) {
engine_policies |= engine_policy_external_gravity;
if (with_cosmology) engine_policies |= engine_policy_cosmology;
if (with_cooling) engine_policies |= engine_policy_cooling;
if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
if (with_stars) engine_policies |= engine_policy_stars;
if (with_star_formation) engine_policies |= engine_policy_star_formation;
if (with_feedback) engine_policies |= engine_policy_feedback;
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],
engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
&hydro_properties, &gravity_properties, &stars_properties,
&mesh, &potential, &cooling_func, &chemistry, &sourceterms);
&mesh, &potential, &cooling_func, &chemistry);
engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
talking, restart_file);
......
......@@ -63,6 +63,9 @@ Scheduler:
cell_sub_size_self_stars: 32000 # (Optional) Maximal number of interactions per sub-self stars task (this is the default value).
cell_split_size: 400 # (Optional) Maximal number of particles per cell (this is the default value).
cell_subdepth_diff_grav: 4 # (Optional) Maximal depth difference between leaves and a cell that gravity tasks can be pushed down to (this is the default value).
cell_extra_parts: 0 # (Optional) Number of spare parts per top-level allocated at rebuild time for on-the-fly creation.
cell_extra_gparts: 0 # (Optional) Number of spare gparts per top-level allocated at rebuild time for on-the-fly creation.
cell_extra_sparts: 400 # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
max_top_level_cells: 12 # (Optional) Maximal number of top-level cells in any dimension. The number of top-level cells will be the cube of this (this is the default value).
tasks_per_cell: 0 # (Optional) The average number of tasks per cell. If not large enough the simulation will fail (means guess...).
mpi_message_limit: 4096 # (Optional) Maximum MPI task message size to send non-buffered, KB.
......
......@@ -44,7 +44,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \
physical_constants.h physical_constants_cgs.h potential.h version.h \
hydro_properties.h riemann.h threadpool.h cooling_io.h cooling.h cooling_struct.h \
sourceterms.h sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \
chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \
......@@ -63,7 +63,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c engine_maketasks.c
proxy.c parallel_io.c units.c common_io.c single_io.c multipole.c version.c map.c \
kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \
physical_constants.c potential.c hydro_properties.c \
threadpool.c cooling.c sourceterms.c \
threadpool.c cooling.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
part_type.c xmf.c gravity_properties.c gravity.c \
collectgroup.c hydro_space.c equation_of_state.c \
......@@ -81,7 +81,6 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
gravity/Potential/gravity.h gravity/Potential/gravity_iact.h gravity/Potential/gravity_io.h \
gravity/Potential/gravity_debug.h gravity/Potential/gravity_part.h \
sourceterms.h \
equation_of_state.h \
equation_of_state/ideal_gas/equation_of_state.h equation_of_state/isothermal/equation_of_state.h \
hydro.h hydro_io.h \
......
......@@ -972,6 +972,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
/* Store the counts and offsets. */
for (int k = 0; k < 8; k++) {
c->progeny[k]->hydro.count = bucket_count[k];
c->progeny[k]->hydro.count_total = c->progeny[k]->hydro.count;
c->progeny[k]->hydro.parts = &c->hydro.parts[bucket_offset[k]];
c->progeny[k]->hydro.xparts = &c->hydro.xparts[bucket_offset[k]];
}
......@@ -1089,6 +1090,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
/* Store the counts and offsets. */
for (int k = 0; k < 8; k++) {
c->progeny[k]->stars.count = bucket_count[k];
c->progeny[k]->stars.count_total = c->progeny[k]->stars.count;
c->progeny[k]->stars.parts = &c->stars.parts[bucket_offset[k]];
}
......@@ -1151,6 +1153,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
/* Store the counts and offsets. */
for (int k = 0; k < 8; k++) {
c->progeny[k]->grav.count = bucket_count[k];
c->progeny[k]->grav.count_total = c->progeny[k]->grav.count;
c->progeny[k]->grav.parts = &c->grav.parts[bucket_offset[k]];
}
}
......@@ -3813,6 +3816,190 @@ void cell_check_timesteps(struct cell *c) {
#endif
}
void cell_check_spart_pos(const struct cell *c,
const struct spart *global_sparts) {
#ifdef SWIFT_DEBUG_CHECKS
/* Recurse */
if (c->split) {
for (int k = 0; k < 8; ++k)
if (c->progeny[k] != NULL)
cell_check_spart_pos(c->progeny[k], global_sparts);
}
struct spart *sparts = c->stars.parts;
const int count = c->stars.count;
for (int i = 0; i < count; ++i) {
const struct spart *sp = &sparts[i];
if ((sp->x[0] < c->loc[0] / space_stretch) ||
(sp->x[1] < c->loc[1] / space_stretch) ||
(sp->x[2] < c->loc[2] / space_stretch) ||
(sp->x[0] >= (c->loc[0] + c->width[0]) * space_stretch) ||
(sp->x[1] >= (c->loc[1] + c->width[1]) * space_stretch) ||
(sp->x[2] >= (c->loc[2] + c->width[2]) * space_stretch))
error("spart not in its cell!");
if (sp->time_bin != time_bin_not_created &&
sp->time_bin != time_bin_inhibited) {
const struct gpart *gp = sp->gpart;
if (gp == NULL && sp->time_bin != time_bin_not_created)
error("Unlinked spart!");
if (&global_sparts[-gp->id_or_neg_offset] != sp)
error("Incorrectly linked spart!");
}
}
#else
error("Calling a degugging function outside debugging mode.");
#endif
}
/**
* @brief Recursively update the pointer and counter for #spart after the
* addition of a new particle.
*
* @param c The cell we are working on.
* @param progeny_list The list of the progeny index at each level for the
* leaf-cell where the particle was added.
* @param main_branch Are we in a cell directly above the leaf where the new
* particle was added?
*/
void cell_recursively_shift_sparts(struct cell *c,
const int progeny_list[space_cell_maxdepth],
const int main_branch) {
if (c->split) {
/* No need to recurse in progenies located before the insestion point */
const int first_progeny = main_branch ? progeny_list[(int)c->depth] : 0;
for (int k = first_progeny; k < 8; ++k) {
if (c->progeny[k] != NULL)
cell_recursively_shift_sparts(c->progeny[k], progeny_list,
main_branch && (k == first_progeny));
}
}
/* When directly above the leaf with the new particle: increase the particle
* count */
/* When after the leaf with the new particle: shift by one position */
if (main_branch)
c->stars.count++;
else
c->stars.parts++;
}
/**
* @brief "Add" a #spart in a given #cell.
*
* This function will a a #spart at the start of the current cell's array by
* shifting all the #spart in the top-level cell by one position. All the
* pointers and cell counts are updated accordingly.
*
* @param e The #engine.
* @param c The leaf-cell in which to add the #spart.
*
* @return A pointer to the newly added #spart. The spart has a been zeroed and
* given a position within the cell as well as set to the minimal active time
* bin.
*/
struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
/* Perform some basic consitency checks */
if (c->nodeID != engine_rank) error("Adding spart on a foreign node");
if (c->grav.ti_old_part != e->ti_current) error("Undrifted cell!");
if (c->split) error("Addition of spart performed above the leaf level");
/* Progeny number at each level */
int progeny[space_cell_maxdepth];
#ifdef SWIFT_DEBUG_CHECKS
for (int i = 0; i < space_cell_maxdepth; ++i) progeny[i] = -1;
#endif
/* Get the top-level this leaf cell is in and compute the progeny indices at
each level */
struct cell *top = c;
while (top->parent != NULL) {
for (int k = 0; k < 8; ++k) {
if (top->parent->progeny[k] == top) {
progeny[(int)top->parent->depth] = k;
}
}
top = top->parent;
}
/* Are there any extra particles left? */
if (top->stars.count == top->stars.count_total - 1) {
message("We ran out of star particles!");
atomic_inc(&e->forcerebuild);
return NULL;
}
/* Number of particles to shift in order to get a free space. */
const size_t n_copy = &top->stars.parts[top->stars.count] - c->stars.parts;
#ifdef SWIFT_DEBUG_CHECKS
if (c->stars.parts + n_copy > top->stars.parts + top->stars.count)
error("Copying beyond the allowed range");
#endif
if (n_copy > 0) {
// MATTHIEU: This can be improved. We don't need to copy everything, just
// need to swap a few particles.
memmove(&c->stars.parts[1], &c->stars.parts[0],
n_copy * sizeof(struct spart));
/* Update the gpart->spart links (shift by 1) */
for (size_t i = 0; i < n_copy; ++i) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->stars.parts[i + 1].gpart == NULL) {
error("Incorrectly linked spart!");
}
#endif
c->stars.parts[i + 1].gpart->id_or_neg_offset--;
}
}
/* Recursively shift all the stars to get a free spot at the start of the
* current cell*/
cell_recursively_shift_sparts(top, progeny, /* main_branch=*/1);
/* We now have an empty spart as the first particle in that cell */
struct spart *sp = &c->stars.parts[0];
bzero(sp, sizeof(struct spart));
/* Give it a decent position */
sp->x[0] = c->loc[0] + 0.5 * c->width[0];
sp->x[1] = c->loc[1] + 0.5 * c->width[1];
sp->x[2] = c->loc[2] + 0.5 * c->width[2];
/* Set it to the current time-bin */
sp->time_bin = e->min_active_bin;
top = c;
while (top->parent != NULL) {
top->grav.ti_end_min = e->ti_current;
top = top->parent;
}
top->grav.ti_end_min = e->ti_current;
#ifdef SWIFT_DEBUG_CHECKS
/* Specify it was drifted to this point */
sp->ti_drift = e->ti_current;
#endif
/* Register that we used one of the free slots. */
const size_t one = 1;
atomic_sub(&e->s->nr_extra_sparts, one);
return sp;
}
/**
* @brief "Remove" a gas particle from the calculation.
*
......@@ -3899,15 +4086,21 @@ void cell_remove_spart(const struct engine *e, struct cell *c,
* @brief "Remove" a gas particle from the calculation and convert its gpart
* friend to a dark matter particle.
*
* Note that the #part is not destroyed. The pointer is still valid
* after this call and the properties of the #part are not altered
* apart from the time-bin and #gpart pointer.
* 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.
*
* @return Pointer to the #gpart the #part has become. It carries the
* ID of the #part and has a dark matter type.
*/
void cell_convert_part_to_gpart(const struct engine *e, struct cell *c,
struct part *p, struct xpart *xp) {
struct gpart *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)
......@@ -3932,20 +4125,28 @@ void cell_convert_part_to_gpart(const struct engine *e, struct cell *c,
#ifdef SWIFT_DEBUG_CHECKS
gp->ti_kick = p->ti_kick;
#endif
return gp;
}
/**
* @brief "Remove" a spart particle from the calculation and convert its gpart
* friend to a dark matter particle.
*
* Note that the #spart is not destroyed. The pointer is still valid
* after this call and the properties of the #spart are not altered
* apart from the time-bin and #gpart pointer.
* 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.
*
* @return Pointer to the #gpart the #spart has become. It carries the
* ID of the #spart and has a dark matter type.
*/
void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c,
struct spart *sp) {
struct gpart *cell_convert_spart_to_gpart(const struct engine *e,
struct cell *c, struct spart *sp) {
/* Quick cross-check */
if (c->nodeID != e->nodeID)
......@@ -3970,6 +4171,210 @@ void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c,
#ifdef SWIFT_DEBUG_CHECKS
gp->ti_kick = sp->ti_kick;
#endif
return gp;
}
/**
* @brief "Remove" a #part from a #cell and replace it with a #spart
* connected to the same #gpart.
*
* Note that the #part is not destroyed. The pointer is still valid
* after this call and the properties of the #part are not altered
* apart from the time-bin and #gpart pointer.
* The particle is inhibited and will officially be removed at the next rebuild.
*
* @param e The #engine.
* @param c The #cell from which to remove the #part.
* @param p The #part to remove (must be inside c).
* @param xp The extended data of the #part.
*
* @return A fresh #spart with the same ID, position, velocity and
* time-bin as the original #part.
*/
struct spart *cell_convert_part_to_spart(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.");
if (p->gpart == NULL)
error("Trying to convert part without gpart friend to star!");
/* Create a fresh (empty) spart */
struct spart *sp = cell_add_spart(e, c);
/* Did we run out of free spart slots? */
if (sp == NULL) return NULL;
/* Destroy the gas particle and get it's gpart friend */
struct gpart *gp = cell_convert_part_to_gpart(e, c, p, xp);
/* Assign the ID back */
sp->id = gp->id_or_neg_offset;
gp->type = swift_type_stars;
/* Re-link things */
sp->gpart = gp;
gp->id_or_neg_offset = -(sp - e->s->sparts);
/* Synchronize clocks */
gp->time_bin = sp->time_bin;
/* Synchronize masses, positions and velocities */
sp->mass = gp->mass;
sp->x[0] = gp->x[0];
sp->x[1] = gp->x[1];
sp->x[2] = gp->x[2];
sp->v[0] = gp->v_full[0];
sp->v[1] = gp->v_full[1];
sp->v[2] = gp->v_full[2];
#ifdef SWIFT_DEBUG_CHECKS
sp->ti_kick = gp->ti_kick;
gp->ti_drift = sp->ti_drift;
#endif
/* Set a smoothing length */
sp->h = max(c->stars.h_max, c->hydro.h_max);
/* Here comes the Sun! */
return sp;
}
/**
* @brief Re-arrange the #part in a top-level cell such that all the extra ones
* for on-the-fly creation are located at the end of the array.
*
* @param c The #cell to sort.
* @param parts_offset The offset between the first #part in the array and the
* first #part in the global array in the space structure (for re-linking).
*/
void cell_reorder_extra_parts(struct cell *c, const ptrdiff_t parts_offset) {
struct part *parts = c->hydro.parts;
struct xpart *xparts = c->hydro.xparts;
const int count_real = c->hydro.count;
if (c->depth != 0 || c->nodeID != engine_rank)
error("This function should only be called on local top-level cells!");
int first_not_extra = count_real;
/* Find extra particles */
for (int i = 0; i < count_real; ++i) {
if (parts[i].time_bin == time_bin_not_created) {
/* Find the first non-extra particle after the end of the
real particles */
while (parts[first_not_extra].time_bin == time_bin_not_created) {
++first_not_extra;
}
#ifdef SWIFT_DEBUG_CHECKS
if (first_not_extra >= count_real + space_extra_parts)
error("Looking for extra particles beyond this cell's range!");
#endif
/* Swap everything, including g-part pointer */
memswap(&parts[i], &parts[first_not_extra], sizeof(struct part));
memswap(&xparts[i], &xparts[first_not_extra], sizeof(struct xpart));
if (parts[i].gpart)
parts[i].gpart->id_or_neg_offset = -(i + parts_offset);
}
}
}
/**
* @brief Re-arrange the #spart in a top-level cell such that all the extra ones
* for on-the-fly creation are located at the end of the array.
*
* @param c The #cell to sort.
* @param sparts_offset The offset between the first #spart in the array and the
* first #spart in the global array in the space structure (for re-linking).
*/
void cell_reorder_extra_sparts(struct cell *c, const ptrdiff_t sparts_offset) {
struct spart *sparts = c->stars.parts;
const int count_real = c->stars.count;
if (c->depth != 0 || c->nodeID != engine_rank)
error("This function should only be called on local top-level cells!");
int first_not_extra = count_real;
/* Find extra particles */
for (int i = 0; i < count_real; ++i) {
if (sparts[i].time_bin == time_bin_not_created) {
/* Find the first non-extra particle after the end of the
real particles */
while (sparts[first_not_extra].time_bin == time_bin_not_created) {
++first_not_extra;
}
#ifdef SWIFT_DEBUG_CHECKS
if (first_not_extra >= count_real + space_extra_sparts)
error("Looking for extra particles beyond this cell's range!");
#endif
/* Swap everything, including g-part pointer */
memswap(&sparts[i], &sparts[first_not_extra], sizeof(struct spart));
if (sparts[i].gpart)
sparts[i].gpart->id_or_neg_offset = -(i + sparts_offset);
sparts[first_not_extra].gpart = NULL;
#ifdef SWIFT_DEBUG_CHECKS
if (sparts[first_not_extra].time_bin != time_bin_not_created)
error("Incorrect swap occured!");
#endif
}
}
}
/**
* @brief Re-arrange the #gpart in a top-level cell such that all the extra ones
* for on-the-fly creation are located at the end of the array.
*
* @param c The #cell to sort.
* @param parts The global array of #part (for re-linking).
* @param sparts The global array of #spart (for re-linking).
*/
void cell_reorder_extra_gparts(struct cell *c, struct part *parts,
struct spart *sparts) {