diff --git a/README b/README
index 97d096ef1804a4348e88dfaf67c22b2c427a3ad8..5b0b5f3ebf8b322e121e0398ca01062e226be649 100644
--- a/README
+++ b/README
@@ -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.
diff --git a/README.md b/README.md
index 94e95776cd80c1bb822f0f68290c2add9d2bb58b..92925017d9a8f44521a542e1df3760cc58ba73cb 100644
--- a/README.md
+++ b/README.md
@@ -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.
diff --git a/doc/RTD/source/CommandLineOptions/index.rst b/doc/RTD/source/CommandLineOptions/index.rst
index 9fb9d784d6057e4d9aa4a923143d622e577f142c..de969730bebc6c20950086cb1ba9a20b94e01af6 100644
--- a/doc/RTD/source/CommandLineOptions/index.rst
+++ b/doc/RTD/source/CommandLineOptions/index.rst
@@ -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.
diff --git a/examples/main.c b/examples/main.c
index 93074a637e7427ac7bec53c51d7e6e608b6cdff2..db7df3b4a129b5469758485c7a042af173e63123 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -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);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 6adccf2963dbeff67755bdac946e7bfb10d4a897..c8a3d8b0e8271d622536e4a9550ea7646406a259 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -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.
diff --git a/src/Makefile.am b/src/Makefile.am
index e0c5dbe48d22d0340e2ab13a2dfdc833dd115908..114a3dc4d794c302b06c4339d35ed88cbfd13035 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -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 \
@@ -64,7 +64,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 \
@@ -82,7 +82,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 \
diff --git a/src/cell.c b/src/cell.c
index 3fe5e21e7c888f2358395a27e13710db460fd74c..ed54163a008f7fad8294f50458ef8d2004be832d 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -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) {
+
+ struct gpart *gparts = c->grav.parts;
+ const int count_real = c->grav.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 (gparts[i].time_bin == time_bin_not_created) {
+
+ /* Find the first non-extra particle after the end of the
+ real particles */
+ while (gparts[first_not_extra].time_bin == time_bin_not_created) {
+ ++first_not_extra;
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (first_not_extra >= count_real + space_extra_gparts)
+ error("Looking for extra particles beyond this cell's range!");
+#endif
+
+ /* Swap everything (including pointers) */
+ memswap(&gparts[i], &gparts[first_not_extra], sizeof(struct gpart));
+ if (gparts[i].type == swift_type_gas) {
+ parts[-gparts[i].id_or_neg_offset].gpart = &gparts[i];
+ } else if (gparts[i].type == swift_type_stars) {
+ sparts[-gparts[i].id_or_neg_offset].gpart = &gparts[i];
+ }
+ }
+ }
}
/**
diff --git a/src/cell.h b/src/cell.h
index 97ca22e584c67de20ca0826425f60523b8158ffa..7178698c425c3791223789930c82761d13fc0116 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -220,12 +220,12 @@ struct cell {
/*! The cell dimensions. */
double width[3];
- /*! Linking pointer for "memory management". */
- struct cell *next;
-
/*! Pointers to the next level of cells. */
struct cell *progeny[8];
+ /*! Linking pointer for "memory management". */
+ struct cell *next;
+
/*! Parent cell. */
struct cell *parent;
@@ -248,18 +248,45 @@ struct cell {
* pair/self tasks */
struct cell *super;
- /*! Last (integer) time the cell's part were drifted forward in time. */
- integertime_t ti_old_part;
+ /*! The task computing this cell's sorts. */
+ struct task *sorts;
- /*! Maximum part movement in this cell since last construction. */
- float dx_max_part;
+ /*! The drift task for parts */
+ struct task *drift;
- /*! Maximum particle movement in this cell since the last sort. */
- float dx_max_sort;
+ /*! Linked list of the tasks computing this cell's hydro density. */
+ struct link *density;
+
+ /* Linked list of the tasks computing this cell's hydro gradients. */
+ struct link *gradient;
+
+ /*! Linked list of the tasks computing this cell's hydro forces. */
+ struct link *force;
+
+ /*! Dependency implicit task for the ghost (in->ghost->out)*/
+ struct task *ghost_in;
+
+ /*! Dependency implicit task for the ghost (in->ghost->out)*/
+ struct task *ghost_out;
+
+ /*! The ghost task itself */
+ struct task *ghost;
+
+ /*! The extra ghost task for complex hydro schemes */
+ struct task *extra_ghost;
+
+ /*! Task for cooling */
+ struct task *cooling;
+
+ /*! Task for star formation */
+ struct task *star_formation;
/*! Max smoothing length in this cell. */
double h_max;
+ /*! Last (integer) time the cell's part were drifted forward in time. */
+ integertime_t ti_old_part;
+
/*! Minimum end of (integer) time step in this cell for hydro tasks. */
integertime_t ti_end_min;
@@ -270,20 +297,14 @@ struct cell {
*/
integertime_t ti_beg_max;
- /*! Nr of #part in this cell. */
- int count;
-
/*! Spin lock for various uses (#part case). */
swift_lock_type lock;
- /*! Number of #part updated in this cell. */
- int updated;
-
- /*! Number of #part inhibited in this cell. */
- int inhibited;
+ /*! Maximum part movement in this cell since last construction. */
+ float dx_max_part;
- /*! Is the #part data of this cell being used in a sub-cell? */
- int hold;
+ /*! Maximum particle movement in this cell since the last sort. */
+ float dx_max_sort;
/*! Values of h_max before the drifts, used for sub-cell tasks. */
float h_max_old;
@@ -294,12 +315,30 @@ struct cell {
/*! Values of dx_max_sort before the drifts, used for sub-cell tasks. */
float dx_max_sort_old;
+ /*! Nr of #part in this cell. */
+ int count;
+
+ /*! Nr of #part this cell can hold after addition of new #part. */
+ int count_total;
+
+ /*! 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;
+
/*! Bit mask of sort directions that will be needed in the next timestep. */
unsigned int requires_sorts;
/*! Bit mask of sorts that need to be computed for this cell. */
unsigned int do_sort;
+ /*! Bit-mask indicating the sorted directions */
+ unsigned int sorted;
+
/*! Does this cell need to be drifted (hydro)? */
char do_drift;
@@ -309,42 +348,6 @@ struct cell {
/*! Do any of this cell's sub-cells need to be sorted? */
char do_sub_sort;
- /*! Bit-mask indicating the sorted directions */
- unsigned int sorted;
-
- /*! The task computing this cell's sorts. */
- struct task *sorts;
-
- /*! The drift task for parts */
- struct task *drift;
-
- /*! Linked list of the tasks computing this cell's hydro density. */
- struct link *density;
-
- /* Linked list of the tasks computing this cell's hydro gradients. */
- struct link *gradient;
-
- /*! Linked list of the tasks computing this cell's hydro forces. */
- struct link *force;
-
- /*! Dependency implicit task for the ghost (in->ghost->out)*/
- struct task *ghost_in;
-
- /*! Dependency implicit task for the ghost (in->ghost->out)*/
- struct task *ghost_out;
-
- /*! The ghost task itself */
- struct task *ghost;
-
- /*! The extra ghost task for complex hydro schemes */
- struct task *extra_ghost;
-
- /*! 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. */
@@ -367,6 +370,36 @@ struct cell {
* tasks */
struct cell *super;
+ /*! The drift task for gparts */
+ struct task *drift;
+
+ /*! Implicit task (going up- and down the tree) for the #gpart drifts */
+ struct task *drift_out;
+
+ /*! Linked list of the tasks computing this cell's gravity forces. */
+ struct link *grav;
+
+ /*! Linked list of the tasks computing this cell's gravity M-M forces. */
+ struct link *mm;
+
+ /*! The multipole initialistation task */
+ struct task *init;
+
+ /*! Implicit task for the gravity initialisation */
+ struct task *init_out;
+
+ /*! Task computing long range non-periodic gravity interactions */
+ struct task *long_range;
+
+ /*! Implicit task for the down propagation */
+ struct task *down_in;
+
+ /*! Task propagating the mesh forces to the particles */
+ struct task *mesh;
+
+ /*! Task propagating the multipole to the particles */
+ struct task *down;
+
/*! Minimum end of (integer) time step in this cell for gravity tasks. */
integertime_t ti_end_min;
@@ -383,15 +416,18 @@ struct cell {
/*! Last (integer) time the cell's multipole was drifted forward in time. */
integertime_t ti_old_multipole;
- /*! Nr of #gpart in this cell. */
- int count;
-
/*! Spin lock for various uses (#gpart case). */
swift_lock_type plock;
/*! Spin lock for various uses (#multipole case). */
swift_lock_type mlock;
+ /*! Nr of #gpart in this cell. */
+ int count;
+
+ /*! Nr of #gpart this cell can hold after addition of new #gpart. */
+ int count_total;
+
/*! Number of #gpart updated in this cell. */
int updated;
@@ -404,58 +440,49 @@ struct cell {
/*! Is the #multipole data of this cell being used in a sub-cell? */
int mhold;
+ /*! Number of M-M tasks that are associated with this cell. */
+ short int nr_mm_tasks;
+
/*! Does this cell need to be drifted (gravity)? */
char do_drift;
/*! Do any of this cell's sub-cells need to be drifted (gravity)? */
char do_sub_drift;
- /*! The drift task for gparts */
- struct task *drift;
-
- /*! Implicit task (going up- and down the tree) for the #gpart drifts */
- struct task *drift_out;
-
- /*! Linked list of the tasks computing this cell's gravity forces. */
- struct link *grav;
-
- /*! Linked list of the tasks computing this cell's gravity M-M forces. */
- struct link *mm;
-
- /*! The multipole initialistation task */
- struct task *init;
+ } grav;
- /*! Implicit task for the gravity initialisation */
- struct task *init_out;
+ /*! Stars variables */
+ struct {
- /*! Task computing long range non-periodic gravity interactions */
- struct task *long_range;
+ /*! Pointer to the #spart data. */
+ struct spart *parts;
- /*! Implicit task for the down propagation */
- struct task *down_in;
+ /*! Dependency implicit task for the star ghost (in->ghost->out)*/
+ struct task *ghost_in;
- /*! Task propagating the mesh forces to the particles */
- struct task *mesh;
+ /*! Dependency implicit task for the star ghost (in->ghost->out)*/
+ struct task *ghost_out;
- /*! Task propagating the multipole to the particles */
- struct task *down;
+ /*! The star ghost task itself */
+ struct task *ghost;
- /*! Number of M-M tasks that are associated with this cell. */
- short int nr_mm_tasks;
+ /*! Linked list of the tasks computing this cell's star density. */
+ struct link *density;
- } grav;
+ /*! The task computing this cell's sorts. */
+ struct task *sorts;
- /*! Stars variables */
- struct {
+ /*! Max smoothing length in this cell. */
+ double h_max;
- /*! Pointer to the #spart data. */
- struct spart *parts;
+ /*! Spin lock for various uses (#spart case). */
+ swift_lock_type lock;
/*! Nr of #spart in this cell. */
int count;
- /*! Max smoothing length in this cell. */
- double h_max;
+ /*! Nr of #spart this cell can hold after addition of new #spart. */
+ int count_total;
/*! Values of h_max before the drifts, used for sub-cell tasks. */
float h_max_old;
@@ -490,21 +517,6 @@ struct cell {
/*! Maximum end of (integer) time step in this cell for gravity tasks. */
integertime_t ti_end_min;
- /*! Dependency implicit task for the star ghost (in->ghost->out)*/
- struct task *ghost_in;
-
- /*! Dependency implicit task for the star ghost (in->ghost->out)*/
- struct task *ghost_out;
-
- /*! The star ghost task itself */
- struct task *ghost;
-
- /*! The task computing this cell's sorts. */
- struct task *sorts;
-
- /*! Linked list of the tasks computing this cell's star density. */
- struct link *density;
-
/*! Number of #spart updated in this cell. */
int updated;
@@ -514,9 +526,6 @@ struct cell {
/*! Is the #spart data of this cell being used in a sub-cell? */
int hold;
- /*! Spin lock for various uses (#spart case). */
- swift_lock_type lock;
-
#ifdef SWIFT_DEBUG_CHECKS
/*! Last (integer) time the cell's sort arrays were updated. */
integertime_t ti_sort;
@@ -695,6 +704,8 @@ void cell_activate_hydro_sorts(struct cell *c, int sid, struct scheduler *s);
void cell_activate_stars_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);
+void cell_check_spart_pos(const struct cell *c,
+ const struct spart *global_sparts);
int cell_has_tasks(struct cell *c);
void cell_remove_part(const struct engine *e, struct cell *c, struct part *p,
struct xpart *xp);
@@ -702,10 +713,17 @@ 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);
+struct spart *cell_add_spart(struct engine *e, struct cell *c);
+struct gpart *cell_convert_part_to_gpart(const struct engine *e, struct cell *c,
+ struct part *p, struct xpart *xp);
+struct gpart *cell_convert_spart_to_gpart(const struct engine *e,
+ struct cell *c, struct spart *sp);
+struct spart *cell_convert_part_to_spart(struct engine *e, struct cell *c,
+ struct part *p, struct xpart *xp);
+void cell_reorder_extra_parts(struct cell *c, const ptrdiff_t parts_offset);
+void cell_reorder_extra_gparts(struct cell *c, struct part *parts,
+ struct spart *sparts);
+void cell_reorder_extra_sparts(struct cell *c, const ptrdiff_t sparts_offset);
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,
diff --git a/src/common_io.c b/src/common_io.c
index 087697b489269d97a268966d341093dd666dd9c9..24e74014fd52936023b5c7a41378faf3268bfdb3 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -875,7 +875,8 @@ void io_collect_parts_to_write(const struct part* restrict 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) {
+ if (parts[i].time_bin != time_bin_inhibited &&
+ parts[i].time_bin != time_bin_not_created) {
parts_written[count] = parts[i];
xparts_written[count] = xparts[i];
@@ -909,7 +910,8 @@ void io_collect_sparts_to_write(const struct spart* restrict sparts,
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) {
+ if (sparts[i].time_bin != time_bin_inhibited &&
+ sparts[i].time_bin != time_bin_not_created) {
sparts_written[count] = sparts[i];
count++;
@@ -943,6 +945,7 @@ void io_collect_gparts_to_write(const struct gpart* restrict gparts,
/* And collect the ones that have not been removed */
if ((gparts[i].time_bin != time_bin_inhibited) &&
+ (gparts[i].time_bin != time_bin_not_created) &&
(gparts[i].type == swift_type_dark_matter)) {
gparts_written[count] = gparts[i];
diff --git a/src/engine.c b/src/engine.c
index 73b7d9b64e744e08e19145a574a3a2b896ab3995..9ae8c8273d690b312b53e76dcacd88f20cb8e9b2 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -83,7 +83,6 @@
#include "serial_io.h"
#include "single_io.h"
#include "sort_part.h"
-#include "sourceterms.h"
#include "stars_io.h"
#include "statistics.h"
#include "timers.h"
@@ -110,7 +109,6 @@ const char *engine_policy_names[] = {"none",
"drift everything",
"reconstruct multi-poles",
"cooling",
- "sourceterms",
"stars",
"structure finding",
"star formation",
@@ -1925,8 +1923,8 @@ int engine_estimate_nr_tasks(struct engine *e) {
/* Cooling task + extra space */
n1 += 2;
}
- if (e->policy & engine_policy_sourceterms) {
- n1 += 2;
+ if (e->policy & engine_policy_star_formation) {
+ n1 += 1;
}
if (e->policy & engine_policy_stars) {
/* 2 self (density, feedback), 1 sort, 26/2 density pairs
@@ -2016,9 +2014,10 @@ void engine_rebuild(struct engine *e, int repartitioned,
const ticks tic2 = getticks();
/* Update the global counters of particles */
- long long num_particles[3] = {(long long)e->s->nr_parts,
- (long long)e->s->nr_gparts,
- (long long)e->s->nr_sparts};
+ long long num_particles[3] = {
+ (long long)(e->s->nr_parts - e->s->nr_extra_parts),
+ (long long)(e->s->nr_gparts - e->s->nr_extra_gparts),
+ (long long)(e->s->nr_sparts - e->s->nr_extra_sparts)};
#ifdef WITH_MPI
MPI_Allreduce(MPI_IN_PLACE, num_particles, 3, MPI_LONG_LONG, MPI_SUM,
MPI_COMM_WORLD);
@@ -2588,8 +2587,7 @@ void engine_skip_force_and_kick(struct engine *e) {
t->type == task_type_timestep || t->subtype == task_subtype_force ||
t->subtype == task_subtype_grav || t->type == task_type_end_force ||
t->type == task_type_grav_long_range || t->type == task_type_grav_mm ||
- t->type == task_type_grav_down || t->type == task_type_cooling ||
- t->type == task_type_sourceterms)
+ t->type == task_type_grav_down || t->type == task_type_cooling)
t->skip = 1;
}
@@ -2815,6 +2813,10 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
double *prev_x = s->parts[0].x;
long long *prev_id = &s->parts[0].id;
for (size_t k = 1; k < s->nr_parts; k++) {
+
+ /* Ignore fake buffer particles for on-the-fly creation */
+ if (s->parts[k].time_bin == time_bin_not_created) continue;
+
if (prev_x[0] == s->parts[k].x[0] && prev_x[1] == s->parts[k].x[1] &&
prev_x[2] == s->parts[k].x[2]) {
if (e->verbose)
@@ -2837,6 +2839,10 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
int failed = 0;
double *prev_x = s->gparts[0].x;
for (size_t k = 1; k < s->nr_gparts; k++) {
+
+ /* Ignore fake buffer particles for on-the-fly creation */
+ if (s->gparts[k].time_bin == time_bin_not_created) continue;
+
if (prev_x[0] == s->gparts[k].x[0] && prev_x[1] == s->gparts[k].x[1] &&
prev_x[2] == s->gparts[k].x[2]) {
if (e->verbose)
@@ -3776,8 +3782,8 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
/* Re-allocate the local parts. */
if (e->verbose)
message("Re-allocating parts array from %zu to %zu.", s->size_parts,
- (size_t)(s->nr_parts * 1.2));
- s->size_parts = s->nr_parts * 1.2;
+ (size_t)(s->nr_parts * engine_redistribute_alloc_margin));
+ s->size_parts = s->nr_parts * engine_redistribute_alloc_margin;
struct part *parts_new = NULL;
struct xpart *xparts_new = NULL;
if (posix_memalign((void **)&parts_new, part_align,
@@ -3801,8 +3807,8 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
/* Re-allocate the local sparts. */
if (e->verbose)
message("Re-allocating sparts array from %zu to %zu.", s->size_sparts,
- (size_t)(s->nr_sparts * 1.2));
- s->size_sparts = s->nr_sparts * 1.2;
+ (size_t)(s->nr_sparts * engine_redistribute_alloc_margin));
+ s->size_sparts = s->nr_sparts * engine_redistribute_alloc_margin;
struct spart *sparts_new = NULL;
if (posix_memalign((void **)&sparts_new, spart_align,
sizeof(struct spart) * s->size_sparts) != 0)
@@ -3819,8 +3825,8 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
/* Re-allocate the local gparts. */
if (e->verbose)
message("Re-allocating gparts array from %zu to %zu.", s->size_gparts,
- (size_t)(s->nr_gparts * 1.2));
- s->size_gparts = s->nr_gparts * 1.2;
+ (size_t)(s->nr_gparts * engine_redistribute_alloc_margin));
+ s->size_gparts = s->nr_gparts * engine_redistribute_alloc_margin;
struct gpart *gparts_new = NULL;
if (posix_memalign((void **)&gparts_new, gpart_align,
sizeof(struct gpart) * s->size_gparts) != 0)
@@ -4034,7 +4040,6 @@ void engine_unpin(void) {
* @param potential The properties of the external potential.
* @param cooling_func The properties of the cooling function.
* @param chemistry The chemistry information.
- * @param sourceterms The properties of the source terms function.
*/
void engine_init(struct engine *e, struct space *s, struct swift_params *params,
long long Ngas, long long Ngparts, long long Nstars,
@@ -4046,8 +4051,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
struct pm_mesh *mesh,
const struct external_potential *potential,
struct cooling_function_data *cooling_func,
- const struct chemistry_global_data *chemistry,
- struct sourceterms *sourceterms) {
+ const struct chemistry_global_data *chemistry) {
/* Clean-up everything */
bzero(e, sizeof(struct engine));
@@ -4110,7 +4114,6 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->external_potential = potential;
e->cooling_func = cooling_func;
e->chemistry = chemistry;
- e->sourceterms = sourceterms;
e->parameter_file = params;
e->cell_loc = NULL;
#ifdef WITH_MPI
@@ -5042,9 +5045,11 @@ void engine_recompute_displacement_constraint(struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check that the minimal mass collection worked */
float min_part_mass_check = FLT_MAX;
- for (size_t i = 0; i < e->s->nr_parts; ++i)
+ for (size_t i = 0; i < e->s->nr_parts; ++i) {
+ if (e->s->parts[i].time_bin >= num_time_bins) continue;
min_part_mass_check =
min(min_part_mass_check, hydro_get_mass(&e->s->parts[i]));
+ }
if (min_part_mass_check != min_mass[swift_type_gas])
error("Error collecting minimal mass of gas particles.");
#endif
@@ -5204,7 +5209,6 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
potential_struct_dump(e->external_potential, stream);
cooling_struct_dump(e->cooling_func, stream);
chemistry_struct_dump(e->chemistry, stream);
- sourceterms_struct_dump(e->sourceterms, stream);
parser_struct_dump(e->parameter_file, stream);
if (e->output_list_snapshots)
output_list_struct_dump(e->output_list_snapshots, stream);
@@ -5301,11 +5305,6 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
chemistry_struct_restore(chemistry, stream);
e->chemistry = chemistry;
- struct sourceterms *sourceterms =
- (struct sourceterms *)malloc(sizeof(struct sourceterms));
- sourceterms_struct_restore(sourceterms, stream);
- e->sourceterms = sourceterms;
-
struct swift_params *parameter_file =
(struct swift_params *)malloc(sizeof(struct swift_params));
parser_struct_restore(parameter_file, stream);
diff --git a/src/engine.h b/src/engine.h
index 9c39ea15a85d19cf26dfd94a9c42897251c6b42a..d7d916c1bc25fc0d2cbac26c95de406accec5a67 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -46,7 +46,6 @@
#include "potential.h"
#include "runner.h"
#include "scheduler.h"
-#include "sourceterms_struct.h"
#include "space.h"
#include "task.h"
#include "units.h"
@@ -71,13 +70,12 @@ enum engine_policy {
engine_policy_drift_all = (1 << 11),
engine_policy_reconstruct_mpoles = (1 << 12),
engine_policy_cooling = (1 << 13),
- engine_policy_sourceterms = (1 << 14),
- engine_policy_stars = (1 << 15),
- engine_policy_structure_finding = (1 << 16),
- engine_policy_star_formation = (1 << 17),
- engine_policy_feedback = (1 << 18)
+ engine_policy_stars = (1 << 14),
+ engine_policy_structure_finding = (1 << 15),
+ engine_policy_star_formation = (1 << 16),
+ engine_policy_feedback = (1 << 17)
};
-#define engine_maxpolicy 19
+#define engine_maxpolicy 18
extern const char *engine_policy_names[engine_maxpolicy + 1];
/**
@@ -360,9 +358,6 @@ struct engine {
/* Properties of the chemistry model */
const struct chemistry_global_data *chemistry;
- /* Properties of source terms */
- struct sourceterms *sourceterms;
-
/* The (parsed) parameter file */
struct swift_params *parameter_file;
@@ -417,8 +412,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
struct pm_mesh *mesh,
const struct external_potential *potential,
struct cooling_function_data *cooling_func,
- const struct chemistry_global_data *chemistry,
- struct sourceterms *sourceterms);
+ const struct chemistry_global_data *chemistry);
void engine_config(int restart, struct engine *e, struct swift_params *params,
int nr_nodes, int nodeID, int nr_threads, int with_aff,
int verbose, const char *restart_file);
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index 68841aa5999441e6a2621f867038a44e9f52794c..2c9fe3cb70cfcf118cdb9d1f6fc36a6f854fd7b5 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -666,7 +666,6 @@ void engine_add_ghosts(struct engine *e, struct cell *c, struct task *ghost_in,
void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
struct scheduler *s = &e->sched;
- const int is_with_sourceterms = (e->policy & engine_policy_sourceterms);
/* Are we in a super-cell ? */
if (c->hydro.super == c) {
@@ -696,12 +695,6 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c) {
c->hydro.extra_ghost = scheduler_addtask(
s, task_type_extra_ghost, task_subtype_none, 0, 0, c, NULL);
#endif
-
- /* add source terms */
- if (is_with_sourceterms) {
- c->sourceterms = scheduler_addtask(s, task_type_sourceterms,
- task_subtype_none, 0, 0, c, NULL);
- }
}
} else { /* We are above the super-cell so need to go deeper */
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 746fd4778563aeaab43bafcc7904683ed5b6811c..a3c5e21dbdf8df60b25b01c0326c33c3a10d1bce 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -55,6 +55,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
float wj, wj_dx;
float dv[3], curlvr[3];
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pi->time_bin >= time_bin_inhibited)
+ error("Inhibited pi in interaction function!");
+ if (pj->time_bin >= time_bin_inhibited)
+ error("Inhibited pj in interaction function!");
+#endif
+
/* Get the masses. */
const float mi = pi->mass;
const float mj = pj->mass;
@@ -145,6 +152,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
float wi, wi_dx;
float dv[3], curlvr[3];
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pi->time_bin >= time_bin_inhibited)
+ error("Inhibited pi in interaction function!");
+ if (pj->time_bin >= time_bin_inhibited)
+ error("Inhibited pj in interaction function!");
+#endif
+
/* Get the masses. */
const float mj = pj->mass;
@@ -436,6 +450,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
float wi, wj, wi_dx, wj_dx;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pi->time_bin >= time_bin_inhibited)
+ error("Inhibited pi in interaction function!");
+ if (pj->time_bin >= time_bin_inhibited)
+ error("Inhibited pj in interaction function!");
+#endif
+
/* Cosmological factors entering the EoMs */
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
@@ -558,6 +579,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float wi, wj, wi_dx, wj_dx;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (pi->time_bin >= time_bin_inhibited)
+ error("Inhibited pi in interaction function!");
+ if (pj->time_bin >= time_bin_inhibited)
+ error("Inhibited pj in interaction function!");
+#endif
+
/* Cosmological factors entering the EoMs */
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index e060cb3562f1b319c64d6f6523b18858662312e7..b29f44588c2e13bb5b7c5c9cd5297205557c3fc9 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -54,9 +54,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
float wi, wj, wi_dx, wj_dx;
#ifdef SWIFT_DEBUG_CHECKS
- if (pi->time_bin == time_bin_inhibited)
+ if (pi->time_bin >= time_bin_inhibited)
error("Inhibited pi in interaction function!");
- if (pj->time_bin == time_bin_inhibited)
+ if (pj->time_bin >= time_bin_inhibited)
error("Inhibited pj in interaction function!");
#endif
@@ -135,9 +135,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
float wi, wi_dx;
#ifdef SWIFT_DEBUG_CHECKS
- if (pi->time_bin == time_bin_inhibited)
+ if (pi->time_bin >= time_bin_inhibited)
error("Inhibited pi in interaction function!");
- if (pj->time_bin == time_bin_inhibited)
+ if (pj->time_bin >= time_bin_inhibited)
error("Inhibited pj in interaction function!");
#endif
@@ -196,9 +196,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
struct part *restrict pj, float a, float H) {
#ifdef SWIFT_DEBUG_CHECKS
- if (pi->time_bin == time_bin_inhibited)
+ if (pi->time_bin >= time_bin_inhibited)
error("Inhibited pi in interaction function!");
- if (pj->time_bin == time_bin_inhibited)
+ if (pj->time_bin >= time_bin_inhibited)
error("Inhibited pj in interaction function!");
#endif
@@ -323,9 +323,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const struct part *restrict pj, float a, float H) {
#ifdef SWIFT_DEBUG_CHECKS
- if (pi->time_bin == time_bin_inhibited)
+ if (pi->time_bin >= time_bin_inhibited)
error("Inhibited pi in interaction function!");
- if (pj->time_bin == time_bin_inhibited)
+ if (pj->time_bin >= time_bin_inhibited)
error("Inhibited pj in interaction function!");
#endif
diff --git a/src/parallel_io.c b/src/parallel_io.c
index febbb4a7db6796fd751cb1eacc174a42936d19a2..7fcc22fcae60a393523315cf4cb928629d41541e 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1219,9 +1219,12 @@ void write_output_parallel(struct engine* e, const char* baseName,
// const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
/* Number of particles that we will write */
- const size_t Ntot_written = e->s->nr_gparts - e->s->nr_inhibited_sparts;
- const size_t Ngas_written = e->s->nr_parts - e->s->nr_inhibited_parts;
- const size_t Nstars_written = e->s->nr_sparts - e->s->nr_inhibited_gparts;
+ const size_t Ntot_written =
+ e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
+ const size_t Ngas_written =
+ e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
+ const size_t Nstars_written =
+ e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
const size_t Nbaryons_written = Ngas_written + Nstars_written;
const size_t Ndm_written =
Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
diff --git a/src/part.c b/src/part.c
index 3a626e652cf28f0376cadc1d9a40ab85b752e6c1..ec3627d728f69f469cc7d75eb2beb9ae39ed107e 100644
--- a/src/part.c
+++ b/src/part.c
@@ -139,8 +139,9 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
for (size_t k = 0; k < nr_gparts; ++k) {
- /* We have a DM particle */
- if (gparts[k].type == swift_type_dark_matter) {
+ /* We have a real DM particle */
+ if (gparts[k].type == swift_type_dark_matter &&
+ gparts[k].time_bin != time_bin_not_created) {
/* Check that it's not linked */
if (gparts[k].id_or_neg_offset <= 0)
diff --git a/src/runner.c b/src/runner.c
index f14f8d8bc2721d0edfaca542b68af645bd5ac1a0..018e35c582efb30d4a0f090b478f8a2adb091d39 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -58,7 +58,6 @@
#include "runner_doiact_vec.h"
#include "scheduler.h"
#include "sort_part.h"
-#include "sourceterms.h"
#include "space.h"
#include "space_getsid.h"
#include "stars.h"
@@ -107,42 +106,6 @@
#include "runner_doiact_stars.h"
#undef FUNCTION
-/**
- * @brief Perform source terms
- *
- * @param r runner task
- * @param c cell
- * @param timer 1 if the time is to be recorded.
- */
-void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
- const int count = c->hydro.count;
- const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
- const double cell_width[3] = {c->width[0], c->width[1], c->width[2]};
- struct sourceterms *sourceterms = r->e->sourceterms;
- const int dimen = 3;
-
- TIMER_TIC;
-
- /* Recurse? */
- if (c->split) {
- for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) runner_do_sourceterms(r, c->progeny[k], 0);
- } else {
-
- if (count > 0) {
-
- /* do sourceterms in this cell? */
- const int incell =
- sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
- if (incell == 1) {
- sourceterms_apply(r, sourceterms, c);
- }
- }
- }
-
- if (timer) TIMER_TOC(timer_dosource);
-}
-
/**
* @brief Intermediate task after the density to check that the smoothing
* lengths are correct.
@@ -499,7 +462,7 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
*/
void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
- const struct engine *e = r->e;
+ struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
const int count = c->hydro.count;
struct part *restrict parts = c->hydro.parts;
@@ -529,9 +492,16 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
// MATTHIEU: Temporary star-formation law
// Do not use this at home.
- if (rho > 1.5e7 && e->step > 2) {
+ if (rho > 1.7e7 && e->step > 2) {
message("Removing particle id=%lld rho=%e", p->id, rho);
- cell_convert_part_to_gpart(e, c, p, xp);
+
+ struct spart *sp = cell_convert_part_to_spart(e, c, p, xp);
+
+ /* Did we run out of fresh particles? */
+ if (sp == NULL) continue;
+
+ /* Set everything to a valid state */
+ stars_init_spart(sp);
}
}
}
@@ -623,9 +593,6 @@ void runner_do_sort_ascending(struct entry *sort, int N) {
* @brief Recursively checks that the flags are consistent in a cell hierarchy.
*
* Debugging function. Exists in two flavours: hydro & stars.
- *
- * @param c The #cell to check.
- * @param flags The sorting flags to check.
*/
#define RUNNER_CHECK_SORTS(TYPE) \
void runner_check_sorts_##TYPE(struct cell *c, int flags) { \
@@ -2913,9 +2880,6 @@ void *runner_main(void *data) {
case task_type_star_formation:
runner_do_star_formation(r, t->ci, 1);
break;
- case task_type_sourceterms:
- runner_do_sourceterms(r, t->ci, 1);
- break;
default:
error("Unknown/invalid task type (%d).", t->type);
}
diff --git a/src/serial_io.c b/src/serial_io.c
index 059318df180e0d06e446f9d3f839b16439dd1b34..8ccd2b5fcffd309eb6b4e465dd935922302f4025 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -792,9 +792,12 @@ void write_output_serial(struct engine* e, const char* baseName,
// const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
/* Number of particles that we will write */
- const size_t Ntot_written = e->s->nr_gparts - e->s->nr_inhibited_sparts;
- const size_t Ngas_written = e->s->nr_parts - e->s->nr_inhibited_parts;
- const size_t Nstars_written = e->s->nr_sparts - e->s->nr_inhibited_gparts;
+ const size_t Ntot_written =
+ e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
+ const size_t Ngas_written =
+ e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
+ const size_t Nstars_written =
+ e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
const size_t Nbaryons_written = Ngas_written + Nstars_written;
const size_t Ndm_written =
Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
diff --git a/src/single_io.c b/src/single_io.c
index 833c4a80cb2f43455d34ad3cd694255b7b19038c..e04e4e9ea9fcea67fb64db6b5756d87d692e6d57 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -654,9 +654,12 @@ void write_output_single(struct engine* e, const char* baseName,
// const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
/* Number of particles that we will write */
- const size_t Ntot_written = e->s->nr_gparts - e->s->nr_inhibited_sparts;
- const size_t Ngas_written = e->s->nr_parts - e->s->nr_inhibited_parts;
- const size_t Nstars_written = e->s->nr_sparts - e->s->nr_inhibited_gparts;
+ const size_t Ntot_written =
+ e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
+ const size_t Ngas_written =
+ e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
+ const size_t Nstars_written =
+ e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
const size_t Nbaryons_written = Ngas_written + Nstars_written;
const size_t Ndm_written =
Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
diff --git a/src/sourceterms.c b/src/sourceterms.c
deleted file mode 100644
index 993045e61503e4e78b855816921bc057706b76d1..0000000000000000000000000000000000000000
--- a/src/sourceterms.c
+++ /dev/null
@@ -1,85 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Tom Theuns (tom.theuns@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"
-
-/* Local includes. */
-#include "const.h"
-#include "hydro.h"
-#include "parser.h"
-#include "units.h"
-
-/* This object's header. */
-#include "sourceterms.h"
-
-/**
- * @brief Initialises the sourceterms
- *
- * @param parameter_file The parsed parameter file
- * @param us The current internal system of units
- * @param source the structure that has all the source term properties
- */
-void sourceterms_init(struct swift_params *parameter_file,
- struct unit_system *us, struct sourceterms *source) {
-#ifdef SOURCETERMS_SN_FEEDBACK
- supernova_init(parameter_file, us, source);
-#endif /* SOURCETERMS_SN_FEEDBACK */
-};
-
-/**
- * @brief Prints the properties of the source terms to stdout
- * @param source the structure that has all the source term properties
- */
-void sourceterms_print(struct sourceterms *source) {
-#ifdef SOURCETERMS_NONE
- error(" no sourceterms defined yet you ran with -F");
-#ifdef SOURCETERMS_SN_FEEDBACK
-#error "can't have sourceterms when defined SOURCETERMS_NONE"
-#endif
-#endif
-#ifdef SOURCETERMS_SN_FEEDBACK
- supernova_print(source);
-#endif /* SOURCETERMS_SN_FEEDBACK */
-};
-
-/**
- * @brief Write a sourceterms struct to the given FILE as a stream of bytes.
- *
- * @param sourceterms the struct
- * @param stream the file stream
- */
-void sourceterms_struct_dump(const struct sourceterms *sourceterms,
- FILE *stream) {
- restart_write_blocks((void *)sourceterms, sizeof(struct sourceterms), 1,
- stream, "sourceterms", "sourceterms");
-}
-
-/**
- * @brief Restore a sourceterms struct from the given FILE as a stream of
- * bytes.
- *
- * @param sourceterms the struct
- * @param stream the file stream
- */
-void sourceterms_struct_restore(const struct sourceterms *sourceterms,
- FILE *stream) {
- restart_read_blocks((void *)sourceterms, sizeof(struct sourceterms), 1,
- stream, NULL, "sourceterms");
-}
diff --git a/src/sourceterms.h b/src/sourceterms.h
deleted file mode 100644
index 407d2f19362531a3fd3537889593c484319919b5..0000000000000000000000000000000000000000
--- a/src/sourceterms.h
+++ /dev/null
@@ -1,78 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_SOURCETERMS_H
-#define SWIFT_SOURCETERMS_H
-
-/**
- * @file src/sourceterms.h
- * @brief Branches between the different sourceterms functions.
- */
-
-#include "./const.h"
-#include "runner.h"
-
-#ifdef SOURCETERMS_SN_FEEDBACK
-#include "sourceterms/sn_feedback/sn_feedback_struct.h"
-#endif
-
-/* So far only one model here */
-struct sourceterms {
-#ifdef SOURCETERMS_SN_FEEDBACK
- struct supernova_struct supernova;
-#endif
-};
-#ifdef SOURCETERMS_SN_FEEDBACK
-#include "sourceterms/sn_feedback/sn_feedback.h"
-#endif
-
-void sourceterms_init(struct swift_params* parameter_file,
- struct unit_system* us, struct sourceterms* source);
-void sourceterms_print(struct sourceterms* source);
-
-/* Dump/restore. */
-void sourceterms_struct_dump(const struct sourceterms* source, FILE* stream);
-void sourceterms_struct_restore(const struct sourceterms* source, FILE* stream);
-
-/**
- * @brief Routines related to source terms
- * @param cell_min: corner of cell to test
- * @param cell_width: width of cell to test
- * @param sourceterms: properties of source terms to test
- * @param dimen: dimensionality of the problem
- *
- * This routine tests whether a source term should be applied to this cell
- * return: 1 if yes, return: 0 if no
- */
-
-__attribute__((always_inline)) INLINE static int sourceterms_test_cell(
- const double cell_min[], const double cell_width[],
- struct sourceterms* sourceterms, const int dimen) {
-#ifdef SOURCETERMS_SN_FEEDBACK
- return supernova_feedback_test_cell(cell_min, cell_width, sourceterms, dimen);
-#endif
- return 0;
-};
-
-__attribute__((always_inline)) INLINE static void sourceterms_apply(
- struct runner* r, struct sourceterms* sourceterms, struct cell* c) {
-#ifdef SOURCETERMS_SN_FEEDBACK
- supernova_feedback_apply(r, sourceterms, c);
-#endif
-};
-#endif /* SWIFT_SOURCETERMS_H */
diff --git a/src/sourceterms/sn_feedback/sn_feedback.h b/src/sourceterms/sn_feedback/sn_feedback.h
deleted file mode 100644
index 411673c37e82ff89d906425d1cadaa135c46a38d..0000000000000000000000000000000000000000
--- a/src/sourceterms/sn_feedback/sn_feedback.h
+++ /dev/null
@@ -1,192 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_SN_FEEDBACK_H
-#define SWIFT_SN_FEEDBACK_H
-#include <float.h>
-/* Config parameters. */
-#include "../config.h"
-
-#include "engine.h"
-#include "equation_of_state.h"
-#include "hydro.h"
-#include "runner.h"
-#include "timestep.h"
-
-/**
- * @file src/sourceterms/sn_feedback.h
- *
- * @brief Routines related to sourceterms (supernova feedback): determine if
- * feedback occurs in this cell
- *
- * @param cell_min: corner of cell to test
- * @param cell_width: width of cell to test
- * @param sourceterms: properties of source terms to test
- * @param dimen: dimensionality of the problem
- *
- * This routine tests whether a source term should be applied to this cell
- * return: 1 if yes, return: 0 if no
- */
-__attribute__((always_inline)) INLINE static int supernova_feedback_test_cell(
- const double cell_min[], const double cell_width[],
- struct sourceterms* sourceterms, const int dimen) {
- if (sourceterms->supernova.status == supernova_is_done) return 0;
-
- const double location[3] = {sourceterms->supernova.x,
- sourceterms->supernova.y,
- sourceterms->supernova.z};
- for (int i = 0; i < dimen; i++) {
- if (cell_min[i] > location[i]) return 0;
- if ((cell_min[i] + cell_width[i]) <= location[i]) return 0;
- };
- return 1;
-};
-
-/**
- * @file src/sourceterms/sn_feedback.h
- *
- * @brief Routines related to source terms (supernova feedback): perform
- * feedback in this cell
- * @param r: the runner
- * @param sourceterms the structure describing the source terms properties
- * @param c the cell to apply feedback to
- *
- * This routine heats an individual particle (p), increasing its thermal energy
- * per unit mass
- * by supernova energy / particle mass.
- */
-__attribute__((always_inline)) INLINE static void supernova_feedback_apply(
- struct runner* restrict r, struct sourceterms* restrict sourceterms,
- struct cell* restrict c) {
-
- const int count = c->count;
- struct part* restrict parts = c->parts;
- struct xpart* restrict xparts = c->xparts;
- const double timeBase = r->e->timeBase;
- const int ti_current = r->e->ti_current;
-
- /* inject SN energy into the particle with highest id in this cell if it is
- * active */
- int imax = 0;
- struct part* restrict p_sn = NULL;
- struct xpart* restrict xp_sn = NULL;
-
- for (int i = 0; i < count; i++) {
-
- /* Get a direct pointer on the part. */
- struct part* restrict p = &parts[i];
- if (p->id > imax) {
- imax = p->id;
- p_sn = p;
- xp_sn = &xparts[i];
- }
- }
-
- /* Is this part within the time step? */
- if (p_sn->ti_begin == ti_current) {
-
- /* Does this time step straddle the feedback injection time? */
- const float t_begin = p_sn->ti_begin * timeBase;
- const float t_end = p_sn->ti_end * timeBase;
- if (t_begin <= sourceterms->supernova.time &&
- t_end > sourceterms->supernova.time) {
-
- /* store old time step */
- const int dti_old = p_sn->ti_end - p_sn->ti_begin;
-
- /* add supernova feedback */
- const float u_old = hydro_get_internal_energy(p_sn, 0);
- const float ent_old = hydro_get_entropy(p_sn, 0.0);
- const float u_new =
- u_old + sourceterms->supernova.energy / hydro_get_mass(p_sn);
- hydro_set_internal_energy(p_sn, u_new);
- const float u_set = hydro_get_internal_energy(p_sn, 0.0);
- const float ent_set = hydro_get_entropy(p_sn, 0.0);
- message(
- " applied super nova, time = %e, location= %e %e %e velocity= %e %e "
- "%e",
- ti_current * timeBase, p_sn->x[0], p_sn->x[1], p_sn->x[2], p_sn->v[0],
- p_sn->v[1], p_sn->v[2]);
- message(
- " injected SN energy in particle = %lld, increased energy from %e to "
- "%e and is notw %e, entropy from %e to %e",
- p_sn->id, u_old, u_new, u_set, ent_old, ent_set);
-
- /* label supernova as done */
- sourceterms->supernova.status = supernova_is_done;
-
- /* update timestep if new time step shorter than old time step */
- const int dti = get_part_timestep(p_sn, xp_sn, r->e);
- if (dti < dti_old) {
- p_sn->ti_end = p_sn->ti_begin + dti;
- message(" changed timestep from %d to %d", dti_old, dti);
-
- /* apply simple time-step limiter on all particles in same cell:
- */
- int i_limit = 0;
- for (int i = 0; i < count; i++) {
- struct part* restrict p = &parts[i];
- const int dti_old = p->ti_end - p->ti_begin;
- if (dti_old > 2 * dti) {
- i_limit++;
- const int dti_new = 2 * dti;
- p->ti_end = p->ti_begin + dti_new;
- message(" old step = %d new step = %d", dti_old, dti_new);
- } else
- message(" old step = %d", dti_old);
- }
- message(" count= %d limited timestep of %d particles ", count, i_limit);
- } /* end of limiter */
- error("end");
- }
- }
-};
-
-/**
- * @file src/sourceterms/sn_feedback.h
- *
- * @brief Routine to initialise supernova feedback
- * @param parameterfile: the parse parmeter file
- * @param us: the unit system in use
- * @param sourceterms the structure describing the source terms properties
- *
- * This routine heats an individual particle (p), increasing its thermal energy
- * per unit mass
- * by supernova energy / particle mass.
- */
-
-__attribute__((always_inline)) INLINE static void supernova_init(
- struct swift_params* parameter_file, struct unit_system* us,
- struct sourceterms* source) {
- source->supernova.time = parser_get_param_double(parameter_file, "SN:time");
- source->supernova.energy =
- parser_get_param_double(parameter_file, "SN:energy");
- source->supernova.x = parser_get_param_double(parameter_file, "SN:x");
- source->supernova.y = parser_get_param_double(parameter_file, "SN:y");
- source->supernova.z = parser_get_param_double(parameter_file, "SN:z");
- source->supernova.status = supernova_is_not_done;
-}
-__attribute__((always_inline)) INLINE static void supernova_print(
- struct sourceterms* source) {
- message(
- " Single SNe of energy= %e will explode at time= %e at location "
- "(%e,%e,%e)",
- source->supernova.energy, source->supernova.time, source->supernova.x,
- source->supernova.y, source->supernova.z);
-}
-#endif /* SWIFT_SN_FEEDBACK_H */
diff --git a/src/sourceterms/sn_feedback/sn_feedback_struct.h b/src/sourceterms/sn_feedback/sn_feedback_struct.h
deleted file mode 100644
index dd1842a6717c6c5a20352324cbe6b018c73e7b3e..0000000000000000000000000000000000000000
--- a/src/sourceterms/sn_feedback/sn_feedback_struct.h
+++ /dev/null
@@ -1,45 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Tom Theuns (tom.theuns@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/>.
- *
- ******************************************************************************/
-/**
- * @file src/sourceterms/sn_feedback_struct.h
- * @brief Routines related to source terms (feedback)
- *
- * enumeration type that sets if supernova explosion is done (is_done) or still
- * needs doing (is_not_done)
- */
-#ifndef SWIFT_SN_FEEDBACK_STRUCT_H
-#define SWIFT_SN_FEEDBACK_STRUCT_H
-enum supernova_status { supernova_is_done, supernova_is_not_done };
-
-/**
- * @file src/sourceterms/sn_feedback_struct.h
- * @brief Routines related to source terms (feedback)
- *
- * The structure that describes the source term (supernova feedback)
- * It specifies the time, energy and location of the desired supernova
- * explosion, and a status (supernova_is_done/supernova_is_not_done)
- * that records the status of the supernova
- */
-struct supernova_struct {
- double time;
- double energy;
- double x, y, z;
- enum supernova_status status;
-};
-#endif /* SWIFT_SN_FEEDBACK_STRUCT_H */
diff --git a/src/sourceterms_struct.h b/src/sourceterms_struct.h
deleted file mode 100644
index b3c38986db52d72df825fda97b36c985dff922b6..0000000000000000000000000000000000000000
--- a/src/sourceterms_struct.h
+++ /dev/null
@@ -1,26 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_SOURCETERMS_STRUCT_H
-#define SWIFT_SOURCETERMS_STRUCT_H
-#include "./const.h"
-#ifdef SOURCETERMS_SN_FEEDBACK
-#include "sourceterms/sn_feedback/sn_feedback_struct.h"
-#endif
-
-#endif /* SWIFT_SOURCETERMS_STRUCT_H */
diff --git a/src/space.c b/src/space.c
index 82f369a501bc1be13d27d8096acc8a57a004a580..fa34146d80071fc0c581675ec259b2b139c9e25e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -70,6 +70,18 @@ int space_subsize_pair_stars = space_subsize_pair_stars_default;
int space_subsize_self_stars = space_subsize_self_stars_default;
int space_subdepth_diff_grav = space_subdepth_diff_grav_default;
int space_maxsize = space_maxsize_default;
+
+/*! Number of extra #part we allocate memory for per top-level cell */
+int space_extra_parts = space_extra_parts_default;
+
+/*! Number of extra #spart we allocate memory for per top-level cell */
+int space_extra_sparts = space_extra_sparts_default;
+
+/*! Number of extra #gpart we allocate memory for per top-level cell */
+int space_extra_gparts = space_extra_gparts_default;
+
+/*! Expected maximal number of strays received at a rebuild */
+int space_expected_max_nr_strays = space_expected_max_nr_strays_default;
#ifdef SWIFT_DEBUG_CHECKS
int last_cell_id;
#endif
@@ -104,9 +116,12 @@ struct index_data {
struct space *s;
int *ind;
int *cell_counts;
- int count_inhibited_part;
- int count_inhibited_gpart;
- int count_inhibited_spart;
+ size_t count_inhibited_part;
+ size_t count_inhibited_gpart;
+ size_t count_inhibited_spart;
+ size_t count_extra_part;
+ size_t count_extra_gpart;
+ size_t count_extra_spart;
};
/**
@@ -136,13 +151,13 @@ void space_rebuild_recycle_rec(struct space *s, struct cell *c,
c->progeny[k]->next = *cell_rec_begin;
*cell_rec_begin = c->progeny[k];
- if (s->gravity) {
+ if (s->with_self_gravity) {
c->progeny[k]->grav.multipole->next = *multipole_rec_begin;
*multipole_rec_begin = c->progeny[k]->grav.multipole;
}
if (*cell_rec_end == NULL) *cell_rec_end = *cell_rec_begin;
- if (s->gravity && *multipole_rec_end == NULL)
+ if (s->with_self_gravity && *multipole_rec_end == NULL)
*multipole_rec_end = *multipole_rec_begin;
c->progeny[k]->grav.multipole = NULL;
@@ -182,12 +197,15 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->hydro.sorted = 0;
c->stars.sorted = 0;
c->hydro.count = 0;
+ c->hydro.count_total = 0;
c->hydro.updated = 0;
c->hydro.inhibited = 0;
c->grav.count = 0;
+ c->grav.count_total = 0;
c->grav.updated = 0;
c->grav.inhibited = 0;
c->stars.count = 0;
+ c->stars.count_total = 0;
c->stars.updated = 0;
c->stars.inhibited = 0;
c->grav.init = NULL;
@@ -232,7 +250,8 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
#ifdef SWIFT_DEBUG_CHECKS
c->cellID = 0;
#endif
- if (s->gravity) bzero(c->grav.multipole, sizeof(struct gravity_tensors));
+ if (s->with_self_gravity)
+ bzero(c->grav.multipole, sizeof(struct gravity_tensors));
for (int i = 0; i < 13; i++) {
if (c->hydro.sort[i] != NULL) {
free(c->hydro.sort[i]);
@@ -445,7 +464,7 @@ void space_regrid(struct space *s, int verbose) {
bzero(s->cells_top, s->nr_cells * sizeof(struct cell));
/* Allocate the multipoles for the top-level cells. */
- if (s->gravity) {
+ if (s->with_self_gravity) {
if (posix_memalign((void **)&s->multipoles_top, multipole_align,
s->nr_cells * sizeof(struct gravity_tensors)) != 0)
error("Failed to allocate top-level multipoles.");
@@ -525,7 +544,7 @@ void space_regrid(struct space *s, int verbose) {
c->mpi.grav.recv = NULL;
c->mpi.grav.send = NULL;
#endif // WITH_MPI
- if (s->gravity) c->grav.multipole = &s->multipoles_top[cid];
+ if (s->with_self_gravity) c->grav.multipole = &s->multipoles_top[cid];
#ifdef SWIFT_DEBUG_CHECKS
c->cellID = -last_cell_id;
last_cell_id++;
@@ -602,6 +621,341 @@ void space_regrid(struct space *s, int verbose) {
clocks_getunit());
}
+/**
+ * @brief Allocate memory for the extra particles used for on-the-fly creation.
+ *
+ * This rarely actually allocates memory. Most of the time, we convert
+ * pre-allocated memory inot extra particles.
+ *
+ * This function also sets the extra particles' location to their top-level
+ * cells. They can then be sorted into their correct memory position later on.
+ *
+ * @param s The current #space.
+ * @param verbose Are we talkative?
+ */
+void space_allocate_extras(struct space *s, int verbose) {
+
+ const int local_nodeID = s->e->nodeID;
+
+ /* Anything to do here? (Abort if we don't want extras)*/
+ if (space_extra_parts == 0 && space_extra_gparts == 0 &&
+ space_extra_sparts == 0)
+ return;
+
+ /* The top-level cells */
+ const struct cell *cells = s->cells_top;
+ const double half_cell_width[3] = {0.5 * cells[0].width[0],
+ 0.5 * cells[0].width[1],
+ 0.5 * cells[0].width[2]};
+
+ /* The current number of particles (including spare ones) */
+ size_t nr_parts = s->nr_parts;
+ size_t nr_gparts = s->nr_gparts;
+ size_t nr_sparts = s->nr_sparts;
+
+ /* The current number of actual particles */
+ size_t nr_actual_parts = nr_parts - s->nr_extra_parts;
+ size_t nr_actual_gparts = nr_gparts - s->nr_extra_gparts;
+ size_t nr_actual_sparts = nr_sparts - s->nr_extra_sparts;
+
+ /* The number of particles we allocated memory for (MPI overhead) */
+ size_t size_parts = s->size_parts;
+ size_t size_gparts = s->size_gparts;
+ size_t size_sparts = s->size_sparts;
+
+ int local_cells = 0;
+ for (int i = 0; i < s->nr_cells; ++i)
+ if (s->cells_top[i].nodeID == local_nodeID) local_cells++;
+
+ /* Number of extra particles we want for each type */
+ const size_t expected_num_extra_parts = local_cells * space_extra_parts;
+ const size_t expected_num_extra_gparts = local_cells * space_extra_gparts;
+ const size_t expected_num_extra_sparts = local_cells * space_extra_sparts;
+
+ if (verbose) {
+ message("Currently have %zd/%zd/%zd real particles.", nr_actual_parts,
+ nr_actual_gparts, nr_actual_sparts);
+ message("Currently have %zd/%zd/%zd spaces for extra particles.",
+ s->nr_extra_parts, s->nr_extra_gparts, s->nr_extra_sparts);
+ message("Requesting space for future %zd/%zd/%zd part/gpart/sparts.",
+ expected_num_extra_parts, expected_num_extra_gparts,
+ expected_num_extra_sparts);
+ }
+
+ if (expected_num_extra_parts < s->nr_extra_parts)
+ error("Reduction in top-level cells number not handled.");
+ if (expected_num_extra_gparts < s->nr_extra_gparts)
+ error("Reduction in top-level cells number not handled.");
+ if (expected_num_extra_sparts < s->nr_extra_sparts)
+ error("Reduction in top-level cells number not handled.");
+
+ /* Do we have enough space for the extra gparts (i.e. we haven't used up any)
+ * ? */
+ if (nr_gparts + expected_num_extra_gparts > size_gparts) {
+
+ /* Ok... need to put some more in the game */
+
+ /* Do we need to reallocate? */
+ if (nr_actual_gparts + expected_num_extra_gparts > size_gparts) {
+
+ size_gparts = (nr_actual_gparts + expected_num_extra_gparts) *
+ engine_redistribute_alloc_margin;
+
+ if (verbose)
+ message("Re-allocating gparts array from %zd to %zd", s->size_gparts,
+ size_gparts);
+
+ /* Create more space for parts */
+ struct gpart *gparts_new = NULL;
+ if (posix_memalign((void **)&gparts_new, gpart_align,
+ sizeof(struct gpart) * size_gparts) != 0)
+ error("Failed to allocate new gpart data");
+ const ptrdiff_t delta = gparts_new - s->gparts;
+ memcpy(gparts_new, s->gparts, sizeof(struct gpart) * s->size_gparts);
+ free(s->gparts);
+ s->gparts = gparts_new;
+
+ /* Update the counter */
+ s->size_gparts = size_gparts;
+
+ /* We now need to reset all the part and spart pointers */
+ for (size_t i = 0; i < nr_parts; ++i) {
+ if (s->parts[i].time_bin != time_bin_not_created)
+ s->parts[i].gpart += delta;
+ }
+ for (size_t i = 0; i < nr_sparts; ++i) {
+ if (s->sparts[i].time_bin != time_bin_not_created)
+ s->sparts[i].gpart += delta;
+ }
+ }
+
+ /* Turn some of the allocated spares into particles we can use */
+ for (size_t i = nr_gparts; i < nr_actual_gparts + expected_num_extra_gparts;
+ ++i) {
+ bzero(&s->gparts[i], sizeof(struct gpart));
+ s->gparts[i].time_bin = time_bin_not_created;
+ s->gparts[i].type = swift_type_dark_matter;
+ s->gparts[i].id_or_neg_offset = -1;
+ }
+
+ /* Put the spare particles in their correct cell */
+#ifdef WITH_MPI
+ error("Need to do this correctly over MPI for only the local cells.");
+#endif
+ int count_in_cell = 0, current_cell = 0;
+ size_t count_extra_gparts = 0;
+ for (size_t i = 0; i < nr_actual_gparts + expected_num_extra_gparts; ++i) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (current_cell == s->nr_cells)
+ error("Cell counter beyond the maximal nr. cells.");
+#endif
+
+ if (s->gparts[i].time_bin == time_bin_not_created) {
+
+ /* We want the extra particles to be at the centre of their cell */
+ s->gparts[i].x[0] = cells[current_cell].loc[0] + half_cell_width[0];
+ s->gparts[i].x[1] = cells[current_cell].loc[1] + half_cell_width[1];
+ s->gparts[i].x[2] = cells[current_cell].loc[2] + half_cell_width[2];
+ ++count_in_cell;
+ count_extra_gparts++;
+ }
+
+ /* Once we have reached the number of extra gpart per cell, we move to the
+ * next */
+ if (count_in_cell == space_extra_gparts) {
+ ++current_cell;
+ count_in_cell = 0;
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (count_extra_gparts != expected_num_extra_gparts)
+ error("Constructed the wrong number of extra gparts (%zd vs. %zd)",
+ count_extra_gparts, expected_num_extra_gparts);
+#endif
+
+ /* Update the counters */
+ s->nr_gparts = nr_actual_gparts + expected_num_extra_gparts;
+ s->nr_extra_gparts = expected_num_extra_gparts;
+ }
+
+ /* Do we have enough space for the extra parts (i.e. we haven't used up any) ?
+ */
+ if (expected_num_extra_parts > s->nr_extra_parts) {
+
+ /* Ok... need to put some more in the game */
+
+ /* Do we need to reallocate? */
+ if (nr_actual_parts + expected_num_extra_parts > size_parts) {
+
+ size_parts = (nr_actual_parts + expected_num_extra_parts) *
+ engine_redistribute_alloc_margin;
+
+ if (verbose)
+ message("Re-allocating parts array from %zd to %zd", s->size_parts,
+ size_parts);
+
+ /* Create more space for parts */
+ struct part *parts_new = NULL;
+ if (posix_memalign((void **)&parts_new, part_align,
+ sizeof(struct part) * size_parts) != 0)
+ error("Failed to allocate new part data");
+ memcpy(parts_new, s->parts, sizeof(struct part) * s->size_parts);
+ free(s->parts);
+ s->parts = parts_new;
+
+ /* Same for xparts */
+ struct xpart *xparts_new = NULL;
+ if (posix_memalign((void **)&xparts_new, xpart_align,
+ sizeof(struct xpart) * size_parts) != 0)
+ error("Failed to allocate new xpart data");
+ memcpy(xparts_new, s->xparts, sizeof(struct xpart) * s->size_parts);
+ free(s->xparts);
+ s->xparts = xparts_new;
+
+ /* Update the counter */
+ s->size_parts = size_parts;
+ }
+
+ /* Turn some of the allocated spares into particles we can use */
+ for (size_t i = nr_parts; i < nr_actual_parts + expected_num_extra_parts;
+ ++i) {
+ bzero(&s->parts[i], sizeof(struct part));
+ bzero(&s->xparts[i], sizeof(struct xpart));
+ s->parts[i].time_bin = time_bin_not_created;
+ s->parts[i].id = -1;
+ }
+
+ /* Put the spare particles in their correct cell */
+#ifdef WITH_MPI
+ error("Need to do this correctly over MPI for only the local cells.");
+#endif
+ int count_in_cell = 0, current_cell = 0;
+ size_t count_extra_parts = 0;
+ for (size_t i = 0; i < nr_actual_parts + expected_num_extra_parts; ++i) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (current_cell == s->nr_cells)
+ error("Cell counter beyond the maximal nr. cells.");
+#endif
+
+ if (s->parts[i].time_bin == time_bin_not_created) {
+
+ /* We want the extra particles to be at the centre of their cell */
+ s->parts[i].x[0] = cells[current_cell].loc[0] + half_cell_width[0];
+ s->parts[i].x[1] = cells[current_cell].loc[1] + half_cell_width[1];
+ s->parts[i].x[2] = cells[current_cell].loc[2] + half_cell_width[2];
+ ++count_in_cell;
+ count_extra_parts++;
+ }
+
+ /* Once we have reached the number of extra part per cell, we move to the
+ * next */
+ if (count_in_cell == space_extra_parts) {
+ ++current_cell;
+ count_in_cell = 0;
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (count_extra_parts != expected_num_extra_parts)
+ error("Constructed the wrong number of extra parts (%zd vs. %zd)",
+ count_extra_parts, expected_num_extra_parts);
+#endif
+
+ /* Update the counters */
+ s->nr_parts = nr_actual_parts + expected_num_extra_parts;
+ s->nr_extra_parts = expected_num_extra_parts;
+ }
+
+ /* Do we have enough space for the extra sparts (i.e. we haven't used up any)
+ * ? */
+ if (nr_actual_sparts + expected_num_extra_sparts > nr_sparts) {
+
+ /* Ok... need to put some more in the game */
+
+ /* Do we need to reallocate? */
+ if (nr_actual_sparts + expected_num_extra_sparts > size_sparts) {
+
+ size_sparts = (nr_actual_sparts + expected_num_extra_sparts) *
+ engine_redistribute_alloc_margin;
+
+ if (verbose)
+ message("Re-allocating sparts array from %zd to %zd", s->size_sparts,
+ size_sparts);
+
+ /* Create more space for parts */
+ struct spart *sparts_new = NULL;
+ if (posix_memalign((void **)&sparts_new, spart_align,
+ sizeof(struct spart) * size_sparts) != 0)
+ error("Failed to allocate new spart data");
+ memcpy(sparts_new, s->sparts, sizeof(struct spart) * s->size_sparts);
+ free(s->sparts);
+ s->sparts = sparts_new;
+
+ /* Update the counter */
+ s->size_sparts = size_sparts;
+ }
+
+ /* Turn some of the allocated spares into particles we can use */
+ for (size_t i = nr_sparts; i < nr_actual_sparts + expected_num_extra_sparts;
+ ++i) {
+ bzero(&s->sparts[i], sizeof(struct spart));
+ s->sparts[i].time_bin = time_bin_not_created;
+ s->sparts[i].id = -42;
+ }
+
+ /* Put the spare particles in their correct cell */
+#ifdef WITH_MPI
+ error("Need to do this correctly over MPI for only the local cells.");
+#endif
+ int count_in_cell = 0, current_cell = 0;
+ size_t count_extra_sparts = 0;
+ for (size_t i = 0; i < nr_actual_sparts + expected_num_extra_sparts; ++i) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (current_cell == s->nr_cells)
+ error("Cell counter beyond the maximal nr. cells.");
+#endif
+
+ if (s->sparts[i].time_bin == time_bin_not_created) {
+
+ /* We want the extra particles to be at the centre of their cell */
+ s->sparts[i].x[0] = cells[current_cell].loc[0] + half_cell_width[0];
+ s->sparts[i].x[1] = cells[current_cell].loc[1] + half_cell_width[1];
+ s->sparts[i].x[2] = cells[current_cell].loc[2] + half_cell_width[2];
+ ++count_in_cell;
+ count_extra_sparts++;
+ }
+
+ /* Once we have reached the number of extra spart per cell, we move to the
+ * next */
+ if (count_in_cell == space_extra_sparts) {
+ ++current_cell;
+ count_in_cell = 0;
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (count_extra_sparts != expected_num_extra_sparts)
+ error("Constructed the wrong number of extra sparts (%zd vs. %zd)",
+ count_extra_sparts, expected_num_extra_sparts);
+#endif
+
+ /* Update the counters */
+ s->nr_sparts = nr_actual_sparts + expected_num_extra_sparts;
+ s->nr_extra_sparts = expected_num_extra_sparts;
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Verify that the links are correct */
+ if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0))
+ part_verify_links(s->parts, s->gparts, s->sparts, nr_parts, nr_gparts,
+ nr_sparts, verbose);
+#endif
+}
+
/**
* @brief Re-build the cells as well as the tasks.
*
@@ -622,67 +976,105 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Re-grid if necessary, or just re-set the cell data. */
space_regrid(s, verbose);
+ /* Allocate extra space for particles that will be created */
+ if (s->with_star_formation) space_allocate_extras(s, verbose);
+
+ struct cell *cells_top = s->cells_top;
+ const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
+ const int local_nodeID = s->e->nodeID;
+
+ /* The current number of particles */
size_t nr_parts = s->nr_parts;
size_t nr_gparts = s->nr_gparts;
size_t nr_sparts = s->nr_sparts;
- int count_inhibited_parts = 0;
- int count_inhibited_gparts = 0;
- int count_inhibited_sparts = 0;
- struct cell *restrict cells_top = s->cells_top;
- const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
- /* Run through the particles and get their cell index. Allocates
- an index that is larger than the number of particles to avoid
- re-allocating after shuffling. */
- const size_t ind_size = s->size_parts + 100;
- int *ind = (int *)malloc(sizeof(int) * ind_size);
- if (ind == NULL) error("Failed to allocate temporary particle indices.");
- int *cell_part_counts = (int *)calloc(sizeof(int), s->nr_cells);
- if (cell_part_counts == NULL)
- error("Failed to allocate cell part count buffer.");
- if (s->size_parts > 0)
- space_parts_get_cell_index(s, ind, cell_part_counts, &count_inhibited_parts,
- verbose);
+ /* The number of particles we allocated memory for */
+ size_t size_parts = s->size_parts;
+ size_t size_gparts = s->size_gparts;
+ size_t size_sparts = s->size_sparts;
+
+ /* Counter for the number of inhibited particles found on the node */
+ size_t count_inhibited_parts = 0;
+ size_t count_inhibited_gparts = 0;
+ size_t count_inhibited_sparts = 0;
+
+ /* Counter for the number of extra particles found on the node */
+ size_t count_extra_parts = 0;
+ size_t count_extra_gparts = 0;
+ size_t count_extra_sparts = 0;
+
+ /* Number of particles we expect to have after strays exchange */
+ const size_t h_index_size = size_parts + space_expected_max_nr_strays;
+ const size_t g_index_size = size_gparts + space_expected_max_nr_strays;
+ const size_t s_index_size = size_sparts + space_expected_max_nr_strays;
+
+ /* Allocate arrays to store the indices of the cells where particles
+ belong. We allocate extra space to allow for particles we may
+ receive from other nodes */
+ int *h_index = (int *)malloc(sizeof(int) * h_index_size);
+ int *g_index = (int *)malloc(sizeof(int) * g_index_size);
+ int *s_index = (int *)malloc(sizeof(int) * s_index_size);
+ if (h_index == NULL || g_index == NULL || s_index == NULL)
+ error("Failed to allocate temporary particle indices.");
+
+ /* Allocate counters of particles that will land in each cell */
+ int *cell_part_counts = (int *)malloc(sizeof(int) * s->nr_cells);
+ int *cell_gpart_counts = (int *)malloc(sizeof(int) * s->nr_cells);
+ int *cell_spart_counts = (int *)malloc(sizeof(int) * s->nr_cells);
+ if (cell_part_counts == NULL || cell_gpart_counts == NULL ||
+ cell_spart_counts == NULL)
+ error("Failed to allocate cell particle count buffer.");
+
+ /* Initialise the counters, including buffer space for future particles */
+ for (int i = 0; i < s->nr_cells; ++i) {
+ cell_part_counts[i] = 0;
+ cell_gpart_counts[i] = 0;
+ cell_spart_counts[i] = 0;
+ }
- /* Run through the gravity particles and get their cell index. */
- const size_t gind_size = s->size_gparts + 100;
- int *gind = (int *)malloc(sizeof(int) * gind_size);
- if (gind == NULL) error("Failed to allocate temporary g-particle indices.");
- int *cell_gpart_counts = (int *)calloc(sizeof(int), s->nr_cells);
- if (cell_gpart_counts == NULL)
- error("Failed to allocate cell gpart count buffer.");
- if (s->size_gparts > 0)
- space_gparts_get_cell_index(s, gind, cell_gpart_counts,
- &count_inhibited_gparts, verbose);
-
- /* Run through the star particles and get their cell index. */
- const size_t sind_size = s->size_sparts + 100;
- int *sind = (int *)malloc(sizeof(int) * sind_size);
- if (sind == NULL) error("Failed to allocate temporary s-particle indices.");
- int *cell_spart_counts = (int *)calloc(sizeof(int), s->nr_cells);
- if (cell_spart_counts == NULL)
- error("Failed to allocate cell gpart count buffer.");
- if (s->size_sparts > 0)
- space_sparts_get_cell_index(s, sind, cell_spart_counts,
- &count_inhibited_sparts, verbose);
+ /* Run through the particles and get their cell index. */
+ if (nr_parts > 0)
+ space_parts_get_cell_index(s, h_index, cell_part_counts,
+ &count_inhibited_parts, &count_extra_parts,
+ verbose);
+ if (nr_gparts > 0)
+ space_gparts_get_cell_index(s, g_index, cell_gpart_counts,
+ &count_inhibited_gparts, &count_extra_gparts,
+ verbose);
+ if (nr_sparts > 0)
+ space_sparts_get_cell_index(s, s_index, cell_spart_counts,
+ &count_inhibited_sparts, &count_extra_sparts,
+ verbose);
#ifdef SWIFT_DEBUG_CHECKS
+ /* Some safety checks */
if (repartitioned && count_inhibited_parts)
error("We just repartitioned but still found inhibited parts.");
if (repartitioned && count_inhibited_sparts)
error("We just repartitioned but still found inhibited sparts.");
if (repartitioned && count_inhibited_gparts)
error("We just repartitioned but still found inhibited gparts.");
-#endif
- const int local_nodeID = s->e->nodeID;
+ if (count_extra_parts != s->nr_extra_parts)
+ error(
+ "Number of extra parts in the part array not matching the space "
+ "counter.");
+ if (count_extra_gparts != s->nr_extra_gparts)
+ error(
+ "Number of extra gparts in the gpart array not matching the space "
+ "counter.");
+ if (count_extra_sparts != s->nr_extra_sparts)
+ error(
+ "Number of extra sparts in the spart array not matching the space "
+ "counter.");
+#endif
/* Move non-local parts and inhibited parts to the end of the list. */
if (!repartitioned && (s->e->nr_nodes > 1 || count_inhibited_parts > 0)) {
for (size_t k = 0; k < nr_parts; /* void */) {
/* Inhibited particle or foreign particle */
- if (ind[k] == -1 || cells_top[ind[k]].nodeID != local_nodeID) {
+ if (h_index[k] == -1 || cells_top[h_index[k]].nodeID != local_nodeID) {
/* One fewer particle */
nr_parts -= 1;
@@ -701,7 +1093,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Swap the xpart */
memswap(&s->xparts[k], &s->xparts[nr_parts], sizeof(struct xpart));
/* Swap the index */
- memswap(&ind[k], &ind[nr_parts], sizeof(int));
+ memswap(&h_index[k], &h_index[nr_parts], sizeof(int));
} else {
/* Increment when not exchanging otherwise we need to retest "k".*/
@@ -712,17 +1104,17 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all parts are in the correct places. */
- int check_count_inhibited_part = 0;
+ size_t check_count_inhibited_part = 0;
for (size_t k = 0; k < nr_parts; k++) {
- if (ind[k] == -1 || cells_top[ind[k]].nodeID != local_nodeID) {
+ if (h_index[k] == -1 || cells_top[h_index[k]].nodeID != local_nodeID) {
error("Failed to move all non-local parts to send list");
}
}
for (size_t k = nr_parts; k < s->nr_parts; k++) {
- if (ind[k] != -1 && cells_top[ind[k]].nodeID == local_nodeID) {
+ if (h_index[k] != -1 && cells_top[h_index[k]].nodeID == local_nodeID) {
error("Failed to remove local parts from send list");
}
- if (ind[k] == -1) ++check_count_inhibited_part;
+ if (h_index[k] == -1) ++check_count_inhibited_part;
}
if (check_count_inhibited_part != count_inhibited_parts)
error("Counts of inhibited particles do not match!");
@@ -733,7 +1125,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
for (size_t k = 0; k < nr_sparts; /* void */) {
/* Inhibited particle or foreign particle */
- if (sind[k] == -1 || cells_top[sind[k]].nodeID != local_nodeID) {
+ if (s_index[k] == -1 || cells_top[s_index[k]].nodeID != local_nodeID) {
/* One fewer particle */
nr_sparts -= 1;
@@ -750,7 +1142,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
/* Swap the index */
- memswap(&sind[k], &sind[nr_sparts], sizeof(int));
+ memswap(&s_index[k], &s_index[nr_sparts], sizeof(int));
} else {
/* Increment when not exchanging otherwise we need to retest "k".*/
@@ -761,17 +1153,17 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all sparts are in the correct place. */
- int check_count_inhibited_spart = 0;
+ size_t check_count_inhibited_spart = 0;
for (size_t k = 0; k < nr_sparts; k++) {
- if (sind[k] == -1 || cells_top[sind[k]].nodeID != local_nodeID) {
+ if (s_index[k] == -1 || cells_top[s_index[k]].nodeID != local_nodeID) {
error("Failed to move all non-local sparts to send list");
}
}
for (size_t k = nr_sparts; k < s->nr_sparts; k++) {
- if (sind[k] != -1 && cells_top[sind[k]].nodeID == local_nodeID) {
+ if (s_index[k] != -1 && cells_top[s_index[k]].nodeID == local_nodeID) {
error("Failed to remove local sparts from send list");
}
- if (sind[k] == -1) ++check_count_inhibited_spart;
+ if (s_index[k] == -1) ++check_count_inhibited_spart;
}
if (check_count_inhibited_spart != count_inhibited_sparts)
error("Counts of inhibited s-particles do not match!");
@@ -782,7 +1174,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
for (size_t k = 0; k < nr_gparts; /* void */) {
/* Inhibited particle or foreign particle */
- if (gind[k] == -1 || cells_top[gind[k]].nodeID != local_nodeID) {
+ if (g_index[k] == -1 || cells_top[g_index[k]].nodeID != local_nodeID) {
/* One fewer particle */
nr_gparts -= 1;
@@ -805,7 +1197,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
/* Swap the index */
- memswap(&gind[k], &gind[nr_gparts], sizeof(int));
+ memswap(&g_index[k], &g_index[nr_gparts], sizeof(int));
} else {
/* Increment when not exchanging otherwise we need to retest "k".*/
k++;
@@ -815,17 +1207,17 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all gparts are in the correct place. */
- int check_count_inhibited_gpart = 0;
+ size_t check_count_inhibited_gpart = 0;
for (size_t k = 0; k < nr_gparts; k++) {
- if (gind[k] == -1 || cells_top[gind[k]].nodeID != local_nodeID) {
+ if (g_index[k] == -1 || cells_top[g_index[k]].nodeID != local_nodeID) {
error("Failed to move all non-local gparts to send list");
}
}
for (size_t k = nr_gparts; k < s->nr_gparts; k++) {
- if (gind[k] != -1 && cells_top[gind[k]].nodeID == local_nodeID) {
+ if (g_index[k] != -1 && cells_top[g_index[k]].nodeID == local_nodeID) {
error("Failed to remove local gparts from send list");
}
- if (gind[k] == -1) ++check_count_inhibited_gpart;
+ if (g_index[k] == -1) ++check_count_inhibited_gpart;
}
if (check_count_inhibited_gpart != count_inhibited_gparts)
error("Counts of inhibited g-particles do not match!");
@@ -834,21 +1226,23 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
#ifdef WITH_MPI
/* Exchange the strays, note that this potentially re-allocates
- the parts arrays. This can be skipped if we just repartitioned aspace
- there should be no strays */
+ the parts arrays. This can be skipped if we just repartitioned space
+ as there should be no strays in that case */
if (!repartitioned) {
size_t nr_parts_exchanged = s->nr_parts - nr_parts;
size_t nr_gparts_exchanged = s->nr_gparts - nr_gparts;
size_t nr_sparts_exchanged = s->nr_sparts - nr_sparts;
- engine_exchange_strays(s->e, nr_parts, &ind[nr_parts], &nr_parts_exchanged,
- nr_gparts, &gind[nr_gparts], &nr_gparts_exchanged,
- nr_sparts, &sind[nr_sparts], &nr_sparts_exchanged);
+ engine_exchange_strays(s->e, nr_parts, &h_index[nr_parts],
+ &nr_parts_exchanged, nr_gparts, &g_index[nr_gparts],
+ &nr_gparts_exchanged, nr_sparts, &s_index[nr_sparts],
+ &nr_sparts_exchanged);
/* Set the new particle counts. */
s->nr_parts = nr_parts + nr_parts_exchanged;
s->nr_gparts = nr_gparts + nr_gparts_exchanged;
s->nr_sparts = nr_sparts + nr_sparts_exchanged;
+
} else {
#ifdef SWIFT_DEBUG_CHECKS
if (s->nr_parts != nr_parts)
@@ -870,23 +1264,23 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
/* Re-allocate the index array for the parts if needed.. */
- if (s->nr_parts + 1 > ind_size) {
+ if (s->nr_parts + 1 > h_index_size) {
int *ind_new;
if ((ind_new = (int *)malloc(sizeof(int) * (s->nr_parts + 1))) == NULL)
error("Failed to allocate temporary particle indices.");
- memcpy(ind_new, ind, sizeof(int) * nr_parts);
- free(ind);
- ind = ind_new;
+ memcpy(ind_new, h_index, sizeof(int) * nr_parts);
+ free(h_index);
+ h_index = ind_new;
}
/* Re-allocate the index array for the sparts if needed.. */
- if (s->nr_sparts + 1 > sind_size) {
+ if (s->nr_sparts + 1 > s_index_size) {
int *sind_new;
if ((sind_new = (int *)malloc(sizeof(int) * (s->nr_sparts + 1))) == NULL)
error("Failed to allocate temporary s-particle indices.");
- memcpy(sind_new, sind, sizeof(int) * nr_sparts);
- free(sind);
- sind = sind_new;
+ memcpy(sind_new, s_index, sizeof(int) * nr_sparts);
+ free(s_index);
+ s_index = sind_new;
}
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
@@ -895,13 +1289,13 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Assign each received part to its cell. */
for (size_t k = nr_parts; k < s->nr_parts; k++) {
const struct part *const p = &s->parts[k];
- ind[k] =
+ h_index[k] =
cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
- cell_part_counts[ind[k]]++;
+ cell_part_counts[h_index[k]]++;
#ifdef SWIFT_DEBUG_CHECKS
- if (cells_top[ind[k]].nodeID != local_nodeID)
+ if (cells_top[h_index[k]].nodeID != local_nodeID)
error("Received part that does not belong to me (nodeID=%i).",
- cells_top[ind[k]].nodeID);
+ cells_top[h_index[k]].nodeID);
#endif
}
nr_parts = s->nr_parts;
@@ -909,13 +1303,13 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Assign each received spart to its cell. */
for (size_t k = nr_sparts; k < s->nr_sparts; k++) {
const struct spart *const sp = &s->sparts[k];
- sind[k] =
+ s_index[k] =
cell_getid(cdim, sp->x[0] * ih[0], sp->x[1] * ih[1], sp->x[2] * ih[2]);
- cell_spart_counts[sind[k]]++;
+ cell_spart_counts[s_index[k]]++;
#ifdef SWIFT_DEBUG_CHECKS
- if (cells_top[sind[k]].nodeID != local_nodeID)
+ if (cells_top[s_index[k]].nodeID != local_nodeID)
error("Received s-part that does not belong to me (nodeID=%i).",
- cells_top[sind[k]].nodeID);
+ cells_top[s_index[k]].nodeID);
#endif
}
nr_sparts = s->nr_sparts;
@@ -930,8 +1324,8 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Sort the parts according to their cells. */
if (nr_parts > 0)
- space_parts_sort(s->parts, s->xparts, ind, cell_part_counts, s->nr_cells,
- 0);
+ space_parts_sort(s->parts, s->xparts, h_index, cell_part_counts,
+ s->nr_cells, 0);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the part have been sorted correctly. */
@@ -949,7 +1343,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* New cell of this part */
const struct cell *c = &s->cells_top[new_ind];
- if (ind[k] != new_ind)
+ if (h_index[k] != new_ind)
error("part's new cell index not matching sorted index.");
if (p->x[0] < c->loc[0] || p->x[0] > c->loc[0] + c->width[0] ||
@@ -961,7 +1355,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Sort the sparts according to their cells. */
if (nr_sparts > 0)
- space_sparts_sort(s->sparts, sind, cell_spart_counts, s->nr_cells, 0);
+ space_sparts_sort(s->sparts, s_index, cell_spart_counts, s->nr_cells, 0);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the spart have been sorted correctly. */
@@ -979,7 +1373,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* New cell of this spart */
const struct cell *c = &s->cells_top[new_sind];
- if (sind[k] != new_sind)
+ if (s_index[k] != new_sind)
error("spart's new cell index not matching sorted index.");
if (sp->x[0] < c->loc[0] || sp->x[0] > c->loc[0] + c->width[0] ||
@@ -989,54 +1383,58 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
#endif /* SWIFT_DEBUG_CHECKS */
- /* Extract the cell counts from the sorted indices. */
+ /* Extract the cell counts from the sorted indices. Deduct the extra
+ * particles. */
size_t last_index = 0;
- ind[nr_parts] = s->nr_cells; // sentinel.
+ h_index[nr_parts] = s->nr_cells; // sentinel.
for (size_t k = 0; k < nr_parts; k++) {
- if (ind[k] < ind[k + 1]) {
- cells_top[ind[k]].hydro.count = k - last_index + 1;
+ if (h_index[k] < h_index[k + 1]) {
+ cells_top[h_index[k]].hydro.count =
+ k - last_index + 1 - space_extra_parts;
last_index = k + 1;
}
}
- /* Extract the cell counts from the sorted indices. */
+ /* Extract the cell counts from the sorted indices. Deduct the extra
+ * particles. */
size_t last_sindex = 0;
- sind[nr_sparts] = s->nr_cells; // sentinel.
+ s_index[nr_sparts] = s->nr_cells; // sentinel.
for (size_t k = 0; k < nr_sparts; k++) {
- if (sind[k] < sind[k + 1]) {
- cells_top[sind[k]].stars.count = k - last_sindex + 1;
+ if (s_index[k] < s_index[k + 1]) {
+ cells_top[s_index[k]].stars.count =
+ k - last_sindex + 1 - space_extra_sparts;
last_sindex = k + 1;
}
}
/* We no longer need the indices as of here. */
- free(ind);
+ free(h_index);
free(cell_part_counts);
- free(sind);
+ free(s_index);
free(cell_spart_counts);
#ifdef WITH_MPI
/* Re-allocate the index array for the gparts if needed.. */
- if (s->nr_gparts + 1 > gind_size) {
+ if (s->nr_gparts + 1 > g_index_size) {
int *gind_new;
if ((gind_new = (int *)malloc(sizeof(int) * (s->nr_gparts + 1))) == NULL)
error("Failed to allocate temporary g-particle indices.");
- memcpy(gind_new, gind, sizeof(int) * nr_gparts);
- free(gind);
- gind = gind_new;
+ memcpy(gind_new, g_index, sizeof(int) * nr_gparts);
+ free(g_index);
+ g_index = gind_new;
}
/* Assign each received gpart to its cell. */
for (size_t k = nr_gparts; k < s->nr_gparts; k++) {
const struct gpart *const p = &s->gparts[k];
- gind[k] =
+ g_index[k] =
cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
- cell_gpart_counts[gind[k]]++;
+ cell_gpart_counts[g_index[k]]++;
#ifdef SWIFT_DEBUG_CHECKS
- if (cells_top[gind[k]].nodeID != s->e->nodeID)
+ if (cells_top[g_index[k]].nodeID != s->e->nodeID)
error("Received g-part that does not belong to me (nodeID=%i).",
- cells_top[gind[k]].nodeID);
+ cells_top[g_index[k]].nodeID);
#endif
}
nr_gparts = s->nr_gparts;
@@ -1055,8 +1453,8 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Sort the gparts according to their cells. */
if (nr_gparts > 0)
- space_gparts_sort(s->gparts, s->parts, s->sparts, gind, cell_gpart_counts,
- s->nr_cells);
+ space_gparts_sort(s->gparts, s->parts, s->sparts, g_index,
+ cell_gpart_counts, s->nr_cells);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the gpart have been sorted correctly. */
@@ -1074,7 +1472,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* New cell of this gpart */
const struct cell *c = &s->cells_top[new_gind];
- if (gind[k] != new_gind)
+ if (g_index[k] != new_gind)
error("gpart's new cell index not matching sorted index.");
if (gp->x[0] < c->loc[0] || gp->x[0] > c->loc[0] + c->width[0] ||
@@ -1084,18 +1482,20 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
}
#endif /* SWIFT_DEBUG_CHECKS */
- /* Extract the cell counts from the sorted indices. */
+ /* Extract the cell counts from the sorted indices. Deduct the extra
+ * particles. */
size_t last_gindex = 0;
- gind[nr_gparts] = s->nr_cells;
+ g_index[nr_gparts] = s->nr_cells;
for (size_t k = 0; k < nr_gparts; k++) {
- if (gind[k] < gind[k + 1]) {
- cells_top[gind[k]].grav.count = k - last_gindex + 1;
+ if (g_index[k] < g_index[k + 1]) {
+ cells_top[g_index[k]].grav.count =
+ k - last_gindex + 1 - space_extra_gparts;
last_gindex = k + 1;
}
}
/* We no longer need the indices as of here. */
- free(gind);
+ free(g_index);
free(cell_gpart_counts);
#ifdef SWIFT_DEBUG_CHECKS
@@ -1134,10 +1534,15 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
c->hydro.xparts = xfinger;
c->grav.parts = gfinger;
c->stars.parts = sfinger;
- finger = &finger[c->hydro.count];
- xfinger = &xfinger[c->hydro.count];
- gfinger = &gfinger[c->grav.count];
- sfinger = &sfinger[c->stars.count];
+
+ c->hydro.count_total = c->hydro.count + space_extra_parts;
+ c->grav.count_total = c->grav.count + space_extra_gparts;
+ c->stars.count_total = c->stars.count + space_extra_sparts;
+
+ finger = &finger[c->hydro.count_total];
+ xfinger = &xfinger[c->hydro.count_total];
+ gfinger = &gfinger[c->grav.count_total];
+ sfinger = &sfinger[c->stars.count_total];
/* Add this cell to the list of local cells */
s->local_cells_top[s->nr_local_cells] = k;
@@ -1160,13 +1565,17 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
clocks_from_ticks(getticks() - tic2), clocks_getunit());
}
- /* At this point, we have the upper-level cells, old or new. Now make
- sure that the parts in each cell are ok. */
+ /* Re-order the extra particles such that they are at the end of their cell's
+ memory pool. */
+ if (s->with_star_formation) space_reorder_extras(s, verbose);
+
+ /* At this point, we have the upper-level cells. Now recursively split each
+ cell to get the full AMR grid. */
space_split(s, verbose);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that the multipole construction went OK */
- if (s->gravity)
+ if (s->with_self_gravity)
for (int k = 0; k < s->nr_cells; k++)
cell_check_multipole(&s->cells_top[k]);
#endif
@@ -1201,6 +1610,75 @@ void space_split(struct space *s, int verbose) {
clocks_getunit());
}
+void space_reorder_extra_parts_mapper(void *map_data, int num_cells,
+ void *extra_data) {
+
+ struct cell *cells_top = (struct cell *)map_data;
+ struct space *s = (struct space *)extra_data;
+
+ for (int ind = 0; ind < num_cells; ind++) {
+ struct cell *c = &cells_top[ind];
+ cell_reorder_extra_parts(c, c->hydro.parts - s->parts);
+ }
+}
+
+void space_reorder_extra_gparts_mapper(void *map_data, int num_cells,
+ void *extra_data) {
+
+ struct cell *cells_top = (struct cell *)map_data;
+ struct space *s = (struct space *)extra_data;
+
+ for (int ind = 0; ind < num_cells; ind++) {
+ struct cell *c = &cells_top[ind];
+ cell_reorder_extra_gparts(c, s->parts, s->sparts);
+ }
+}
+
+void space_reorder_extra_sparts_mapper(void *map_data, int num_cells,
+ void *extra_data) {
+
+ struct cell *cells_top = (struct cell *)map_data;
+ struct space *s = (struct space *)extra_data;
+
+ for (int ind = 0; ind < num_cells; ind++) {
+ struct cell *c = &cells_top[ind];
+ cell_reorder_extra_sparts(c, c->stars.parts - s->sparts);
+ }
+}
+
+/**
+ * @brief Re-orders the particles in each cell such that the extra particles
+ * for on-the-fly creation are located at the end of their respective cells.
+ *
+ * This assumes that all the particles (real and extra) have already been sorted
+ * in their correct top-level cell.
+ *
+ * @param s The #space to act upon.
+ * @param verbose Are we talkative?
+ */
+void space_reorder_extras(struct space *s, int verbose) {
+
+#ifdef WITH_MPI
+ if (space_extra_parts || space_extra_gparts || space_extra_sparts)
+ error("Need an MPI-proof version of this.");
+#endif
+
+ /* Re-order the gas particles */
+ if (space_extra_parts)
+ threadpool_map(&s->e->threadpool, space_reorder_extra_parts_mapper,
+ s->cells_top, s->nr_cells, sizeof(struct cell), 0, s);
+
+ /* Re-order the gravity particles */
+ if (space_extra_gparts)
+ threadpool_map(&s->e->threadpool, space_reorder_extra_gparts_mapper,
+ s->cells_top, s->nr_cells, sizeof(struct cell), 0, s);
+
+ /* Re-order the star particles */
+ if (space_extra_sparts)
+ threadpool_map(&s->e->threadpool, space_reorder_extra_sparts_mapper,
+ s->cells_top, s->nr_cells, sizeof(struct cell), 0, s);
+}
+
/**
* @brief #threadpool mapper function to sanitize the cells
*
@@ -1264,7 +1742,8 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts,
/* Init the local collectors */
float min_mass = FLT_MAX;
float sum_vel_norm = 0.f;
- int count_inhibited_part = 0;
+ size_t count_inhibited_part = 0;
+ size_t count_extra_part = 0;
/* Loop over the parts. */
for (int k = 0; k < nr_parts; k++) {
@@ -1296,26 +1775,31 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts,
pos_z);
#endif
- /* Is this particle to be removed? */
if (p->time_bin == time_bin_inhibited) {
+ /* Is this particle to be removed? */
ind[k] = -1;
++count_inhibited_part;
+ } else if (p->time_bin == time_bin_not_created) {
+ /* Is this a place-holder for on-the-fly creation? */
+ ind[k] = index;
+ cell_counts[index]++;
+ ++count_extra_part;
} else {
- /* List its top-level cell index */
+ /* Normal case: list its top-level cell index */
ind[k] = index;
cell_counts[index]++;
- }
- /* Compute minimal mass */
- min_mass = min(min_mass, hydro_get_mass(p));
+ /* Compute minimal mass */
+ min_mass = min(min_mass, hydro_get_mass(p));
- /* Compute sum of velocity norm */
- sum_vel_norm += p->v[0] * p->v[0] + p->v[1] * p->v[1] + p->v[2] * p->v[2];
+ /* Compute sum of velocity norm */
+ sum_vel_norm += p->v[0] * p->v[0] + p->v[1] * p->v[1] + p->v[2] * p->v[2];
- /* Update the position */
- p->x[0] = pos_x;
- p->x[1] = pos_y;
- p->x[2] = pos_z;
+ /* Update the position */
+ p->x[0] = pos_x;
+ p->x[1] = pos_y;
+ p->x[2] = pos_z;
+ }
}
/* Write the counts back to the global array. */
@@ -1323,8 +1807,10 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts,
if (cell_counts[k]) atomic_add(&data->cell_counts[k], cell_counts[k]);
free(cell_counts);
- /* Write the count of inhibited parts */
- atomic_add(&data->count_inhibited_part, count_inhibited_part);
+ /* Write the count of inhibited and extra parts */
+ if (count_inhibited_part)
+ atomic_add(&data->count_inhibited_part, count_inhibited_part);
+ if (count_extra_part) atomic_add(&data->count_extra_part, count_extra_part);
/* Write back the minimal part mass and velocity sum */
atomic_min_f(&s->min_part_mass, min_mass);
@@ -1364,7 +1850,8 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
/* Init the local collectors */
float min_mass = FLT_MAX;
float sum_vel_norm = 0.f;
- int count_inhibited_gpart = 0;
+ size_t count_inhibited_gpart = 0;
+ size_t count_extra_gpart = 0;
for (int k = 0; k < nr_gparts; k++) {
@@ -1395,31 +1882,36 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
pos_z);
#endif
- /* Is this particle to be removed? */
if (gp->time_bin == time_bin_inhibited) {
+ /* Is this particle to be removed? */
ind[k] = -1;
++count_inhibited_gpart;
+ } else if (gp->time_bin == time_bin_not_created) {
+ /* Is this a place-holder for on-the-fly creation? */
+ ind[k] = index;
+ cell_counts[index]++;
+ ++count_extra_gpart;
} else {
/* List its top-level cell index */
ind[k] = index;
cell_counts[index]++;
- }
- if (gp->type == swift_type_dark_matter) {
+ if (gp->type == swift_type_dark_matter) {
- /* Compute minimal mass */
- min_mass = min(min_mass, gp->mass);
+ /* Compute minimal mass */
+ min_mass = min(min_mass, gp->mass);
- /* Compute sum of velocity norm */
- sum_vel_norm += gp->v_full[0] * gp->v_full[0] +
- gp->v_full[1] * gp->v_full[1] +
- gp->v_full[2] * gp->v_full[2];
- }
+ /* Compute sum of velocity norm */
+ sum_vel_norm += gp->v_full[0] * gp->v_full[0] +
+ gp->v_full[1] * gp->v_full[1] +
+ gp->v_full[2] * gp->v_full[2];
+ }
- /* Update the position */
- gp->x[0] = pos_x;
- gp->x[1] = pos_y;
- gp->x[2] = pos_z;
+ /* Update the position */
+ gp->x[0] = pos_x;
+ gp->x[1] = pos_y;
+ gp->x[2] = pos_z;
+ }
}
/* Write the counts back to the global array. */
@@ -1427,8 +1919,11 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
if (cell_counts[k]) atomic_add(&data->cell_counts[k], cell_counts[k]);
free(cell_counts);
- /* Write the count of inhibited gparts */
- atomic_add(&data->count_inhibited_gpart, count_inhibited_gpart);
+ /* Write the count of inhibited and extra gparts */
+ if (count_inhibited_gpart)
+ atomic_add(&data->count_inhibited_gpart, count_inhibited_gpart);
+ if (count_extra_gpart)
+ atomic_add(&data->count_extra_gpart, count_extra_gpart);
/* Write back the minimal part mass and velocity sum */
atomic_min_f(&s->min_gpart_mass, min_mass);
@@ -1468,7 +1963,8 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
/* Init the local collectors */
float min_mass = FLT_MAX;
float sum_vel_norm = 0.f;
- int count_inhibited_spart = 0;
+ size_t count_inhibited_spart = 0;
+ size_t count_extra_spart = 0;
for (int k = 0; k < nr_sparts; k++) {
@@ -1503,23 +1999,28 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
if (sp->time_bin == time_bin_inhibited) {
ind[k] = -1;
++count_inhibited_spart;
+ } else if (sp->time_bin == time_bin_not_created) {
+ /* Is this a place-holder for on-the-fly creation? */
+ ind[k] = index;
+ cell_counts[index]++;
+ ++count_extra_spart;
} else {
/* List its top-level cell index */
ind[k] = index;
cell_counts[index]++;
- }
- /* Compute minimal mass */
- min_mass = min(min_mass, sp->mass);
+ /* Compute minimal mass */
+ min_mass = min(min_mass, sp->mass);
- /* Compute sum of velocity norm */
- sum_vel_norm +=
- sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2];
+ /* Compute sum of velocity norm */
+ sum_vel_norm +=
+ sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2];
- /* Update the position */
- sp->x[0] = pos_x;
- sp->x[1] = pos_y;
- sp->x[2] = pos_z;
+ /* Update the position */
+ sp->x[0] = pos_x;
+ sp->x[1] = pos_y;
+ sp->x[2] = pos_z;
+ }
}
/* Write the counts back to the global array. */
@@ -1527,8 +2028,11 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
if (cell_counts[k]) atomic_add(&data->cell_counts[k], cell_counts[k]);
free(cell_counts);
- /* Write the count of inhibited parts */
- atomic_add(&data->count_inhibited_spart, count_inhibited_spart);
+ /* Write the count of inhibited and extra sparts */
+ if (count_inhibited_spart)
+ atomic_add(&data->count_inhibited_spart, count_inhibited_spart);
+ if (count_extra_spart)
+ atomic_add(&data->count_extra_spart, count_extra_spart);
/* Write back the minimal part mass and velocity sum */
atomic_min_f(&s->min_spart_mass, min_mass);
@@ -1544,10 +2048,13 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
* @param ind The array of indices to fill.
* @param cell_counts The cell counters to update.
* @param count_inhibited_parts (return) The number of #part to remove.
+ * @param count_extra_parts (return) The number of #part for on-the-fly
+ * creation.
* @param verbose Are we talkative ?
*/
void space_parts_get_cell_index(struct space *s, int *ind, int *cell_counts,
- int *count_inhibited_parts, int verbose) {
+ size_t *count_inhibited_parts,
+ size_t *count_extra_parts, int verbose) {
const ticks tic = getticks();
@@ -1563,11 +2070,15 @@ void space_parts_get_cell_index(struct space *s, int *ind, int *cell_counts,
data.count_inhibited_part = 0;
data.count_inhibited_gpart = 0;
data.count_inhibited_spart = 0;
+ data.count_extra_part = 0;
+ data.count_extra_gpart = 0;
+ data.count_extra_spart = 0;
threadpool_map(&s->e->threadpool, space_parts_get_cell_index_mapper, s->parts,
s->nr_parts, sizeof(struct part), 0, &data);
*count_inhibited_parts = data.count_inhibited_part;
+ *count_extra_parts = data.count_extra_part;
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -1583,10 +2094,13 @@ void space_parts_get_cell_index(struct space *s, int *ind, int *cell_counts,
* @param gind The array of indices to fill.
* @param cell_counts The cell counters to update.
* @param count_inhibited_gparts (return) The number of #gpart to remove.
+ * @param count_extra_gparts (return) The number of #gpart for on-the-fly
+ * creation.
* @param verbose Are we talkative ?
*/
void space_gparts_get_cell_index(struct space *s, int *gind, int *cell_counts,
- int *count_inhibited_gparts, int verbose) {
+ size_t *count_inhibited_gparts,
+ size_t *count_extra_gparts, int verbose) {
const ticks tic = getticks();
@@ -1602,11 +2116,15 @@ void space_gparts_get_cell_index(struct space *s, int *gind, int *cell_counts,
data.count_inhibited_part = 0;
data.count_inhibited_gpart = 0;
data.count_inhibited_spart = 0;
+ data.count_extra_part = 0;
+ data.count_extra_gpart = 0;
+ data.count_extra_spart = 0;
threadpool_map(&s->e->threadpool, space_gparts_get_cell_index_mapper,
s->gparts, s->nr_gparts, sizeof(struct gpart), 0, &data);
*count_inhibited_gparts = data.count_inhibited_gpart;
+ *count_extra_gparts = data.count_extra_gpart;
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -1622,10 +2140,13 @@ void space_gparts_get_cell_index(struct space *s, int *gind, int *cell_counts,
* @param sind The array of indices to fill.
* @param cell_counts The cell counters to update.
* @param count_inhibited_sparts (return) The number of #spart to remove.
+ * @param count_extra_sparts (return) The number of #spart for on-the-fly
+ * creation.
* @param verbose Are we talkative ?
*/
void space_sparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
- int *count_inhibited_sparts, int verbose) {
+ size_t *count_inhibited_sparts,
+ size_t *count_extra_sparts, int verbose) {
const ticks tic = getticks();
@@ -1641,11 +2162,15 @@ void space_sparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
data.count_inhibited_part = 0;
data.count_inhibited_gpart = 0;
data.count_inhibited_spart = 0;
+ data.count_extra_part = 0;
+ data.count_extra_gpart = 0;
+ data.count_extra_spart = 0;
threadpool_map(&s->e->threadpool, space_sparts_get_cell_index_mapper,
s->sparts, s->nr_sparts, sizeof(struct spart), 0, &data);
*count_inhibited_sparts = data.count_inhibited_spart;
+ *count_extra_sparts = data.count_extra_spart;
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -2018,7 +2543,7 @@ void space_split_recursive(struct space *s, struct cell *c,
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
- const int with_gravity = s->gravity;
+ const int with_self_gravity = s->with_self_gravity;
const int depth = c->depth;
int maxdepth = 0;
float h_max = 0.0f;
@@ -2046,6 +2571,8 @@ void space_split_recursive(struct space *s, struct cell *c,
#ifdef SWIFT_DEBUG_CHECKS
if (parts[k].time_bin == time_bin_inhibited)
error("Inhibited particle present in space_split()");
+ if (parts[k].time_bin == time_bin_not_created)
+ error("Extra particle present in space_split()");
#endif
buff[k].x[0] = parts[k].x[0];
buff[k].x[1] = parts[k].x[1];
@@ -2060,6 +2587,8 @@ void space_split_recursive(struct space *s, struct cell *c,
#ifdef SWIFT_DEBUG_CHECKS
if (gparts[k].time_bin == time_bin_inhibited)
error("Inhibited particle present in space_split()");
+ if (gparts[k].time_bin == time_bin_not_created)
+ error("Extra particle present in space_split()");
#endif
gbuff[k].x[0] = gparts[k].x[0];
gbuff[k].x[1] = gparts[k].x[1];
@@ -2074,6 +2603,8 @@ void space_split_recursive(struct space *s, struct cell *c,
#ifdef SWIFT_DEBUG_CHECKS
if (sparts[k].time_bin == time_bin_inhibited)
error("Inhibited particle present in space_split()");
+ if (sparts[k].time_bin == time_bin_not_created)
+ error("Extra particle present in space_split()");
#endif
sbuff[k].x[0] = sparts[k].x[0];
sbuff[k].x[1] = sparts[k].x[1];
@@ -2094,8 +2625,8 @@ void space_split_recursive(struct space *s, struct cell *c,
}
/* Split or let it be? */
- if ((with_gravity && gcount > space_splitsize) ||
- (!with_gravity &&
+ if ((with_self_gravity && gcount > space_splitsize) ||
+ (!with_self_gravity &&
(count > space_splitsize || scount > space_splitsize))) {
/* No longer just a leaf. */
@@ -2108,6 +2639,9 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->hydro.count = 0;
cp->grav.count = 0;
cp->stars.count = 0;
+ cp->hydro.count_total = 0;
+ cp->grav.count_total = 0;
+ cp->stars.count_total = 0;
cp->hydro.ti_old_part = c->hydro.ti_old_part;
cp->grav.ti_old_part = c->grav.ti_old_part;
cp->grav.ti_old_multipole = c->grav.ti_old_multipole;
@@ -2193,7 +2727,7 @@ void space_split_recursive(struct space *s, struct cell *c,
}
/* Deal with the multipole */
- if (s->gravity) {
+ if (s->with_self_gravity) {
/* Reset everything */
gravity_reset(c->grav.multipole);
@@ -2317,6 +2851,8 @@ void space_split_recursive(struct space *s, struct cell *c,
/* parts: Get dt_min/dt_max and h_max. */
for (int k = 0; k < count; k++) {
#ifdef SWIFT_DEBUG_CHECKS
+ if (parts[k].time_bin == time_bin_not_created)
+ error("Extra particle present in space_split()");
if (parts[k].time_bin == time_bin_inhibited)
error("Inhibited particle present in space_split()");
#endif
@@ -2335,6 +2871,8 @@ void space_split_recursive(struct space *s, struct cell *c,
/* gparts: Get dt_min/dt_max. */
for (int k = 0; k < gcount; k++) {
#ifdef SWIFT_DEBUG_CHECKS
+ if (gparts[k].time_bin == time_bin_not_created)
+ error("Extra g-particle present in space_split()");
if (gparts[k].time_bin == time_bin_inhibited)
error("Inhibited g-particle present in space_split()");
#endif
@@ -2345,6 +2883,8 @@ void space_split_recursive(struct space *s, struct cell *c,
/* sparts: Get dt_min/dt_max */
for (int k = 0; k < scount; k++) {
#ifdef SWIFT_DEBUG_CHECKS
+ if (sparts[k].time_bin == time_bin_not_created)
+ error("Extra s-particle present in space_split()");
if (sparts[k].time_bin == time_bin_inhibited)
error("Inhibited s-particle present in space_split()");
#endif
@@ -2372,7 +2912,7 @@ void space_split_recursive(struct space *s, struct cell *c,
ti_stars_end_min = get_integer_time_end(ti_current, stars_time_bin_min);
/* Construct the multipole and the centre of mass*/
- if (s->gravity) {
+ if (s->with_self_gravity) {
if (gcount > 0) {
gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count);
@@ -2480,7 +3020,7 @@ void space_recycle(struct space *s, struct cell *c) {
lock_lock(&s->lock);
/* Hook the multipole back in the buffer */
- if (s->gravity) {
+ if (s->with_self_gravity) {
c->grav.multipole->next = s->multipoles_sub;
s->multipoles_sub = c->grav.multipole;
}
@@ -2538,7 +3078,7 @@ void space_recycle_list(struct space *s, struct cell *cell_list_begin,
s->tot_cells -= count;
/* Hook the multipoles into the buffer. */
- if (s->gravity) {
+ if (s->with_self_gravity) {
multipole_list_end->next = s->multipoles_sub;
s->multipoles_sub = multipole_list_begin;
}
@@ -2582,7 +3122,7 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
}
/* Is the multipole buffer empty? */
- if (s->gravity && s->multipoles_sub == NULL) {
+ if (s->with_self_gravity && s->multipoles_sub == NULL) {
if (posix_memalign(
(void **)&s->multipoles_sub, multipole_align,
space_cellallocchunk * sizeof(struct gravity_tensors)) != 0)
@@ -2600,7 +3140,7 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
s->tot_cells += 1;
/* Hook the multipole */
- if (s->gravity) {
+ if (s->with_self_gravity) {
cells[j]->grav.multipole = s->multipoles_sub;
s->multipoles_sub = cells[j]->grav.multipole->next;
}
@@ -3045,8 +3585,11 @@ void space_convert_quantities_mapper(void *restrict map_data, int count,
const ptrdiff_t index = parts - s->parts;
struct xpart *restrict xparts = s->xparts + index;
+ /* Loop over all the particles ignoring the extra buffer ones for on-the-fly
+ * creation */
for (int k = 0; k < count; k++)
- hydro_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props);
+ if (parts[k].time_bin <= num_time_bins)
+ hydro_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props);
}
/**
@@ -3087,6 +3630,7 @@ void space_convert_quantities(struct space *s, int verbose) {
* @param generate_gas_in_ics Are we generating gas particles from the gparts?
* @param hydro flag whether we are doing hydro or not?
* @param self_gravity flag whether we are doing gravity or not?
+ * @param star_formation flag whether we are doing star formation or not?
* @param verbose Print messages to stdout or not.
* @param dry_run If 1, just initialise stuff, don't do anything with the parts.
*
@@ -3100,7 +3644,8 @@ void space_init(struct space *s, struct swift_params *params,
struct part *parts, struct gpart *gparts, struct spart *sparts,
size_t Npart, size_t Ngpart, size_t Nspart, int periodic,
int replicate, int generate_gas_in_ics, int hydro,
- int self_gravity, int verbose, int dry_run) {
+ int self_gravity, int star_formation, int verbose,
+ int dry_run) {
/* Clean-up everything */
bzero(s, sizeof(struct space));
@@ -3110,16 +3655,23 @@ void space_init(struct space *s, struct swift_params *params,
s->dim[1] = dim[1];
s->dim[2] = dim[2];
s->periodic = periodic;
- s->gravity = self_gravity;
- s->hydro = hydro;
+ s->with_self_gravity = self_gravity;
+ s->with_hydro = hydro;
+ s->with_star_formation = star_formation;
s->nr_parts = Npart;
- s->size_parts = Npart;
- s->parts = parts;
s->nr_gparts = Ngpart;
- s->size_gparts = Ngpart;
- s->gparts = gparts;
s->nr_sparts = Nspart;
+ s->size_parts = Npart;
+ s->size_gparts = Ngpart;
s->size_sparts = Nspart;
+ s->nr_inhibited_parts = 0;
+ s->nr_inhibited_gparts = 0;
+ s->nr_inhibited_sparts = 0;
+ s->nr_extra_parts = 0;
+ s->nr_extra_gparts = 0;
+ s->nr_extra_sparts = 0;
+ s->parts = parts;
+ s->gparts = gparts;
s->sparts = sparts;
s->min_part_mass = FLT_MAX;
s->min_gpart_mass = FLT_MAX;
@@ -3202,6 +3754,12 @@ void space_init(struct space *s, struct swift_params *params,
space_subdepth_diff_grav =
parser_get_opt_param_int(params, "Scheduler:cell_subdepth_diff_grav",
space_subdepth_diff_grav_default);
+ space_extra_parts = parser_get_opt_param_int(
+ params, "Scheduler:cell_extra_parts", space_extra_parts_default);
+ space_extra_sparts = parser_get_opt_param_int(
+ params, "Scheduler:cell_extra_sparts", space_extra_sparts_default);
+ space_extra_gparts = parser_get_opt_param_int(
+ params, "Scheduler:cell_extra_gparts", space_extra_gparts_default);
if (verbose) {
message("max_size set to %d split_size set to %d", space_maxsize,
@@ -3308,6 +3866,9 @@ void space_init(struct space *s, struct swift_params *params,
last_cell_id = 1;
#endif
+ /* Do we want any spare particles for on the fly cration? */
+ if (!star_formation) space_extra_sparts = 0;
+
/* Build the cells recursively. */
if (!dry_run) space_regrid(s, verbose);
}
@@ -3453,7 +4014,7 @@ void space_generate_gas(struct space *s, const struct cosmology *cosmo,
int periodic, const double dim[3], int verbose) {
/* Check that this is a sensible ting to do */
- if (!s->hydro)
+ if (!s->with_hydro)
error(
"Cannot generate gas from ICs if we are running without "
"hydrodynamics. Need to run with -s and the corresponding "
diff --git a/src/space.h b/src/space.h
index e6d774200be1a31d622419dceafb16b3826ce177..a1280945d2aa232cbb5e5b519266bc7058e5dc57 100644
--- a/src/space.h
+++ b/src/space.h
@@ -44,6 +44,10 @@ struct cosmology;
#define space_cellallocchunk 1000
#define space_splitsize_default 400
#define space_maxsize_default 8000000
+#define space_extra_parts_default 0
+#define space_extra_gparts_default 0
+#define space_extra_sparts_default 100
+#define space_expected_max_nr_strays_default 100
#define space_subsize_pair_hydro_default 256000000
#define space_subsize_self_hydro_default 32000
#define space_subsize_pair_grav_default 256000000
@@ -68,6 +72,9 @@ extern int space_subsize_self_grav;
extern int space_subsize_pair_stars;
extern int space_subsize_self_stars;
extern int space_subdepth_diff_grav;
+extern int space_extra_parts;
+extern int space_extra_gparts;
+extern int space_extra_sparts;
/**
* @brief The space in which the cells and particles reside.
@@ -84,10 +91,13 @@ struct space {
struct hydro_space hs;
/*! Are we doing hydrodynamics? */
- int hydro;
+ int with_hydro;
/*! Are we doing gravity? */
- int gravity;
+ int with_self_gravity;
+
+ /*! Are we doing star formation? */
+ int with_star_formation;
/*! Width of the top-level cells. */
double width[3];
@@ -149,14 +159,23 @@ struct space {
/*! The indices of the top-level cells that have >0 particles (of any kind) */
int *local_cells_with_particles_top;
- /*! The total number of parts in the space. */
- size_t nr_parts, size_parts;
+ /*! The total number of #part in the space. */
+ size_t nr_parts;
+
+ /*! The total number of #gpart in the space. */
+ size_t nr_gparts;
+
+ /*! The total number of #spart in the space. */
+ size_t nr_sparts;
+
+ /*! The total number of #part we allocated memory for */
+ size_t size_parts;
- /*! The total number of g-parts in the space. */
- size_t nr_gparts, size_gparts;
+ /*! The total number of #gpart we allocated memory for */
+ size_t size_gparts;
- /*! The total number of g-parts in the space. */
- size_t nr_sparts, size_sparts;
+ /*! The total number of #spart we allocated memory for */
+ size_t size_sparts;
/*! Number of inhibted gas particles in the space */
size_t nr_inhibited_parts;
@@ -167,6 +186,15 @@ struct space {
/*! Number of inhibted star particles in the space */
size_t nr_inhibited_sparts;
+ /*! Number of extra #part we allocated (for on-the-fly creation) */
+ size_t nr_extra_parts;
+
+ /*! Number of extra #gpart we allocated (for on-the-fly creation) */
+ size_t nr_extra_gparts;
+
+ /*! Number of extra #spart we allocated (for on-the-fly creation) */
+ size_t nr_extra_sparts;
+
/*! The particle data (cells have pointers to this). */
struct part *parts;
@@ -241,7 +269,7 @@ void space_init(struct space *s, struct swift_params *params,
struct part *parts, struct gpart *gparts, struct spart *sparts,
size_t Npart, size_t Ngpart, size_t Nspart, int periodic,
int replicate, int generate_gas_in_ics, int hydro, int gravity,
- int verbose, int dry_run);
+ int star_formation, int verbose, int dry_run);
void space_sanitize(struct space *s);
void space_map_cells_pre(struct space *s, int full,
void (*fun)(struct cell *c, void *data), void *data);
@@ -260,14 +288,18 @@ void space_recycle_list(struct space *s, struct cell *cell_list_begin,
struct gravity_tensors *multipole_list_begin,
struct gravity_tensors *multipole_list_end);
void space_split(struct space *s, int verbose);
+void space_reorder_extras(struct space *s, int verbose);
void space_split_mapper(void *map_data, int num_elements, void *extra_data);
void space_list_useful_top_level_cells(struct space *s);
void space_parts_get_cell_index(struct space *s, int *ind, int *cell_counts,
- int *count_inibibited_parts, int verbose);
+ size_t *count_inhibited_parts,
+ size_t *count_extra_parts, int verbose);
void space_gparts_get_cell_index(struct space *s, int *gind, int *cell_counts,
- int *count_inibibited_gparts, int verbose);
+ size_t *count_inhibited_gparts,
+ size_t *count_extra_gparts, int verbose);
void space_sparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
- int *count_inibibited_sparts, int verbose);
+ size_t *count_inhibited_sparts,
+ size_t *count_extra_sparts, int verbose);
void space_synchronize_particle_positions(struct space *s);
void space_do_parts_sort(void);
void space_do_gparts_sort(void);
diff --git a/src/swift.h b/src/swift.h
index 153c4ae0d4440d083f1b0c9850e1f2649c0df6fb..7279fc58b3b4ac0227c11b8a835c829df2ceb49c 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -65,7 +65,6 @@
#include "scheduler.h"
#include "serial_io.h"
#include "single_io.h"
-#include "sourceterms.h"
#include "space.h"
#include "stars.h"
#include "stars_io.h"
diff --git a/src/task.c b/src/task.c
index 3918dad3b713c6c226e5dacf3e38756910c1dd27..f68d10f7cfbcfd512151e5bd4c21089128a6b804 100644
--- a/src/task.c
+++ b/src/task.c
@@ -75,7 +75,6 @@ const char *taskID_names[task_type_count] = {"none",
"grav_mesh",
"cooling",
"star_formation",
- "sourceterms",
"logger",
"stars_ghost_in",
"stars_ghost",
@@ -141,7 +140,6 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_ghost:
case task_type_extra_ghost:
case task_type_cooling:
- case task_type_sourceterms:
return task_action_part;
break;
@@ -391,6 +389,12 @@ void task_unlock(struct task *t) {
cell_munlocktree(cj);
break;
+ case task_type_star_formation:
+ cell_unlocktree(ci);
+ cell_sunlocktree(ci);
+ cell_gunlocktree(ci);
+ break;
+
default:
break;
}
@@ -545,6 +549,21 @@ int task_lock(struct task *t) {
cell_munlocktree(ci);
return 0;
}
+ break;
+
+ case task_type_star_formation:
+ /* Lock the gas, gravity and star particles */
+ if (ci->hydro.hold || ci->stars.hold || ci->grav.phold) return 0;
+ if (cell_locktree(ci) != 0) return 0;
+ if (cell_slocktree(ci) != 0) {
+ cell_unlocktree(ci);
+ return 0;
+ }
+ if (cell_glocktree(ci) != 0) {
+ cell_unlocktree(ci);
+ cell_sunlocktree(ci);
+ return 0;
+ }
default:
break;
diff --git a/src/task.h b/src/task.h
index 994b2b14c05965b71e877feac5cb9827a1d1b4bb..5bce55d6a28fae7ad99f966267b14e926b7fd924 100644
--- a/src/task.h
+++ b/src/task.h
@@ -67,7 +67,6 @@ enum task_types {
task_type_grav_mesh,
task_type_cooling,
task_type_star_formation,
- task_type_sourceterms,
task_type_logger,
task_type_stars_ghost_in,
task_type_stars_ghost,
diff --git a/src/timeline.h b/src/timeline.h
index 38727def50b5a81c073ab23375f0c548ca096b66..09161b3ec68e9f7c19bf812103384257473f5719 100644
--- a/src/timeline.h
+++ b/src/timeline.h
@@ -40,6 +40,9 @@ typedef char timebin_t;
/*! Fictious time-bin to hold inhibited particles */
#define time_bin_inhibited (num_time_bins + 2)
+/*! Fictious time-bin to hold particles not yet created */
+#define time_bin_not_created (num_time_bins + 3)
+
/*! Fictitious time-bin for particles not awaken */
#define time_bin_not_awake (0)
diff --git a/src/timers.c b/src/timers.c
index d7523edaa42570f8e5f3c01516075267988dfd9c..9ede0320e49c70b2488cd5ceb3e4b6965659aa74 100644
--- a/src/timers.c
+++ b/src/timers.c
@@ -61,7 +61,6 @@ const char* timers_names[timer_count] = {
"dograv_mesh",
"dograv_top_level",
"dograv_long_range",
- "dosource",
"dosub_self_density",
"dosub_self_gradient",
"dosub_self_force",
diff --git a/src/timers.h b/src/timers.h
index 412f36d7de547040e29efdede9cc6826358929bc..3a2a939339e6d08b43836d4f5ca213af0822c2b2 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -62,7 +62,6 @@ enum {
timer_dograv_mesh,
timer_dograv_top_level,
timer_dograv_long_range,
- timer_dosource,
timer_dosub_self_density,
timer_dosub_self_gradient,
timer_dosub_self_force,
diff --git a/tools/task_plots/analyse_tasks.py b/tools/task_plots/analyse_tasks.py
index ca41970c683a1680e9d1054c9d70d6370992a05e..f79a0090b04e013c9e10c53e6a4796e46579e1b9 100755
--- a/tools/task_plots/analyse_tasks.py
+++ b/tools/task_plots/analyse_tasks.py
@@ -91,7 +91,6 @@ TASKTYPES = [
"grav_mesh",
"cooling",
"star_formation",
- "sourceterms",
"logger",
"stars_ghost_in",
"stars_ghost",
diff --git a/tools/task_plots/plot_tasks.py b/tools/task_plots/plot_tasks.py
index 1fe7bcbd11f30ff17051bc9a7ae789439df8b9e9..1ff722a6079883a043109e26dea8d6b3fb405000 100755
--- a/tools/task_plots/plot_tasks.py
+++ b/tools/task_plots/plot_tasks.py
@@ -176,7 +176,6 @@ TASKTYPES = [
"grav_mesh",
"cooling",
"star_formation",
- "sourceterms",
"logger",
"stars_ghost_in",
"stars_ghost",