From e5b7288a777686a5294cfdb4188da90e69977fd0 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Sun, 11 Nov 2018 13:41:57 +0100 Subject: [PATCH 01/52] Re-ordered the content of the cell structure to save memory. --- src/cell.h | 210 ++++++++++++++++++++++++++--------------------------- 1 file changed, 105 insertions(+), 105 deletions(-) diff --git a/src/cell.h b/src/cell.h index 89f7c954c2..9a9f0a1987 100644 --- a/src/cell.h +++ b/src/cell.h @@ -211,12 +211,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; @@ -239,18 +239,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; @@ -261,20 +288,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; @@ -285,12 +306,27 @@ 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; + + /*! 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; @@ -300,42 +336,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. */ @@ -358,6 +358,33 @@ struct cell { * tasks */ struct cell *super; + /*! The drift task for gparts */ + struct task *drift; + + /*! 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; @@ -374,15 +401,15 @@ 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; + /*! Number of #gpart updated in this cell. */ int updated; @@ -395,42 +422,15 @@ 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; - - /*! 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; - - /*! Number of M-M tasks that are associated with this cell. */ - short int nr_mm_tasks; - } grav; /*! Stars variables */ @@ -439,12 +439,27 @@ struct cell { /*! Pointer to the #spart data. */ struct spart *parts; - /*! Nr of #spart in this cell. */ - int count; + /*! 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; + + /*! Linked list of the tasks computing this cell's star density. */ + struct link *density; /*! Max smoothing length in this cell. */ double h_max; + /*! Spin lock for various uses (#spart case). */ + swift_lock_type lock; + + /*! Nr of #spart in this cell. */ + int count; + /*! Values of h_max before the drifts, used for sub-cell tasks. */ float h_max_old; @@ -454,18 +469,6 @@ struct cell { /*! Values of dx_max before the drifts, used for sub-cell tasks. */ float dx_max_part_old; - /*! 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; - - /*! Linked list of the tasks computing this cell's star density. */ - struct link *density; - /*! Number of #spart updated in this cell. */ int updated; @@ -475,9 +478,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; - } stars; #ifdef WITH_MPI -- GitLab From 9e77d206892f2e1dbfb41bd9b36e225b0a800f4b Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Sun, 11 Nov 2018 15:40:46 +0100 Subject: [PATCH 02/52] Clarify some magical constants involved in the allocation and resizing of particle arrays. --- src/engine.c | 12 ++-- src/space.c | 184 ++++++++++++++++++++++++++++----------------------- src/space.h | 2 + 3 files changed, 110 insertions(+), 88 deletions(-) diff --git a/src/engine.c b/src/engine.c index f728862487..459ebea37c 100644 --- a/src/engine.c +++ b/src/engine.c @@ -3871,8 +3871,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, @@ -3896,8 +3896,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) @@ -3914,8 +3914,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) diff --git a/src/space.c b/src/space.c index bda2f48825..411d24a926 100644 --- a/src/space.c +++ b/src/space.c @@ -70,6 +70,12 @@ int space_subsize_pair_stars = space_subsize_pair_stars_default; int space_subsize_self_stars = space_subsize_self_stars_default; int space_subdepth_grav = space_subdepth_grav_default; int space_maxsize = space_maxsize_default; + +/*! Number of extra #spart we allocate memory for per top-level cell */ +int space_extra_sparts = space_extra_sparts_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 @@ -603,48 +609,61 @@ 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); + struct cell *cells_top = s->cells_top; + const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0; + + /* 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; + + /* 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; + + /* Number of inhibited particles found on the node */ 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; + + /* 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; /* 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 *h_index = (int *)malloc(sizeof(int) * h_index_size); + if (h_index == 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); + if (size_parts > 0) + space_parts_get_cell_index(s, h_index, cell_part_counts, + &count_inhibited_parts, verbose); /* 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 *g_index = (int *)malloc(sizeof(int) * g_index_size); + if (g_index == 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, + if (size_gparts > 0) + space_gparts_get_cell_index(s, g_index, 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 *s_index = (int *)malloc(sizeof(int) * s_index_size); + if (s_index == 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, + if (size_sparts > 0) + space_sparts_get_cell_index(s, s_index, cell_spart_counts, &count_inhibited_sparts, verbose); #ifdef SWIFT_DEBUG_CHECKS @@ -663,7 +682,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { 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; @@ -682,7 +701,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".*/ @@ -695,15 +714,15 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { /* Check that all parts are in the correct places. */ int 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!"); @@ -714,7 +733,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; @@ -731,7 +750,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".*/ @@ -744,15 +763,15 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { /* Check that all sparts are in the correct place. */ int 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!"); @@ -763,7 +782,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; @@ -786,7 +805,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++; @@ -798,15 +817,15 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { /* Check that all gparts are in the correct place. */ int 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!"); @@ -815,16 +834,17 @@ 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; @@ -851,23 +871,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]}; @@ -876,13 +896,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; @@ -890,13 +910,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; @@ -911,8 +931,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. */ @@ -930,7 +950,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] || @@ -942,7 +962,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. */ @@ -960,7 +980,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] || @@ -972,52 +992,52 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { /* Extract the cell counts from the sorted indices. */ 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; last_index = k + 1; } } /* Extract the cell counts from the sorted indices. */ 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; 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; @@ -1036,8 +1056,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. */ @@ -1055,7 +1075,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] || @@ -1067,16 +1087,16 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { /* Extract the cell counts from the sorted indices. */ 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; 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 diff --git a/src/space.h b/src/space.h index d0ec372d6f..93a8dc83b7 100644 --- a/src/space.h +++ b/src/space.h @@ -44,6 +44,8 @@ struct cosmology; #define space_cellallocchunk 1000 #define space_splitsize_default 400 #define space_maxsize_default 8000000 +#define space_extra_sparts_default 200 +#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 -- GitLab From e8faa2fde812707a1f070dc874ed025acbdaf47c Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Sun, 11 Nov 2018 16:11:33 +0100 Subject: [PATCH 03/52] Get the correct number of tasks when running with star formation. --- src/engine.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/engine.c b/src/engine.c index 459ebea37c..a69c13a161 100644 --- a/src/engine.c +++ b/src/engine.c @@ -1909,6 +1909,9 @@ int engine_estimate_nr_tasks(struct engine *e) { if (e->policy & engine_policy_cooling) { n1 += 2; } + if (e->policy & engine_policy_star_formation) { + n1 += 1; + } if (e->policy & engine_policy_sourceterms) { n1 += 2; } -- GitLab From e6f99a0662ebf7f11f26a961e9cf467198b4988e Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Sun, 11 Nov 2018 16:21:36 +0100 Subject: [PATCH 04/52] Added debugging check to prevent inhibited/not-created particles in interactions. --- src/hydro/Gadget2/hydro_iact.h | 28 ++++++++++++++++++++++++++++ src/hydro/Minimal/hydro_iact.h | 16 ++++++++-------- 2 files changed, 36 insertions(+), 8 deletions(-) diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h index 746fd47785..a3c5e21dbd 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 e060cb3562..b29f44588c 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 -- GitLab From 860c721de6bc93848a25b39f7d92a74e22f11df8 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Sun, 11 Nov 2018 16:22:24 +0100 Subject: [PATCH 05/52] Added a new fake time-bin for particles not yet created/added to the calculation. --- src/timeline.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/timeline.h b/src/timeline.h index 38727def50..09161b3ec6 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) -- GitLab From a15b59a0fb242981ae87b98d4f9aa486cf83437d Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 12 Nov 2018 10:05:48 +0100 Subject: [PATCH 06/52] Infrastructure in place for adding future particles at rebuild time. --- examples/main.c | 2 +- src/cell.c | 2 + src/cell.h | 3 ++ src/space.c | 97 ++++++++++++++++++++++++++++++++++++------------- src/space.h | 4 +- 5 files changed, 81 insertions(+), 27 deletions(-) diff --git a/examples/main.c b/examples/main.c index 5ae01a9280..049d49fa32 100644 --- a/examples/main.c +++ b/examples/main.c @@ -938,7 +938,7 @@ int main(int argc, char *argv[]) { engine_policies |= engine_policy_structure_finding; // MATTHIEU: Temporary star formation law - // engine_policies |= engine_policy_star_formation; + engine_policies |= engine_policy_star_formation; /* Initialize the engine with the space and policies. */ if (myrank == 0) clocks_gettime(&tic); diff --git a/src/cell.c b/src/cell.c index 37b7be42f3..1c4d1f9e32 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3633,6 +3633,8 @@ void cell_check_timesteps(struct cell *c) { #endif } +void cell_add_spart(const struct engine *e, struct cell *c) {} + /** * @brief "Remove" a gas particle from the calculation. * diff --git a/src/cell.h b/src/cell.h index 9a9f0a1987..5d542e5f2d 100644 --- a/src/cell.h +++ b/src/cell.h @@ -460,6 +460,9 @@ struct cell { /*! Nr of #spart in this cell. */ int count; + /*! 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; diff --git a/src/space.c b/src/space.c index 411d24a926..fd0bd9e752 100644 --- a/src/space.c +++ b/src/space.c @@ -71,9 +71,15 @@ int space_subsize_self_stars = space_subsize_self_stars_default; int space_subdepth_grav = space_subdepth_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 @@ -285,6 +291,8 @@ void space_regrid(struct space *s, int verbose) { const ticks tic = getticks(); const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0; + message("REGRID!!"); + /* Run through the cells and get the current h_max. */ // tic = getticks(); float h_max = s->cell_min / kernel_gamma / space_stretch; @@ -589,6 +597,41 @@ void space_regrid(struct space *s, int verbose) { clocks_getunit()); } +void space_allocate_extras(struct space *s, int verbose) { + + 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; + + /* 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++; + + const size_t num_extra_parts = local_cells * space_extra_parts; + const size_t num_extra_gparts = local_cells * space_extra_parts; + const size_t num_extra_sparts = local_cells * space_extra_parts; + + if (verbose) + message("Requesting space for future %zd/%zd/%zd part/gpart/sparts.", + num_extra_parts, num_extra_gparts, num_extra_sparts); + + /* Do we need to reallocate? */ + if (nr_parts + num_extra_parts > size_parts) { + } + if (nr_gparts + num_extra_gparts > size_gparts) { + } + if (nr_sparts + num_extra_sparts > size_sparts) { + } +} + /** * @brief Re-build the cells as well as the tasks. * @@ -609,8 +652,12 @@ 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 */ + 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; @@ -632,36 +679,37 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { 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; - /* 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. */ + /* 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); - if (h_index == 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."); + 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] = space_extra_parts; + cell_gpart_counts[i] = space_extra_gparts; + cell_spart_counts[i] = space_extra_sparts; + } + + /* Run through the particles and get their cell index. */ if (size_parts > 0) space_parts_get_cell_index(s, h_index, cell_part_counts, &count_inhibited_parts, verbose); - - /* Run through the gravity particles and get their cell index. */ - int *g_index = (int *)malloc(sizeof(int) * g_index_size); - if (g_index == 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 (size_gparts > 0) space_gparts_get_cell_index(s, g_index, cell_gpart_counts, &count_inhibited_gparts, verbose); - - /* Run through the star particles and get their cell index. */ - int *s_index = (int *)malloc(sizeof(int) * s_index_size); - if (s_index == 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 (size_sparts > 0) space_sparts_get_cell_index(s, s_index, cell_spart_counts, &count_inhibited_sparts, verbose); @@ -675,8 +723,6 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { error("We just repartitioned but still found inhibited gparts."); #endif - const int local_nodeID = s->e->nodeID; - /* 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 */) { @@ -850,6 +896,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { 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) diff --git a/src/space.h b/src/space.h index 93a8dc83b7..2444a2ed88 100644 --- a/src/space.h +++ b/src/space.h @@ -44,7 +44,9 @@ struct cosmology; #define space_cellallocchunk 1000 #define space_splitsize_default 400 #define space_maxsize_default 8000000 -#define space_extra_sparts_default 200 +#define space_extra_parts_default 10 +#define space_extra_gparts_default 10 +#define space_extra_sparts_default 10 #define space_expected_max_nr_strays_default 100 #define space_subsize_pair_hydro_default 256000000 #define space_subsize_self_hydro_default 32000 -- GitLab From 1b829c0be4ab59879df13d3f9ee98bdf1a96ab6a Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 12 Nov 2018 10:43:55 +0100 Subject: [PATCH 07/52] Clarified some space variables. Added space-wide counters of extra particles. --- src/space.c | 51 +++++++++++++++++++++++++++++---------------------- src/space.h | 34 ++++++++++++++++++++++++++-------- 2 files changed, 55 insertions(+), 30 deletions(-) diff --git a/src/space.c b/src/space.c index fd0bd9e752..db50196b00 100644 --- a/src/space.c +++ b/src/space.c @@ -148,13 +148,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; @@ -238,7 +238,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]); @@ -448,7 +449,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."); @@ -520,7 +521,7 @@ void space_regrid(struct space *s, int verbose) { #ifdef WITH_MPI c->mpi.tag = -1; #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++; @@ -1214,7 +1215,7 @@ void space_rebuild(struct space *s, int repartitioned, int 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 @@ -2061,7 +2062,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; @@ -2136,8 +2137,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. */ @@ -2232,7 +2233,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); @@ -2407,7 +2408,7 @@ void space_split_recursive(struct space *s, struct cell *c, get_integer_time_begin(ti_current + 1, gravity_time_bin_max); /* 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); @@ -2514,7 +2515,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; } @@ -2572,7 +2573,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; } @@ -2616,7 +2617,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) @@ -2634,7 +2635,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; } @@ -3137,16 +3138,22 @@ 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->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; @@ -3479,7 +3486,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 2444a2ed88..5bbf548454 100644 --- a/src/space.h +++ b/src/space.h @@ -88,10 +88,10 @@ 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; /*! Width of the top-level cells. */ double width[3]; @@ -153,14 +153,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 g-parts in the space. */ - size_t nr_gparts, size_gparts; + /*! The total number of #gpart in the space. */ + size_t nr_gparts; - /*! The total number of g-parts in the space. */ - size_t nr_sparts, size_sparts; + /*! 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 #gpart we allocated memory for */ + size_t size_gparts; + + /*! 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; @@ -171,6 +180,15 @@ struct space { /*! Number of inhibted star particles in the space */ size_t nr_inhibited_sparts; + /*! Number of extra #part (for on-the-fly creation) we allocated */ + size_t nr_extra_parts; + + /*! Number of extra #gpart (for on-the-fly creation) we allocated */ + size_t nr_extra_gparts; + + /*! Number of extra #spart (for on-the-fly creation) we allocated */ + size_t nr_extra_sparts; + /*! The particle data (cells have pointers to this). */ struct part *parts; -- GitLab From 19a96c98529190c795189d5a62d746a466f7d582 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 12 Nov 2018 15:13:07 +0100 Subject: [PATCH 08/52] Count the extra particles in space_parts_get_cell_index to make sure everything went smoothly. --- src/space.c | 189 +++++++++++++++++++++++++++++++++++++++++----------- src/space.h | 13 ++-- 2 files changed, 159 insertions(+), 43 deletions(-) diff --git a/src/space.c b/src/space.c index db50196b00..1ec6cad4c5 100644 --- a/src/space.c +++ b/src/space.c @@ -116,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; }; /** @@ -602,11 +605,16 @@ void space_allocate_extras(struct space *s, int verbose) { const int local_nodeID = s->e->nodeID; - /* The current number of particles */ + /* 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; @@ -616,20 +624,73 @@ void space_allocate_extras(struct space *s, int verbose) { for (int i = 0; i < s->nr_cells; ++i) if (s->cells_top[i].nodeID == local_nodeID) local_cells++; - const size_t num_extra_parts = local_cells * space_extra_parts; - const size_t num_extra_gparts = local_cells * space_extra_parts; - const size_t num_extra_sparts = local_cells * space_extra_parts; + /* 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_parts; + const size_t expected_num_extra_sparts = local_cells * space_extra_parts; - if (verbose) + 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 future particles.", + s->nr_extra_parts, s->nr_extra_gparts, s->nr_extra_sparts); message("Requesting space for future %zd/%zd/%zd part/gpart/sparts.", - num_extra_parts, num_extra_gparts, num_extra_sparts); + expected_num_extra_parts, expected_num_extra_gparts, + expected_num_extra_sparts); + } + + /* Do we have enough space for the extras ? */ + if (expected_num_extra_parts > s->nr_extra_parts) { + + /* 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; + } - /* Do we need to reallocate? */ - if (nr_parts + num_extra_parts > size_parts) { + /* Do we need to get 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; + } + + /* Update the counters */ + s->nr_parts = nr_actual_parts + expected_num_extra_parts; + s->nr_extra_parts = expected_num_extra_parts; } - if (nr_gparts + num_extra_gparts > size_gparts) { + + if (nr_gparts + expected_num_extra_gparts > size_gparts) { } - if (nr_sparts + num_extra_sparts > size_sparts) { + if (nr_sparts + expected_num_extra_sparts > size_sparts) { } } @@ -671,9 +732,14 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { size_t size_sparts = s->size_sparts; /* Number of inhibited particles found on the node */ - int count_inhibited_parts = 0; - int count_inhibited_gparts = 0; - int count_inhibited_sparts = 0; + size_t count_inhibited_parts = 0; + size_t count_inhibited_gparts = 0; + size_t count_inhibited_sparts = 0; + + /* 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; @@ -699,29 +765,46 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { /* Initialise the counters, including buffer space for future particles */ for (int i = 0; i < s->nr_cells; ++i) { - cell_part_counts[i] = space_extra_parts; - cell_gpart_counts[i] = space_extra_gparts; - cell_spart_counts[i] = space_extra_sparts; + cell_part_counts[i] = 0; // space_extra_parts; + cell_gpart_counts[i] = 0; // space_extra_gparts; + cell_spart_counts[i] = 0; // space_extra_sparts; } /* Run through the particles and get their cell index. */ - if (size_parts > 0) + if (nr_parts > 0) space_parts_get_cell_index(s, h_index, cell_part_counts, - &count_inhibited_parts, verbose); - if (size_gparts > 0) + &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, verbose); - if (size_sparts > 0) + &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, verbose); + &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."); + + 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. */ @@ -759,7 +842,7 @@ 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 (h_index[k] == -1 || cells_top[h_index[k]].nodeID != local_nodeID) { error("Failed to move all non-local parts to send list"); @@ -808,7 +891,7 @@ 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 (s_index[k] == -1 || cells_top[s_index[k]].nodeID != local_nodeID) { error("Failed to move all non-local sparts to send list"); @@ -862,7 +945,7 @@ 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 (g_index[k] == -1 || cells_top[g_index[k]].nodeID != local_nodeID) { error("Failed to move all non-local gparts to send list"); @@ -1310,10 +1393,13 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts, if (cell_counts == NULL) error("Failed to allocate temporary cell count buffer."); + // message("s->nr_cells=%d", s->nr_cells); + /* 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++) { @@ -1345,12 +1431,16 @@ 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] = -2; + ++count_extra_part; } else { - /* List its top-level cell index */ + /* Normal case: list its top-level cell index */ ind[k] = index; cell_counts[index]++; } @@ -1372,8 +1462,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); @@ -1413,7 +1505,7 @@ 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; for (int k = 0; k < nr_gparts; k++) { @@ -1517,7 +1609,7 @@ 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; for (int k = 0; k < nr_sparts; k++) { @@ -1593,10 +1685,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(); @@ -1612,11 +1707,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), @@ -1632,10 +1731,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(); @@ -1651,11 +1753,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), @@ -1671,10 +1777,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(); @@ -1690,11 +1799,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), diff --git a/src/space.h b/src/space.h index 5bbf548454..a2ec8e7cfb 100644 --- a/src/space.h +++ b/src/space.h @@ -45,8 +45,8 @@ struct cosmology; #define space_splitsize_default 400 #define space_maxsize_default 8000000 #define space_extra_parts_default 10 -#define space_extra_gparts_default 10 -#define space_extra_sparts_default 10 +#define space_extra_gparts_default 0 +#define space_extra_sparts_default 0 #define space_expected_max_nr_strays_default 100 #define space_subsize_pair_hydro_default 256000000 #define space_subsize_self_hydro_default 32000 @@ -285,11 +285,14 @@ void space_split(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); -- GitLab From f6707353d3a00f42e28eaa505c7cfcb0b2f28eb4 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 12 Nov 2018 16:01:14 +0100 Subject: [PATCH 09/52] Put the extra particles in the centre of their top-level cells. --- src/space.c | 70 ++++++++++++++++++++++++++++++++++++++++++++--------- src/space.h | 6 ++--- 2 files changed, 62 insertions(+), 14 deletions(-) diff --git a/src/space.c b/src/space.c index 1ec6cad4c5..a27280d411 100644 --- a/src/space.c +++ b/src/space.c @@ -605,6 +605,12 @@ void space_allocate_extras(struct space *s, int verbose) { const int local_nodeID = s->e->nodeID; + /* 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; @@ -639,7 +645,7 @@ void space_allocate_extras(struct space *s, int verbose) { expected_num_extra_sparts); } - /* Do we have enough space for the extras ? */ + /* Do we have enough space for the extras (i.e. we haven't used up any) ? */ if (expected_num_extra_parts > s->nr_extra_parts) { /* Do we need to reallocate? */ @@ -674,8 +680,7 @@ void space_allocate_extras(struct space *s, int verbose) { s->size_parts = size_parts; } - /* Do we need to get some of the allocated spares into particles we can use - */ + /* 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)); @@ -683,6 +688,37 @@ void space_allocate_extras(struct space *s, int verbose) { s->parts[i].time_bin = time_bin_not_created; } + /* 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) { + 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; @@ -1393,8 +1429,6 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts, if (cell_counts == NULL) error("Failed to allocate temporary cell count buffer."); - // message("s->nr_cells=%d", s->nr_cells); - /* Init the local collectors */ float min_mass = FLT_MAX; float sum_vel_norm = 0.f; @@ -1437,7 +1471,8 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts, ++count_inhibited_part; } else if (p->time_bin == time_bin_not_created) { /* Is this a place-holder for on-the-fly creation? */ - ind[k] = -2; + ind[k] = index; + cell_counts[index]++; ++count_extra_part; } else { /* Normal case: list its top-level cell index */ @@ -1506,6 +1541,7 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts, float min_mass = FLT_MAX; float sum_vel_norm = 0.f; size_t count_inhibited_gpart = 0; + size_t count_extra_gpart = 0; for (int k = 0; k < nr_gparts; k++) { @@ -1536,10 +1572,15 @@ 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; @@ -1568,8 +1609,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); @@ -1610,6 +1654,7 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts, float min_mass = FLT_MAX; float sum_vel_norm = 0.f; size_t count_inhibited_spart = 0; + size_t count_extra_spart = 0; for (int k = 0; k < nr_sparts; k++) { @@ -1668,8 +1713,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); diff --git a/src/space.h b/src/space.h index a2ec8e7cfb..cdf83e4ada 100644 --- a/src/space.h +++ b/src/space.h @@ -180,13 +180,13 @@ struct space { /*! Number of inhibted star particles in the space */ size_t nr_inhibited_sparts; - /*! Number of extra #part (for on-the-fly creation) we allocated */ + /*! Number of extra #part we allocated (for on-the-fly creation) */ size_t nr_extra_parts; - /*! Number of extra #gpart (for on-the-fly creation) we allocated */ + /*! Number of extra #gpart we allocated (for on-the-fly creation) */ size_t nr_extra_gparts; - /*! Number of extra #spart (for on-the-fly creation) we allocated */ + /*! Number of extra #spart we allocated (for on-the-fly creation) */ size_t nr_extra_sparts; /*! The particle data (cells have pointers to this). */ -- GitLab From 44f618ee742f08d80b51e26cd5b425ae07fef644 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 12 Nov 2018 17:33:18 +0100 Subject: [PATCH 10/52] Sort the extra particles at the end of their top-level cells. --- src/cell.c | 29 +++++++++++++++++++ src/cell.h | 3 ++ src/space.c | 82 +++++++++++++++++++++++++++++++++++++++++++---------- src/space.h | 4 +++ 4 files changed, 103 insertions(+), 15 deletions(-) diff --git a/src/cell.c b/src/cell.c index 1c4d1f9e32..57795efcd3 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3794,6 +3794,35 @@ void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c, #endif } +void cell_reorder_extra_parts(struct cell *c) { + + struct part *parts = c->hydro.parts; + struct xpart *xparts = c->hydro.xparts; + const int count_real = c->hydro.count; + const int count_total = count_real + space_extra_parts; + + 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; + } + + if (first_not_extra >= count_total) + error("Looking for extra particles beyond this cell's range!"); + + memswap(&parts[i], &parts[first_not_extra], sizeof(struct part)); + memswap(&xparts[i], &xparts[first_not_extra], sizeof(struct xpart)); + if (parts[i].gpart) error("Need to handle gpart pointer!"); + } + } +} + /** * @brief Can we use the MM interactions fo a given pair of cells? * diff --git a/src/cell.h b/src/cell.h index 5d542e5f2d..009e6522db 100644 --- a/src/cell.h +++ b/src/cell.h @@ -663,6 +663,9 @@ 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); +void cell_reorder_extra_parts(struct cell *c); +void cell_reorder_extra_gparts(struct cell *c); +void cell_reorder_extra_sparts(struct cell *c); 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/space.c b/src/space.c index a27280d411..6f78ed8674 100644 --- a/src/space.c +++ b/src/space.c @@ -767,12 +767,12 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { size_t size_gparts = s->size_gparts; size_t size_sparts = s->size_sparts; - /* Number of inhibited particles found on the node */ + /* 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; - /* Number of extra particles found on the node */ + /* 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; @@ -1157,22 +1157,26 @@ 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; h_index[nr_parts] = s->nr_cells; // sentinel. for (size_t k = 0; k < nr_parts; k++) { if (h_index[k] < h_index[k + 1]) { - cells_top[h_index[k]].hydro.count = k - last_index + 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; s_index[nr_sparts] = s->nr_cells; // sentinel. for (size_t k = 0; k < nr_sparts; k++) { if (s_index[k] < s_index[k + 1]) { - cells_top[s_index[k]].stars.count = k - last_sindex + 1; + cells_top[s_index[k]].stars.count = + k - last_sindex + 1 - space_extra_sparts; last_sindex = k + 1; } } @@ -1252,12 +1256,14 @@ 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; g_index[nr_gparts] = s->nr_cells; for (size_t k = 0; k < nr_gparts; k++) { if (g_index[k] < g_index[k + 1]) { - cells_top[g_index[k]].grav.count = k - last_gindex + 1; + cells_top[g_index[k]].grav.count = + k - last_gindex + 1 - space_extra_gparts; last_gindex = k + 1; } } @@ -1302,10 +1308,10 @@ 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]; + finger = &finger[c->hydro.count + space_extra_parts]; + xfinger = &xfinger[c->hydro.count + space_extra_parts]; + gfinger = &gfinger[c->grav.count + space_extra_gparts]; + sfinger = &sfinger[c->stars.count + space_extra_sparts]; /* Add this cell to the list of local cells */ s->local_cells_top[s->nr_local_cells] = k; @@ -1328,8 +1334,12 @@ 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. */ + 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 @@ -1369,6 +1379,39 @@ 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; + for (int ind = 0; ind < num_cells; ind++) { + struct cell *c = &cells_top[ind]; + cell_reorder_extra_parts(c); + } +} + +/** + * @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, NULL); +} + /** * @brief #threadpool mapper function to sanitize the cells * @@ -2250,6 +2293,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]; @@ -2264,6 +2309,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]; @@ -2278,6 +2325,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]; @@ -3234,8 +3283,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 + * cretion */ 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); } /** diff --git a/src/space.h b/src/space.h index cdf83e4ada..4f1a39c287 100644 --- a/src/space.h +++ b/src/space.h @@ -72,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_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. @@ -282,6 +285,7 @@ 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, -- GitLab From b9d76ec4bc32e760931e083b3323ffca63a28451 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 12 Nov 2018 17:39:02 +0100 Subject: [PATCH 11/52] Ignore buffer particles in the initial test for identical positions. --- src/engine.c | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/engine.c b/src/engine.c index a69c13a161..9f2a994ea4 100644 --- a/src/engine.c +++ b/src/engine.c @@ -2773,6 +2773,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) @@ -2795,6 +2799,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) -- GitLab From 3fc4650fb5dd1a200168fb16ed9ed8f292caa351 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 12 Nov 2018 17:56:59 +0100 Subject: [PATCH 12/52] Modify the i/o routines to make sure we do not write out particles that do not exist. --- src/common_io.c | 7 +++++-- src/parallel_io.c | 9 ++++++--- src/serial_io.c | 9 ++++++--- src/single_io.c | 9 ++++++--- 4 files changed, 23 insertions(+), 11 deletions(-) diff --git a/src/common_io.c b/src/common_io.c index 5ea80f6519..46cf7077fe 100644 --- a/src/common_io.c +++ b/src/common_io.c @@ -811,7 +811,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]; @@ -845,7 +846,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++; @@ -879,6 +881,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/parallel_io.c b/src/parallel_io.c index 510b637c67..74b06f919d 100644 --- a/src/parallel_io.c +++ b/src/parallel_io.c @@ -1204,9 +1204,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/serial_io.c b/src/serial_io.c index 66d19b2291..fec7dffc24 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -775,9 +775,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 99f016809d..3ba58aa55e 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -639,9 +639,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; -- GitLab From b2a85c8e0bc22aa09f6501cc0697a3165805f2bf Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 12 Nov 2018 18:04:19 +0100 Subject: [PATCH 13/52] Make sure the extra particle ordering only takes place in top-level cells. --- src/cell.c | 3 +++ src/space.c | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/cell.c b/src/cell.c index 57795efcd3..bcf8b064b0 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3801,6 +3801,9 @@ void cell_reorder_extra_parts(struct cell *c) { const int count_real = c->hydro.count; const int count_total = count_real + space_extra_parts; + if (c->depth != 0) + error("This function should only be called on top-level cells!"); + int first_not_extra = count_real; /* Find extra particles */ diff --git a/src/space.c b/src/space.c index 6f78ed8674..9c67681854 100644 --- a/src/space.c +++ b/src/space.c @@ -3284,7 +3284,7 @@ void space_convert_quantities_mapper(void *restrict map_data, int count, struct xpart *restrict xparts = s->xparts + index; /* Loop over all the particles ignoring the extra buffer ones for on-the-fly - * cretion */ + * creation */ for (int k = 0; k < count; k++) if (parts[k].time_bin <= num_time_bins) hydro_convert_quantities(&parts[k], &xparts[k], cosmo, hydro_props); -- GitLab From ea909d292c443e4811c383cb3853bc1ce76a1c4f Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 12 Nov 2018 18:22:15 +0100 Subject: [PATCH 14/52] Crash explicitely when trying to regrid with fewer cells. --- src/space.c | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/src/space.c b/src/space.c index 9c67681854..d6895fb0fc 100644 --- a/src/space.c +++ b/src/space.c @@ -601,6 +601,18 @@ 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; @@ -632,8 +644,8 @@ void space_allocate_extras(struct space *s, int verbose) { /* 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_parts; - const size_t expected_num_extra_sparts = 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, @@ -724,6 +736,11 @@ void space_allocate_extras(struct space *s, int verbose) { s->nr_extra_parts = expected_num_extra_parts; } + /* Did we reduce the number of top-level cells ? */ + if (expected_num_extra_parts < s->nr_extra_parts) { + error("Un-handled case!"); + } + if (nr_gparts + expected_num_extra_gparts > size_gparts) { } if (nr_sparts + expected_num_extra_sparts > size_sparts) { -- GitLab From 9e2aadba0328da234f513b7de225646979a13329 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 12 Nov 2018 23:16:00 +0100 Subject: [PATCH 15/52] Also handle part with gpart friends. --- src/cell.c | 79 ++++++++++++++++++++++++++++++++++++++++-- src/cell.h | 6 ++-- src/engine.c | 4 ++- src/space.c | 98 ++++++++++++++++++++++++++++++++++++++-------------- 4 files changed, 155 insertions(+), 32 deletions(-) diff --git a/src/cell.c b/src/cell.c index bcf8b064b0..141c9a4079 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3794,7 +3794,7 @@ void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c, #endif } -void cell_reorder_extra_parts(struct cell *c) { +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; @@ -3816,12 +3816,87 @@ void cell_reorder_extra_parts(struct cell *c) { ++first_not_extra; } +#ifdef SWIFT_DEBUG_CHECKS if (first_not_extra >= count_total) 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) error("Need to handle gpart pointer!"); + if (parts[i].gpart) + parts[i].gpart->id_or_neg_offset = -(i + parts_offset); + } + } +} + +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; + const int count_total = count_real + space_extra_parts; + + if (c->depth != 0) + error("This function should only be called on 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_total) + 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); + } + } +} + +void cell_reorder_extra_gparts(struct cell *c, const ptrdiff_t sparts_offset) { + + struct gpart *gparts = c->grav.parts; + const int count_real = c->grav.count; + const int count_total = count_real + space_extra_parts; + + if (c->depth != 0) + error("This function should only be called on 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_total) + error("Looking for extra particles beyond this cell's range!"); +#endif + + /* Swap everything (including pointers) */ + memswap(&gparts[i], &gparts[first_not_extra], sizeof(struct spart)); + if (gparts[i].type == swift_type_gas) { + error("Need to handle this."); + } else if (gparts[i].type == swift_type_stars) { + error("Need to handle this."); + } } } } diff --git a/src/cell.h b/src/cell.h index 009e6522db..cf286f02f1 100644 --- a/src/cell.h +++ b/src/cell.h @@ -663,9 +663,9 @@ 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); -void cell_reorder_extra_parts(struct cell *c); -void cell_reorder_extra_gparts(struct cell *c); -void cell_reorder_extra_sparts(struct cell *c); +void cell_reorder_extra_parts(struct cell *c, const ptrdiff_t parts_offset); +void cell_reorder_extra_gparts(struct cell *c, const ptrdiff_t gparts_offset); +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/engine.c b/src/engine.c index 9f2a994ea4..9239d00c27 100644 --- a/src/engine.c +++ b/src/engine.c @@ -5143,9 +5143,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 diff --git a/src/space.c b/src/space.c index d6895fb0fc..9f313a3617 100644 --- a/src/space.c +++ b/src/space.c @@ -607,8 +607,8 @@ void space_regrid(struct space *s, int verbose) { * 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. + * 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? @@ -617,6 +617,11 @@ 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], @@ -650,7 +655,7 @@ void space_allocate_extras(struct space *s, int verbose) { 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 future particles.", + 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, @@ -707,6 +712,12 @@ void space_allocate_extras(struct space *s, int verbose) { 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 */ @@ -736,9 +747,8 @@ void space_allocate_extras(struct space *s, int verbose) { s->nr_extra_parts = expected_num_extra_parts; } - /* Did we reduce the number of top-level cells ? */ if (expected_num_extra_parts < s->nr_extra_parts) { - error("Un-handled case!"); + error("Reduction in top-level cells number not handled."); } if (nr_gparts + expected_num_extra_gparts > size_gparts) { @@ -1400,9 +1410,35 @@ 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); + 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, c->grav.parts - s->gparts); + } +} + +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); } } @@ -1426,7 +1462,17 @@ void space_reorder_extras(struct space *s, int verbose) { /* 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, NULL); + 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); } /** @@ -1538,13 +1584,13 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts, /* 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; @@ -1645,17 +1691,17 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts, /* 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 */ @@ -1753,14 +1799,14 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts, /* 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; -- GitLab From 6b07dd25544994abbdafe7edbdb57e6ade1a2cee Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Tue, 13 Nov 2018 12:03:41 +0100 Subject: [PATCH 16/52] Handle correctly the addition of spare star particles at rebuild time. --- src/cell.c | 2 +- src/space.c | 123 +++++++++++++++++++++++++++++++++++++++++++--------- 2 files changed, 104 insertions(+), 21 deletions(-) diff --git a/src/cell.c b/src/cell.c index 141c9a4079..67d72bf620 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3834,7 +3834,7 @@ 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; - const int count_total = count_real + space_extra_parts; + const int count_total = count_real + space_extra_sparts; if (c->depth != 0) error("This function should only be called on top-level cells!"); diff --git a/src/space.c b/src/space.c index 9f313a3617..98f93d0d0d 100644 --- a/src/space.c +++ b/src/space.c @@ -662,7 +662,15 @@ void space_allocate_extras(struct space *s, int verbose) { expected_num_extra_sparts); } - /* Do we have enough space for the extras (i.e. we haven't used up any) ? */ + 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 parts (i.e. we haven't used up any) ? + */ if (expected_num_extra_parts > s->nr_extra_parts) { /* Do we need to reallocate? */ @@ -747,13 +755,83 @@ void space_allocate_extras(struct space *s, int verbose) { s->nr_extra_parts = expected_num_extra_parts; } - if (expected_num_extra_parts < s->nr_extra_parts) { - error("Reduction in top-level cells number not handled."); - } - if (nr_gparts + expected_num_extra_gparts > size_gparts) { } + + /* Do we have enough space for the extra sparts (i.e. we haven't used up any) + * ? */ if (nr_sparts + expected_num_extra_sparts > size_sparts) { + + /* 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; + } + + /* 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; } } @@ -1590,12 +1668,12 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts, /* 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. */ @@ -1702,12 +1780,12 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts, 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. */ @@ -1795,6 +1873,11 @@ 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; @@ -1806,12 +1889,12 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts, /* 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. */ -- GitLab From 66a86e1059a71671ff5633d8f2457fa29eb550c4 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Tue, 13 Nov 2018 14:59:43 +0100 Subject: [PATCH 17/52] Also handle the creation of extra gpart at rebuild time. --- src/cell.c | 11 +++--- src/cell.h | 3 +- src/engine.c | 7 ++-- src/part.c | 5 ++- src/space.c | 109 +++++++++++++++++++++++++++++++++++++++++++++++++-- src/space.h | 4 +- 6 files changed, 122 insertions(+), 17 deletions(-) diff --git a/src/cell.c b/src/cell.c index 67d72bf620..4c5064e203 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3864,11 +3864,12 @@ void cell_reorder_extra_sparts(struct cell *c, const ptrdiff_t sparts_offset) { } } -void cell_reorder_extra_gparts(struct cell *c, const ptrdiff_t sparts_offset) { +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; - const int count_total = count_real + space_extra_parts; + const int count_total = count_real + space_extra_gparts; if (c->depth != 0) error("This function should only be called on top-level cells!"); @@ -3891,11 +3892,11 @@ void cell_reorder_extra_gparts(struct cell *c, const ptrdiff_t sparts_offset) { #endif /* Swap everything (including pointers) */ - memswap(&gparts[i], &gparts[first_not_extra], sizeof(struct spart)); + memswap(&gparts[i], &gparts[first_not_extra], sizeof(struct gpart)); if (gparts[i].type == swift_type_gas) { - error("Need to handle this."); + parts[-gparts[i].id_or_neg_offset].gpart = &gparts[i]; } else if (gparts[i].type == swift_type_stars) { - error("Need to handle this."); + sparts[-gparts[i].id_or_neg_offset].gpart = &gparts[i]; } } } diff --git a/src/cell.h b/src/cell.h index cf286f02f1..b6bdf1561b 100644 --- a/src/cell.h +++ b/src/cell.h @@ -664,7 +664,8 @@ void cell_convert_part_to_gpart(const struct engine *e, struct cell *c, void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c, struct spart *sp); void cell_reorder_extra_parts(struct cell *c, const ptrdiff_t parts_offset); -void cell_reorder_extra_gparts(struct cell *c, const ptrdiff_t gparts_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); diff --git a/src/engine.c b/src/engine.c index 9239d00c27..0530ae3fb9 100644 --- a/src/engine.c +++ b/src/engine.c @@ -1992,9 +1992,10 @@ void engine_rebuild(struct engine *e, int repartitioned, space_rebuild(e->s, repartitioned, e->verbose); /* 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); diff --git a/src/part.c b/src/part.c index 3a626e652c..ec3627d728 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/space.c b/src/space.c index 98f93d0d0d..88da0b5fed 100644 --- a/src/space.c +++ b/src/space.c @@ -669,10 +669,103 @@ void space_allocate_extras(struct space *s, int verbose) { 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) { @@ -711,6 +804,7 @@ void space_allocate_extras(struct space *s, int verbose) { 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 */ @@ -755,13 +849,12 @@ void space_allocate_extras(struct space *s, int verbose) { s->nr_extra_parts = expected_num_extra_parts; } - if (nr_gparts + expected_num_extra_gparts > size_gparts) { - } - /* Do we have enough space for the extra sparts (i.e. we haven't used up any) * ? */ if (nr_sparts + expected_num_extra_sparts > size_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) { @@ -790,6 +883,7 @@ void space_allocate_extras(struct space *s, int verbose) { ++i) { bzero(&s->sparts[i], sizeof(struct spart)); s->sparts[i].time_bin = time_bin_not_created; + s->sparts[i].id = -1; } /* Put the spare particles in their correct cell */ @@ -833,6 +927,13 @@ void space_allocate_extras(struct space *s, int verbose) { 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 } /** @@ -1504,7 +1605,7 @@ void space_reorder_extra_gparts_mapper(void *map_data, int num_cells, for (int ind = 0; ind < num_cells; ind++) { struct cell *c = &cells_top[ind]; - cell_reorder_extra_gparts(c, c->grav.parts - s->gparts); + cell_reorder_extra_gparts(c, s->parts, s->sparts); } } diff --git a/src/space.h b/src/space.h index 4f1a39c287..efb3ce407f 100644 --- a/src/space.h +++ b/src/space.h @@ -45,8 +45,8 @@ struct cosmology; #define space_splitsize_default 400 #define space_maxsize_default 8000000 #define space_extra_parts_default 10 -#define space_extra_gparts_default 0 -#define space_extra_sparts_default 0 +#define space_extra_gparts_default 30 +#define space_extra_sparts_default 20 #define space_expected_max_nr_strays_default 100 #define space_subsize_pair_hydro_default 256000000 #define space_subsize_self_hydro_default 32000 -- GitLab From b6662088a2e96b7b3b2045bd562bf8c8e86d44db Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Tue, 13 Nov 2018 15:15:31 +0100 Subject: [PATCH 18/52] Better documentation of new functions. --- src/cell.c | 36 ++++++++++++++++++++++++++++++------ 1 file changed, 30 insertions(+), 6 deletions(-) diff --git a/src/cell.c b/src/cell.c index 4c5064e203..3fcae363c7 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3794,6 +3794,14 @@ void cell_convert_spart_to_gpart(const struct engine *e, struct cell *c, #endif } +/** + * @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; @@ -3801,8 +3809,8 @@ void cell_reorder_extra_parts(struct cell *c, const ptrdiff_t parts_offset) { const int count_real = c->hydro.count; const int count_total = count_real + space_extra_parts; - if (c->depth != 0) - error("This function should only be called on top-level cells!"); + 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; @@ -3830,14 +3838,22 @@ void cell_reorder_extra_parts(struct cell *c, const ptrdiff_t 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; const int count_total = count_real + space_extra_sparts; - if (c->depth != 0) - error("This function should only be called on top-level cells!"); + 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; @@ -3864,6 +3880,14 @@ void cell_reorder_extra_sparts(struct cell *c, const ptrdiff_t sparts_offset) { } } +/** + * @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) { @@ -3871,8 +3895,8 @@ void cell_reorder_extra_gparts(struct cell *c, struct part *parts, const int count_real = c->grav.count; const int count_total = count_real + space_extra_gparts; - if (c->depth != 0) - error("This function should only be called on top-level cells!"); + 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; -- GitLab From 42cc0597f4baf03d0c3d87971ed6affdaf2d08cf Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Tue, 13 Nov 2018 15:37:34 +0100 Subject: [PATCH 19/52] Add the ability to have i/o conversion functions for star particles as well. --- src/common_io.c | 40 +++++++++++++++ src/io_properties.h | 121 +++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 159 insertions(+), 2 deletions(-) diff --git a/src/common_io.c b/src/common_io.c index 46cf7077fe..a9845b990a 100644 --- a/src/common_io.c +++ b/src/common_io.c @@ -522,6 +522,46 @@ void io_convert_gpart_d_mapper(void* restrict temp, int N, props.convert_gpart_d(e, gparts + delta + i, &temp_d[i * dim]); } +/** + * @brief Mapper function to copy #spart into a buffer of floats using a + * conversion function. + */ +void io_convert_spart_f_mapper(void* restrict temp, int N, + void* restrict extra_data) { + + const struct io_props props = *((const struct io_props*)extra_data); + const struct spart* restrict sparts = props.sparts; + const struct engine* e = props.e; + const size_t dim = props.dimension; + + /* How far are we with this chunk? */ + float* restrict temp_f = (float*)temp; + const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim; + + for (int i = 0; i < N; i++) + props.convert_spart_f(e, sparts + delta + i, &temp_f[i * dim]); +} + +/** + * @brief Mapper function to copy #spart into a buffer of doubles using a + * conversion function. + */ +void io_convert_spart_d_mapper(void* restrict temp, int N, + void* restrict extra_data) { + + const struct io_props props = *((const struct io_props*)extra_data); + const struct spart* restrict sparts = props.sparts; + const struct engine* e = props.e; + const size_t dim = props.dimension; + + /* How far are we with this chunk? */ + double* restrict temp_d = (double*)temp; + const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim; + + for (int i = 0; i < N; i++) + props.convert_spart_d(e, sparts + delta + i, &temp_d[i * dim]); +} + /** * @brief Copy the particle data into a temporary buffer ready for i/o. * diff --git a/src/io_properties.h b/src/io_properties.h index 037d32338f..9e948fc399 100644 --- a/src/io_properties.h +++ b/src/io_properties.h @@ -47,6 +47,10 @@ typedef void (*conversion_func_gpart_float)(const struct engine*, const struct gpart*, float*); typedef void (*conversion_func_gpart_double)(const struct engine*, const struct gpart*, double*); +typedef void (*conversion_func_spart_float)(const struct engine*, + const struct spart*, float*); +typedef void (*conversion_func_spart_double)(const struct engine*, + const struct spart*, double*); /** * @brief The properties of a given dataset for i/o @@ -86,6 +90,7 @@ struct io_props { const struct part* parts; const struct xpart* xparts; const struct gpart* gparts; + const struct spart* sparts; /* Are we converting? */ int conversion; @@ -97,6 +102,10 @@ struct io_props { /* Conversion function for gpart */ conversion_func_gpart_float convert_gpart_f; conversion_func_gpart_double convert_gpart_d; + + /* Conversion function for spart */ + conversion_func_spart_float convert_spart_f; + conversion_func_spart_double convert_spart_d; }; /** @@ -134,11 +143,14 @@ INLINE static struct io_props io_make_input_field_( r.parts = NULL; r.xparts = NULL; r.gparts = NULL; + r.sparts = NULL; r.conversion = 0; r.convert_part_f = NULL; r.convert_part_d = NULL; r.convert_gpart_f = NULL; r.convert_gpart_d = NULL; + r.convert_spart_f = NULL; + r.convert_spart_d = NULL; return r; } @@ -175,11 +187,14 @@ INLINE static struct io_props io_make_output_field_( r.partSize = partSize; r.parts = NULL; r.gparts = NULL; + r.sparts = NULL; r.conversion = 0; r.convert_part_f = NULL; r.convert_part_d = NULL; r.convert_gpart_f = NULL; r.convert_gpart_d = NULL; + r.convert_spart_f = NULL; + r.convert_spart_d = NULL; return r; } @@ -223,11 +238,14 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT( r.parts = parts; r.xparts = xparts; r.gparts = NULL; + r.sparts = NULL; r.conversion = 1; r.convert_part_f = functionPtr; r.convert_part_d = NULL; r.convert_gpart_f = NULL; r.convert_gpart_d = NULL; + r.convert_spart_f = NULL; + r.convert_spart_d = NULL; return r; } @@ -242,7 +260,7 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT( * @param partSize The size in byte of the particle * @param parts The particle array * @param xparts The xparticle array - * @param functionPtr The function used to convert a particle to a float + * @param functionPtr The function used to convert a particle to a double * * Do not call this function directly. Use the macro defined above. */ @@ -263,11 +281,14 @@ INLINE static struct io_props io_make_output_field_convert_part_DOUBLE( r.parts = parts; r.xparts = xparts; r.gparts = NULL; + r.sparts = NULL; r.conversion = 1; r.convert_part_f = NULL; r.convert_part_d = functionPtr; r.convert_gpart_f = NULL; r.convert_gpart_d = NULL; + r.convert_spart_f = NULL; + r.convert_spart_d = NULL; return r; } @@ -309,11 +330,14 @@ INLINE static struct io_props io_make_output_field_convert_gpart_FLOAT( r.parts = NULL; r.xparts = NULL; r.gparts = gparts; + r.sparts = NULL; r.conversion = 1; r.convert_part_f = NULL; r.convert_part_d = NULL; r.convert_gpart_f = functionPtr; r.convert_gpart_d = NULL; + r.convert_spart_f = NULL; + r.convert_spart_d = NULL; return r; } @@ -327,7 +351,7 @@ INLINE static struct io_props io_make_output_field_convert_gpart_FLOAT( * @param units The units of the dataset * @param gpartSize The size in byte of the particle * @param gparts The particle array - * @param functionPtr The function used to convert a g-particle to a float + * @param functionPtr The function used to convert a g-particle to a double * * Do not call this function directly. Use the macro defined above. */ @@ -347,11 +371,104 @@ INLINE static struct io_props io_make_output_field_convert_gpart_DOUBLE( r.parts = NULL; r.xparts = NULL; r.gparts = gparts; + r.sparts = NULL; r.conversion = 1; r.convert_part_f = NULL; r.convert_part_d = NULL; r.convert_gpart_f = NULL; r.convert_gpart_d = functionPtr; + r.convert_spart_f = NULL; + r.convert_spart_d = NULL; + + return r; +} + +/** + * @brief Constructs an #io_props (with conversion) from its parameters + */ +#define io_make_output_field_convert_spart(name, type, dim, units, spart, \ + convert) \ + io_make_output_field_convert_spart_##type(name, type, dim, units, \ + sizeof(spart[0]), spart, convert) + +/** + * @brief Construct an #io_props from its parameters + * + * @param name Name of the field to read + * @param type The type of the data + * @param dimension Dataset dimension (1D, 3D, ...) + * @param units The units of the dataset + * @param spartSize The size in byte of the particle + * @param sparts The particle array + * @param functionPtr The function used to convert a g-particle to a float + * + * Do not call this function directly. Use the macro defined above. + */ +INLINE static struct io_props io_make_output_field_convert_spart_FLOAT( + const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension, + enum unit_conversion_factor units, size_t spartSize, + const struct spart* sparts, conversion_func_spart_float functionPtr) { + + struct io_props r; + strcpy(r.name, name); + r.type = type; + r.dimension = dimension; + r.importance = UNUSED; + r.units = units; + r.field = NULL; + r.partSize = spartSize; + r.parts = NULL; + r.xparts = NULL; + r.gparts = NULL; + r.sparts = sparts; + r.conversion = 1; + r.convert_part_f = NULL; + r.convert_part_d = NULL; + r.convert_gpart_f = NULL; + r.convert_gpart_d = NULL; + r.convert_spart_f = functionPtr; + r.convert_spart_d = NULL; + + return r; +} + +/** + * @brief Construct an #io_props from its parameters + * + * @param name Name of the field to read + * @param type The type of the data + * @param dimension Dataset dimension (1D, 3D, ...) + * @param units The units of the dataset + * @param spartSize The size in byte of the particle + * @param sparts The particle array + * @param functionPtr The function used to convert a s-particle to a double + * + * Do not call this function directly. Use the macro defined above. + */ +INLINE static struct io_props io_make_output_field_convert_spart_DOUBLE( + const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension, + enum unit_conversion_factor units, size_t spartSize, + const struct spart* sparts, conversion_func_spart_double functionPtr) { + + struct io_props r; + strcpy(r.name, name); + r.type = type; + r.dimension = dimension; + r.importance = UNUSED; + r.units = units; + r.field = NULL; + r.partSize = spartSize; + r.parts = NULL; + r.xparts = NULL; + r.gparts = NULL; + r.sparts = sparts; + r.conversion = 1; + r.convert_part_f = NULL; + r.convert_part_d = NULL; + r.convert_gpart_f = NULL; + r.convert_gpart_d = NULL; + r.convert_spart_f = NULL; + r.convert_spart_d = functionPtr; return r; } -- GitLab From 17d568e21638cea5e5cd2c5bf0345601181b35a5 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Tue, 13 Nov 2018 15:41:30 +0100 Subject: [PATCH 20/52] Add the ability to have i/o conversion functions for star particles as well. --- src/common_io.c | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/common_io.c b/src/common_io.c index a9845b990a..24e74014fd 100644 --- a/src/common_io.c +++ b/src/common_io.c @@ -643,6 +643,30 @@ void io_copy_temp_buffer(void* temp, const struct engine* e, io_convert_gpart_d_mapper, temp_d, N, copySize, 0, (void*)&props); + } else if (props.convert_spart_f != NULL) { + + /* Prepare some parameters */ + float* temp_f = (float*)temp; + props.start_temp_f = (float*)temp; + props.e = e; + + /* Copy the whole thing into a buffer */ + threadpool_map((struct threadpool*)&e->threadpool, + io_convert_spart_f_mapper, temp_f, N, copySize, 0, + (void*)&props); + + } else if (props.convert_spart_d != NULL) { + + /* Prepare some parameters */ + double* temp_d = (double*)temp; + props.start_temp_d = (double*)temp; + props.e = e; + + /* Copy the whole thing into a buffer */ + threadpool_map((struct threadpool*)&e->threadpool, + io_convert_spart_d_mapper, temp_d, N, copySize, 0, + (void*)&props); + } else { error("Missing conversion function"); } -- GitLab From f26d1107c178a5f6c7d17f3bdf96970b74faa9cc Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Tue, 13 Nov 2018 15:44:15 +0100 Subject: [PATCH 21/52] Temporarily write the time-bins to the snapshots to make sure we only write real particles. --- src/gravity/Default/gravity_io.h | 13 ++++++++++++- src/hydro/Gadget2/hydro_io.h | 15 ++++++++++++++- src/stars/Default/stars_io.h | 16 ++++++++++++++-- 3 files changed, 40 insertions(+), 4 deletions(-) diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h index 1ba3899e7e..0c3c643102 100644 --- a/src/gravity/Default/gravity_io.h +++ b/src/gravity/Default/gravity_io.h @@ -92,6 +92,15 @@ INLINE static void darkmatter_read_particles(struct gpart* gparts, UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset); } +INLINE static void test_time_bin_gpart(const struct engine* e, + const struct gpart* gp, float* ret) { + + if (gp->time_bin >= time_bin_inhibited) + error("Writing inhibited or extra particle time_bin=%d", gp->time_bin); + + *ret = gp->time_bin; +} + /** * @brief Specifies which g-particle fields to write to a dataset * @@ -104,7 +113,7 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts, int* num_fields) { /* Say how much we want to write */ - *num_fields = 4; + *num_fields = 5; /* List what we want to write */ list[0] = io_make_output_field_convert_gpart( @@ -115,6 +124,8 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts, io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, gparts, mass); list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset); + list[4] = io_make_output_field_convert_gpart( + "TimeBin", FLOAT, 1, UNIT_CONV_NO_UNITS, gparts, test_time_bin_gpart); } #endif /* SWIFT_DEFAULT_GRAVITY_IO_H */ diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h index ec7d34f7ad..d776951767 100644 --- a/src/hydro/Gadget2/hydro_io.h +++ b/src/hydro/Gadget2/hydro_io.h @@ -128,6 +128,16 @@ INLINE static void convert_part_potential(const struct engine* e, ret[0] = 0.f; } +INLINE static void test_time_bin_part(const struct engine* e, + const struct part* p, + const struct xpart* xp, float* ret) { + + if (p->time_bin >= time_bin_inhibited) + error("Writing inhibited or extra particle time_bin=%d", p->time_bin); + + *ret = p->time_bin; +} + /** * @brief Specifies which particle fields to write to a dataset * @@ -141,7 +151,7 @@ INLINE static void hydro_write_particles(const struct part* parts, struct io_props* list, int* num_fields) { - *num_fields = 10; + *num_fields = 11; /* List what we want to write */ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3, @@ -168,6 +178,9 @@ INLINE static void hydro_write_particles(const struct part* parts, list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1, UNIT_CONV_POTENTIAL, parts, xparts, convert_part_potential); + list[10] = + io_make_output_field_convert_part("TimeBin", FLOAT, 1, UNIT_CONV_NO_UNITS, + parts, xparts, test_time_bin_part); #ifdef DEBUG_INTERACTIONS_SPH diff --git a/src/stars/Default/stars_io.h b/src/stars/Default/stars_io.h index a6c2768f71..82d06d292f 100644 --- a/src/stars/Default/stars_io.h +++ b/src/stars/Default/stars_io.h @@ -22,6 +22,15 @@ #include "io_properties.h" #include "stars_part.h" +INLINE static void test_time_bin_spart(const struct engine *e, + const struct spart *sp, float *ret) { + + if (sp->time_bin >= time_bin_inhibited) + error("Writing inhibited or extra particle time_bin=%d", sp->time_bin); + + *ret = sp->time_bin; +} + /** * @brief Specifies which s-particle fields to read from a dataset * @@ -60,8 +69,8 @@ INLINE static void stars_write_particles(const struct spart *sparts, struct io_props *list, int *num_fields) { - /* Say how much we want to read */ - *num_fields = 5; + /* Say how much we want to write */ + *num_fields = 6; /* List what we want to read */ list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, @@ -74,6 +83,9 @@ INLINE static void stars_write_particles(const struct spart *sparts, sparts, id); list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH, sparts, h); + + list[5] = io_make_output_field_convert_spart( + "TimeBin", FLOAT, 1, UNIT_CONV_NO_UNITS, sparts, test_time_bin_spart); } /** -- GitLab From fc928650db6cf17191007ae2e697268d0a5947b2 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Tue, 13 Nov 2018 15:58:21 +0100 Subject: [PATCH 22/52] Use the correct array of star particles (the selected ones) when writing a snaphot --- src/parallel_io.c | 2 +- src/serial_io.c | 2 +- src/single_io.c | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/parallel_io.c b/src/parallel_io.c index 74b06f919d..7e34152f10 100644 --- a/src/parallel_io.c +++ b/src/parallel_io.c @@ -1456,7 +1456,7 @@ void write_output_parallel(struct engine* e, const char* baseName, Nstars_written); /* Select the fields to write */ - stars_write_particles(sparts, list, &num_fields); + stars_write_particles(sparts_written, list, &num_fields); } } break; diff --git a/src/serial_io.c b/src/serial_io.c index fec7dffc24..48d9838954 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -1132,7 +1132,7 @@ void write_output_serial(struct engine* e, const char* baseName, Nstars_written); /* Select the fields to write */ - stars_write_particles(sparts, list, &num_fields); + stars_write_particles(sparts_written, list, &num_fields); } } break; diff --git a/src/single_io.c b/src/single_io.c index 3ba58aa55e..780978fe94 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -952,7 +952,7 @@ void write_output_single(struct engine* e, const char* baseName, Nstars_written); /* Select the fields to write */ - stars_write_particles(sparts, list, &num_fields); + stars_write_particles(sparts_written, list, &num_fields); } } break; -- GitLab From 48f5ecb709098a144375005c0c62c9ef5b18ae98 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Tue, 13 Nov 2018 17:39:22 +0100 Subject: [PATCH 23/52] Make the star formation task lock it's part of the tree. --- src/task.c | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/task.c b/src/task.c index 9d4be3aaa5..6bc29730ed 100644 --- a/src/task.c +++ b/src/task.c @@ -380,6 +380,11 @@ void task_unlock(struct task *t) { cell_munlocktree(cj); break; + case task_type_star_formation: + cell_unlocktree(ci); + cell_sunlocktree(ci); + cell_gunlocktree(ci); + default: break; } @@ -518,6 +523,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; -- GitLab From 4feef16693cc6e4f7444840dc6c78e26e0044200 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Tue, 13 Nov 2018 18:36:46 +0100 Subject: [PATCH 24/52] Added counter of total number of particles in the cells not just counter of the real ones. --- src/cell.c | 20 +++++++++++++++++++- src/cell.h | 9 ++++++++- src/runner.c | 1 + src/space.c | 25 +++++++++++++++++++++---- src/space.h | 4 ++-- src/task.c | 1 + 6 files changed, 52 insertions(+), 8 deletions(-) diff --git a/src/cell.c b/src/cell.c index 819b039d14..4bc9ca251a 100644 --- a/src/cell.c +++ b/src/cell.c @@ -967,6 +967,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]]; } @@ -1084,6 +1085,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]]; } @@ -1146,6 +1148,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]]; } } @@ -3645,7 +3648,22 @@ void cell_check_timesteps(struct cell *c) { #endif } -void cell_add_spart(const struct engine *e, struct cell *c) {} +void cell_add_spart(const struct engine *e, struct cell *c) { + + if (c->nodeID != engine_rank) error("Adding spart on a foreign node"); + + /* Get the top-level this leaf cell is in */ + struct cell *top = c; + while (top->parent != NULL) top = top->parent; + +#ifdef SWIFT_DEBUG_CHECKS + if (top->depth != 0) error("Cell-linking issue"); +#endif + + /* Are there any extra particles left? */ + if (top->stars.count == top->stars.count_total) + error("We ran out of star particles!"); +} /** * @brief "Remove" a gas particle from the calculation. diff --git a/src/cell.h b/src/cell.h index 1f9dc110f2..f3fc61308c 100644 --- a/src/cell.h +++ b/src/cell.h @@ -312,6 +312,9 @@ struct cell { /*! 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; @@ -339,7 +342,7 @@ struct cell { /*! Do any of this cell's sub-cells need to be sorted? */ char do_sub_sort; - #ifdef SWIFT_DEBUG_CHECKS +#ifdef SWIFT_DEBUG_CHECKS /*! Last (integer) time the cell's sort arrays were updated. */ integertime_t ti_sort; @@ -416,6 +419,9 @@ struct cell { /*! 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; @@ -659,6 +665,7 @@ void cell_activate_sorts(struct cell *c, int sid, struct scheduler *s); void cell_clear_drift_flags(struct cell *c, void *data); void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data); int cell_has_tasks(struct cell *c); +void cell_add_spart(const struct engine *e, struct cell *c); void cell_remove_part(const struct engine *e, struct cell *c, struct part *p, struct xpart *xp); void cell_remove_gpart(const struct engine *e, struct cell *c, diff --git a/src/runner.c b/src/runner.c index 0026d833b9..4d891063bc 100644 --- a/src/runner.c +++ b/src/runner.c @@ -525,6 +525,7 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { if (rho > 1.5e7 && e->step > 2) { message("Removing particle id=%lld rho=%e", p->id, rho); cell_convert_part_to_gpart(e, c, p, xp); + cell_add_spart(e, c); } } } diff --git a/src/space.c b/src/space.c index 3b3d1927fc..323e9d3a01 100644 --- a/src/space.c +++ b/src/space.c @@ -194,12 +194,15 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements, c->stars.dx_max_part = 0.f; c->hydro.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; @@ -1522,10 +1525,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 + space_extra_parts]; - xfinger = &xfinger[c->hydro.count + space_extra_parts]; - gfinger = &gfinger[c->grav.count + space_extra_gparts]; - sfinger = &sfinger[c->stars.count + space_extra_sparts]; + + 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; @@ -2616,6 +2624,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; @@ -2821,6 +2832,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 @@ -2839,6 +2852,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 @@ -2849,6 +2864,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 diff --git a/src/space.h b/src/space.h index 65a1f4fb90..c06ae26e1a 100644 --- a/src/space.h +++ b/src/space.h @@ -44,8 +44,8 @@ struct cosmology; #define space_cellallocchunk 1000 #define space_splitsize_default 400 #define space_maxsize_default 8000000 -#define space_extra_parts_default 10 -#define space_extra_gparts_default 30 +#define space_extra_parts_default 0 +#define space_extra_gparts_default 0 #define space_extra_sparts_default 20 #define space_expected_max_nr_strays_default 100 #define space_subsize_pair_hydro_default 256000000 diff --git a/src/task.c b/src/task.c index 6bc29730ed..52c62c67f1 100644 --- a/src/task.c +++ b/src/task.c @@ -384,6 +384,7 @@ void task_unlock(struct task *t) { cell_unlocktree(ci); cell_sunlocktree(ci); cell_gunlocktree(ci); + break; default: break; -- GitLab From 3fa75e9fb8f420bbf84b53143f6acb7cfa079a45 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 14 Nov 2018 17:44:43 +0100 Subject: [PATCH 25/52] Version that can survive the addition of one particle. --- src/cell.c | 156 +++++++++++++++++++++++++++++++++++++++++++++++++--- src/cell.h | 4 +- src/space.c | 2 +- 3 files changed, 152 insertions(+), 10 deletions(-) diff --git a/src/cell.c b/src/cell.c index 4bc9ca251a..925b6422a8 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3648,14 +3648,95 @@ void cell_check_timesteps(struct cell *c) { #endif } -void cell_add_spart(const struct engine *e, struct cell *c) { +void cell_check_spart_pos_and_time(const struct cell* c) { + +#ifdef SWIFT_DEBUG_CHECKS + + /* Recurse */ + if (c->split) { + for(int k = 0; k < 8; ++k) + if (c->progeny[k] != NULL) + cell_check_spart_pos_and_time(c->progeny[k]); + } + + 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]) || + (sp->x[1] < c->loc[1]) || + (sp->x[2] < c->loc[2]) || + (sp->x[0] >= c->loc[0] + c->width[0]) || + (sp->x[1] >= c->loc[1] + c->width[1]) || + (sp->x[2] >= c->loc[2] + c->width[2])) + error("spart not in its cell!"); + } + +#else + error("Calling a degugging function outside debugging mode."); +#endif + +} + +void cell_recursively_shift_sparts(struct cell *c, struct cell *const target_cell, struct spart const* start, + const int depth, struct spart *temp, + const int direct_progeny){ + + if( c->stars.parts < start || (c->stars.parts == start && c < target_cell)) { + + if(c->stars.parts + c->stars.count > start) + c->stars.parts++; + + /* Abort if we are in a cell before the point where we want to start + the memory shift */ + return; + } + + message("cellID=%d", c->cellID); + + if (c->split) { + + int first_progeny = -1; + + for(int k = 0; k < 8; ++k) + if (c->progeny[k] != NULL) { + /* if(first_progeny == -1) */ + /* first_progeny = k; */ + + cell_recursively_shift_sparts(c->progeny[k], target_cell, start, depth, temp, k == first_progeny); + } + } + else { + //if (c->stars.count > 0) { + //memswap(temp, &c->stars.parts[0], sizeof(struct spart)); + memswap(temp, &c->stars.parts[c->stars.count], sizeof(struct spart)); + //} + //else + //memswap(temp, &c->stars.parts[0], sizeof(struct spart)); + } + + /* If this is a direct parent, add one element */ + //if((c->stars.parts < start && c->stars.parts + c->stars.count >= start) || direct_progeny) + /* if(direct_progeny) */ + /* c->stars.count++; */ + //else if (c->stars.parts == start) + if(c->stars.parts >= start) + c->stars.parts++; +} + +void cell_add_spart(const struct engine *e, struct cell *const c) { if (c->nodeID != engine_rank) error("Adding spart on a foreign node"); /* Get the top-level this leaf cell is in */ struct cell *top = c; - while (top->parent != NULL) top = top->parent; - + while (top->parent != NULL) { + message("depth=%d cellID=%d size=%d part=%p", top->depth, top->cellID, top->stars.count, top->stars.parts); + top = top->parent; + } + message("depth=%d cellID=%d size=%d part=%p", top->depth, top->cellID, top->stars.count, top->stars.parts); + #ifdef SWIFT_DEBUG_CHECKS if (top->depth != 0) error("Cell-linking issue"); #endif @@ -3663,6 +3744,59 @@ void cell_add_spart(const struct engine *e, struct cell *c) { /* Are there any extra particles left? */ if (top->stars.count == top->stars.count_total) error("We ran out of star particles!"); + + message("cellID=%d top->size=%d c->size=%d c->depth=%d", c->cellID, top->stars.count, c->stars.count, c->depth); + + message("BEFORE: size=%d ptr=%p", c->parent->progeny[3]->progeny[4]->stars.count, c->parent->progeny[3]->progeny[4]->stars.parts); + + /* Recursively shift all the stars to get a free spot at the start of the current cell*/ + cell_recursively_shift_sparts(top, c, c->stars.parts, c->depth, &top->stars.parts[top->stars.count], 1); + c->stars.count++; + c->stars.parts--; + + top = c; + while (top->parent != NULL) { + top = top->parent; + top->stars.count++; + top->stars.parts--; + } + + message("AFTER: size=%d ptr=%p", c->parent->progeny[3]->progeny[4]->stars.count, c->parent->progeny[3]->progeny[4]->stars.parts); + + struct cell *top_new = c; + while (top_new->parent != NULL) { + message("depth=%d cellID=%d size=%d part=%p loc=[%e %e %e] w=[%e %e %e]", + top_new->depth, top_new->cellID, top_new->stars.count, top_new->stars.parts, + top_new->loc[0], top_new->loc[1], top_new->loc[2], + top_new->width[0], top_new->width[1], top_new->width[2]); + top_new = top_new->parent; + } + message("depth=%d cellID=%d size=%d part=%p loc=[%e %e %e] w=[%e %e %e]", + top_new->depth, top_new->cellID, top_new->stars.count, top_new->stars.parts, + top_new->loc[0], top_new->loc[1], top_new->loc[2], + top_new->width[0], top_new->width[1], top_new->width[2]); + + /* 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)); + message("sp->id=%lld sp->x=[%e %e %e]", sp->id, sp->x[0], sp->x[1], sp->x[2]); + message("c->loc=[%e %e %e], c->width=[%e %e %e]", c->loc[0], c->loc[1], c->loc[2], + c->width[0], c->width[1], c->width[2]); + + /* 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]; + + message("sp->id=%lld sp->x=[%e %e %e]", sp->id, sp->x[0], sp->x[1], sp->x[2]); + +#ifdef SWIFT_DEBUG_CHECKS + sp->ti_drift = e->ti_current; +#endif + +#ifdef SWIFT_DEBUG_CHECKS + cell_check_spart_pos_and_time(top); +#endif } /** @@ -3757,9 +3891,11 @@ void cell_remove_spart(const struct engine *e, struct cell *c, * @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. */ -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) @@ -3784,6 +3920,8 @@ 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; } /** @@ -3795,9 +3933,11 @@ void cell_convert_part_to_gpart(const struct engine *e, struct cell *c, * @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. */ -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) @@ -3822,6 +3962,8 @@ 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; } /** diff --git a/src/cell.h b/src/cell.h index f3fc61308c..a78aeb5bec 100644 --- a/src/cell.h +++ b/src/cell.h @@ -672,9 +672,9 @@ 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 gpart* 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 gpart* cell_convert_spart_to_gpart(const struct engine *e, struct cell *c, struct spart *sp); void cell_reorder_extra_parts(struct cell *c, const ptrdiff_t parts_offset); void cell_reorder_extra_gparts(struct cell *c, struct part *parts, diff --git a/src/space.c b/src/space.c index 323e9d3a01..dce5fff8fd 100644 --- a/src/space.c +++ b/src/space.c @@ -894,7 +894,7 @@ void space_allocate_extras(struct space *s, int verbose) { ++i) { bzero(&s->sparts[i], sizeof(struct spart)); s->sparts[i].time_bin = time_bin_not_created; - s->sparts[i].id = -1; + s->sparts[i].id = -42; } /* Put the spare particles in their correct cell */ -- GitLab From e44d0632b706d7c65a077fbd90c9f3b4503ffa51 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 14 Nov 2018 17:46:37 +0100 Subject: [PATCH 26/52] Version that can survive the addition of one particle. --- src/cell.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cell.c b/src/cell.c index 925b6422a8..5873616d74 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3685,8 +3685,8 @@ void cell_recursively_shift_sparts(struct cell *c, struct cell *const target_cel if( c->stars.parts < start || (c->stars.parts == start && c < target_cell)) { - if(c->stars.parts + c->stars.count > start) - c->stars.parts++; + /* if(c->stars.parts + c->stars.count > start) */ + /* c->stars.parts++; */ /* Abort if we are in a cell before the point where we want to start the memory shift */ -- GitLab From 6d6b9bc1ef7ad607eed907f92e5bc02d94dfd57a Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 14 Nov 2018 17:46:46 +0100 Subject: [PATCH 27/52] Revert "Version that can survive the addition of one particle." This reverts commit e44d0632b706d7c65a077fbd90c9f3b4503ffa51. --- src/cell.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cell.c b/src/cell.c index 5873616d74..925b6422a8 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3685,8 +3685,8 @@ void cell_recursively_shift_sparts(struct cell *c, struct cell *const target_cel if( c->stars.parts < start || (c->stars.parts == start && c < target_cell)) { - /* if(c->stars.parts + c->stars.count > start) */ - /* c->stars.parts++; */ + if(c->stars.parts + c->stars.count > start) + c->stars.parts++; /* Abort if we are in a cell before the point where we want to start the memory shift */ -- GitLab From de3421cef5c0b6730f1e99da9f1208882eb86b6f Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 14 Nov 2018 19:04:00 +0100 Subject: [PATCH 28/52] Can now survive 2 full steps. --- src/cell.c | 113 ++++++++++++++++++--------------------------------- src/runner.c | 4 +- src/space.h | 2 +- 3 files changed, 42 insertions(+), 77 deletions(-) diff --git a/src/cell.c b/src/cell.c index 925b6422a8..a7c9f33cad 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3648,7 +3648,7 @@ void cell_check_timesteps(struct cell *c) { #endif } -void cell_check_spart_pos_and_time(const struct cell* c) { +void cell_check_spart_pos(const struct cell* c) { #ifdef SWIFT_DEBUG_CHECKS @@ -3656,7 +3656,7 @@ void cell_check_spart_pos_and_time(const struct cell* c) { if (c->split) { for(int k = 0; k < 8; ++k) if (c->progeny[k] != NULL) - cell_check_spart_pos_and_time(c->progeny[k]); + cell_check_spart_pos(c->progeny[k]); } struct spart *sparts = c->stars.parts; @@ -3683,44 +3683,13 @@ void cell_recursively_shift_sparts(struct cell *c, struct cell *const target_cel const int depth, struct spart *temp, const int direct_progeny){ - if( c->stars.parts < start || (c->stars.parts == start && c < target_cell)) { - - if(c->stars.parts + c->stars.count > start) - c->stars.parts++; - - /* Abort if we are in a cell before the point where we want to start - the memory shift */ - return; - } - - message("cellID=%d", c->cellID); - if (c->split) { - - int first_progeny = -1; - - for(int k = 0; k < 8; ++k) - if (c->progeny[k] != NULL) { - /* if(first_progeny == -1) */ - /* first_progeny = k; */ - - cell_recursively_shift_sparts(c->progeny[k], target_cell, start, depth, temp, k == first_progeny); - } - } - else { - //if (c->stars.count > 0) { - //memswap(temp, &c->stars.parts[0], sizeof(struct spart)); - memswap(temp, &c->stars.parts[c->stars.count], sizeof(struct spart)); - //} - //else - //memswap(temp, &c->stars.parts[0], sizeof(struct spart)); + for(int k = 0; k < 8; ++k) { + if(c->progeny[k] != NULL) + cell_recursively_shift_sparts(c->progeny[k], target_cell, start, depth, temp, 0); + } } - /* If this is a direct parent, add one element */ - //if((c->stars.parts < start && c->stars.parts + c->stars.count >= start) || direct_progeny) - /* if(direct_progeny) */ - /* c->stars.count++; */ - //else if (c->stars.parts == start) if(c->stars.parts >= start) c->stars.parts++; } @@ -3728,74 +3697,70 @@ void cell_recursively_shift_sparts(struct cell *c, struct cell *const target_cel void cell_add_spart(const struct engine *e, struct cell *const c) { if (c->nodeID != engine_rank) error("Adding spart on a foreign node"); - + /* Get the top-level this leaf cell is in */ struct cell *top = c; while (top->parent != NULL) { - message("depth=%d cellID=%d size=%d part=%p", top->depth, top->cellID, top->stars.count, top->stars.parts); top = top->parent; } - message("depth=%d cellID=%d size=%d part=%p", top->depth, top->cellID, top->stars.count, top->stars.parts); #ifdef SWIFT_DEBUG_CHECKS if (top->depth != 0) error("Cell-linking issue"); #endif - + /* Are there any extra particles left? */ if (top->stars.count == top->stars.count_total) error("We ran out of star particles!"); - - message("cellID=%d top->size=%d c->size=%d c->depth=%d", c->cellID, top->stars.count, c->stars.count, c->depth); - - message("BEFORE: size=%d ptr=%p", c->parent->progeny[3]->progeny[4]->stars.count, c->parent->progeny[3]->progeny[4]->stars.parts); - - /* Recursively shift all the stars to get a free spot at the start of the current cell*/ - cell_recursively_shift_sparts(top, c, c->stars.parts, c->depth, &top->stars.parts[top->stars.count], 1); - c->stars.count++; - c->stars.parts--; - - top = c; - while (top->parent != NULL) { - top = top->parent; - top->stars.count++; - top->stars.parts--; - } - message("AFTER: size=%d ptr=%p", c->parent->progeny[3]->progeny[4]->stars.count, c->parent->progeny[3]->progeny[4]->stars.parts); + /* Copy the particles to get a free space */ + const int n_copy = &top->stars.parts[top->stars.count] - c->stars.parts; + + if(n_copy > 0) { + struct spart *temp = NULL; + if (posix_memalign((void**)&temp, spart_align, + n_copy * sizeof(struct spart)) != 0) + error("Impossible to allocate temp buffer"); + + memcpy(temp, c->stars.parts, n_copy * sizeof(struct spart)); + memcpy(c->stars.parts + 1, temp, n_copy * sizeof(struct spart)); + free(temp); + + message("count=%d n_copy=%d", c->stars.count, n_copy); - struct cell *top_new = c; - while (top_new->parent != NULL) { - message("depth=%d cellID=%d size=%d part=%p loc=[%e %e %e] w=[%e %e %e]", - top_new->depth, top_new->cellID, top_new->stars.count, top_new->stars.parts, - top_new->loc[0], top_new->loc[1], top_new->loc[2], - top_new->width[0], top_new->width[1], top_new->width[2]); - top_new = top_new->parent; + /* Recursively shift all the stars to get a free spot at the start of the current cell*/ + cell_recursively_shift_sparts(top, c, c->stars.parts, c->depth, &top->stars.parts[top->stars.count], 1); + } + + /* Recursively shift back the pointers of the direct parents */ + struct cell *up = c; + while (up->parent != NULL) { + up->stars.parts--; + up->stars.count++; + up = up->parent; } - message("depth=%d cellID=%d size=%d part=%p loc=[%e %e %e] w=[%e %e %e]", - top_new->depth, top_new->cellID, top_new->stars.count, top_new->stars.parts, - top_new->loc[0], top_new->loc[1], top_new->loc[2], - top_new->width[0], top_new->width[1], top_new->width[2]); + up->stars.parts--; + up->stars.count++; /* 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)); - message("sp->id=%lld sp->x=[%e %e %e]", sp->id, sp->x[0], sp->x[1], sp->x[2]); - message("c->loc=[%e %e %e], c->width=[%e %e %e]", c->loc[0], c->loc[1], c->loc[2], - c->width[0], c->width[1], c->width[2]); + //message("sp->id=%lld sp->x=[%e %e %e]", sp->id, sp->x[0], sp->x[1], sp->x[2]); + //message("c->loc=[%e %e %e], c->width=[%e %e %e]", c->loc[0], c->loc[1], c->loc[2], + //c->width[0], c->width[1], c->width[2]); /* 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]; - message("sp->id=%lld sp->x=[%e %e %e]", sp->id, sp->x[0], sp->x[1], sp->x[2]); + //message("sp->id=%lld sp->x=[%e %e %e]", sp->id, sp->x[0], sp->x[1], sp->x[2]); #ifdef SWIFT_DEBUG_CHECKS sp->ti_drift = e->ti_current; #endif #ifdef SWIFT_DEBUG_CHECKS - cell_check_spart_pos_and_time(top); + cell_check_spart_pos(top); #endif } diff --git a/src/runner.c b/src/runner.c index 4d891063bc..7a79146609 100644 --- a/src/runner.c +++ b/src/runner.c @@ -522,8 +522,8 @@ 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) { - message("Removing particle id=%lld rho=%e", p->id, rho); + if (rho > 1.7e7 && e->step > 2) { + message("c->cellID=%d Removing particle id=%lld rho=%e", c->cellID, p->id, rho); cell_convert_part_to_gpart(e, c, p, xp); cell_add_spart(e, c); } diff --git a/src/space.h b/src/space.h index c06ae26e1a..c422bce97b 100644 --- a/src/space.h +++ b/src/space.h @@ -46,7 +46,7 @@ struct cosmology; #define space_maxsize_default 8000000 #define space_extra_parts_default 0 #define space_extra_gparts_default 0 -#define space_extra_sparts_default 20 +#define space_extra_sparts_default 40 #define space_expected_max_nr_strays_default 100 #define space_subsize_pair_hydro_default 256000000 #define space_subsize_self_hydro_default 32000 -- GitLab From 14662c0cc4d93461f5c04d53a85ebd01d8327f82 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Fri, 16 Nov 2018 16:23:37 +0100 Subject: [PATCH 29/52] Version can now create multiple particles and put them in the correct cell array. Code can survive multiple rebuilds. --- src/cell.c | 111 +++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 77 insertions(+), 34 deletions(-) diff --git a/src/cell.c b/src/cell.c index a7c9f33cad..31bd4381f9 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3648,7 +3648,7 @@ void cell_check_timesteps(struct cell *c) { #endif } -void cell_check_spart_pos(const struct cell* c) { +void cell_check_spart_pos(const struct cell* c, const struct spart *global_sparts) { #ifdef SWIFT_DEBUG_CHECKS @@ -3656,7 +3656,7 @@ void cell_check_spart_pos(const struct cell* c) { if (c->split) { for(int k = 0; k < 8; ++k) if (c->progeny[k] != NULL) - cell_check_spart_pos(c->progeny[k]); + cell_check_spart_pos(c->progeny[k], global_sparts); } struct spart *sparts = c->stars.parts; @@ -3671,6 +3671,15 @@ void cell_check_spart_pos(const struct cell* c) { (sp->x[1] >= c->loc[1] + c->width[1]) || (sp->x[2] >= c->loc[2] + c->width[2])) 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 @@ -3679,18 +3688,34 @@ void cell_check_spart_pos(const struct cell* c) { } -void cell_recursively_shift_sparts(struct cell *c, struct cell *const target_cell, struct spart const* start, - const int depth, struct spart *temp, - const int direct_progeny){ - +/** + * @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) { - for(int k = 0; k < 8; ++k) { + + /* No need to recurse in progenies located before the insestion point */ + const int first_progeny = main_branch ? progeny_list[c->depth] : 0; + + for(int k = first_progeny; k < 8; ++k) { + if(c->progeny[k] != NULL) - cell_recursively_shift_sparts(c->progeny[k], target_cell, start, depth, temp, 0); + cell_recursively_shift_sparts(c->progeny[k], progeny_list, + main_branch && (k == first_progeny)); } } - if(c->stars.parts >= start) + /* 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++; } @@ -3698,12 +3723,29 @@ void cell_add_spart(const struct engine *e, struct cell *const c) { if (c->nodeID != engine_rank) error("Adding spart on a foreign node"); - /* Get the top-level this leaf cell is in */ + /* Progeny number at each level */ + int progeny[space_cell_maxdepth] = {0}; +#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; } + message("top->cellID=%d", top->cellID); + message("depth=%d progenies=[%d %d %d %d %d %d %d %d]", c->depth, + progeny[0], progeny[1], progeny[2], progeny[3], progeny[4], progeny[5], progeny[6], progeny[7]); + #ifdef SWIFT_DEBUG_CHECKS if (top->depth != 0) error("Cell-linking issue"); #endif @@ -3712,55 +3754,56 @@ void cell_add_spart(const struct engine *e, struct cell *const c) { if (top->stars.count == top->stars.count_total) error("We ran out of star particles!"); - /* Copy the particles to get a free space */ + /* Number of particles to shift in order to get a free space. */ const int n_copy = &top->stars.parts[top->stars.count] - c->stars.parts; + message("top->count=%d c->count=%d n_copy=%d", top->stars.count, c->stars.count, n_copy); + if(n_copy > 0) { + struct spart *temp = NULL; if (posix_memalign((void**)&temp, spart_align, n_copy * sizeof(struct spart)) != 0) error("Impossible to allocate temp buffer"); - + + /* Shift all the spart ahead of the current empty position by 1 */ memcpy(temp, c->stars.parts, n_copy * sizeof(struct spart)); memcpy(c->stars.parts + 1, temp, n_copy * sizeof(struct spart)); free(temp); - - message("count=%d n_copy=%d", c->stars.count, n_copy); - - /* Recursively shift all the stars to get a free spot at the start of the current cell*/ - cell_recursively_shift_sparts(top, c, c->stars.parts, c->depth, &top->stars.parts[top->stars.count], 1); + + /* Update the gpart->spart links (shift by 1) */ + for(int i = 0; i < n_copy; ++i) { + if(c->stars.parts[i+1].gpart != NULL) { + c->stars.parts[i+1].gpart->id_or_neg_offset--; + } else { +#ifdef SWIFT_DEBUG_CHECKS + // MATTHIEU + if(c->stars.parts[i+1].time_bin != time_bin_not_created) + error("Incorrectly linked spart!"); +#endif + } + } } + + /* 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); - /* Recursively shift back the pointers of the direct parents */ - struct cell *up = c; - while (up->parent != NULL) { - up->stars.parts--; - up->stars.count++; - up = up->parent; - } - up->stars.parts--; - up->stars.count++; - /* 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)); - //message("sp->id=%lld sp->x=[%e %e %e]", sp->id, sp->x[0], sp->x[1], sp->x[2]); - //message("c->loc=[%e %e %e], c->width=[%e %e %e]", c->loc[0], c->loc[1], c->loc[2], - //c->width[0], c->width[1], c->width[2]); + 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]; - - //message("sp->id=%lld sp->x=[%e %e %e]", sp->id, sp->x[0], sp->x[1], sp->x[2]); + sp->time_bin = time_bin_not_created; #ifdef SWIFT_DEBUG_CHECKS sp->ti_drift = e->ti_current; #endif #ifdef SWIFT_DEBUG_CHECKS - cell_check_spart_pos(top); + cell_check_spart_pos(top, e->s->sparts); #endif } -- GitLab From 6cbc97fbd52db9f4cbf56c7df792c9d3b1999b7d Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Fri, 16 Nov 2018 17:03:24 +0100 Subject: [PATCH 30/52] Allow the stars to move a bit in their cells in the debugging checks after star formation. --- src/cell.c | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/cell.c b/src/cell.c index 31bd4381f9..864d3f8e64 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3664,12 +3664,12 @@ void cell_check_spart_pos(const struct cell* c, const struct spart *global_spart for (int i = 0; i < count; ++i) { const struct spart *sp = &sparts[i]; - if( (sp->x[0] < c->loc[0]) || - (sp->x[1] < c->loc[1]) || - (sp->x[2] < c->loc[2]) || - (sp->x[0] >= c->loc[0] + c->width[0]) || - (sp->x[1] >= c->loc[1] + c->width[1]) || - (sp->x[2] >= c->loc[2] + c->width[2])) + 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) { @@ -3697,11 +3697,10 @@ void cell_check_spart_pos(const struct cell* c, const struct spart *global_spart */ 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[c->depth] : 0; + const int first_progeny = main_branch ? progeny_list[(int)c->depth] : 0; for(int k = first_progeny; k < 8; ++k) { -- GitLab From f8bf93e97346b3c7b260df49471f11663d459129 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Fri, 16 Nov 2018 17:03:36 +0100 Subject: [PATCH 31/52] Code formatting --- src/cell.c | 137 +++++++++++++++++++++++++++------------------------ src/cell.h | 8 +-- src/runner.c | 3 +- 3 files changed, 79 insertions(+), 69 deletions(-) diff --git a/src/cell.c b/src/cell.c index 864d3f8e64..6839fc765e 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3648,15 +3648,16 @@ void cell_check_timesteps(struct cell *c) { #endif } -void cell_check_spart_pos(const struct cell* c, const struct spart *global_sparts) { +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) + for (int k = 0; k < 8; ++k) if (c->progeny[k] != NULL) - cell_check_spart_pos(c->progeny[k], global_sparts); + cell_check_spart_pos(c->progeny[k], global_sparts); } struct spart *sparts = c->stars.parts; @@ -3664,55 +3665,61 @@ void cell_check_spart_pos(const struct cell* c, const struct spart *global_spart 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)) + 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) { - + 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!"); + 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. + * @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? + * @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) { +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)); + 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 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) + if (main_branch) c->stars.count++; else c->stars.parts++; @@ -3721,48 +3728,49 @@ void cell_recursively_shift_sparts(struct cell *c, const int progeny_list[space_ void cell_add_spart(const struct engine *e, struct cell *const c) { if (c->nodeID != engine_rank) error("Adding spart on a foreign node"); - + /* Progeny number at each level */ int progeny[space_cell_maxdepth] = {0}; #ifdef SWIFT_DEBUG_CHECKS - for(int i = 0; i < space_cell_maxdepth; ++i) - progeny[i] = -1; + 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; + for (int k = 0; k < 8; ++k) { + if (top->parent->progeny[k] == top) { + progeny[(int)top->parent->depth] = k; } } top = top->parent; } - + message("top->cellID=%d", top->cellID); - message("depth=%d progenies=[%d %d %d %d %d %d %d %d]", c->depth, - progeny[0], progeny[1], progeny[2], progeny[3], progeny[4], progeny[5], progeny[6], progeny[7]); - + message("depth=%d progenies=[%d %d %d %d %d %d %d %d]", c->depth, progeny[0], + progeny[1], progeny[2], progeny[3], progeny[4], progeny[5], + progeny[6], progeny[7]); + #ifdef SWIFT_DEBUG_CHECKS if (top->depth != 0) error("Cell-linking issue"); #endif - + /* Are there any extra particles left? */ if (top->stars.count == top->stars.count_total) error("We ran out of star particles!"); - + /* Number of particles to shift in order to get a free space. */ const int n_copy = &top->stars.parts[top->stars.count] - c->stars.parts; - message("top->count=%d c->count=%d n_copy=%d", top->stars.count, c->stars.count, n_copy); - - if(n_copy > 0) { - + message("top->count=%d c->count=%d n_copy=%d", top->stars.count, + c->stars.count, n_copy); + + if (n_copy > 0) { + struct spart *temp = NULL; - if (posix_memalign((void**)&temp, spart_align, - n_copy * sizeof(struct spart)) != 0) + if (posix_memalign((void **)&temp, spart_align, + n_copy * sizeof(struct spart)) != 0) error("Impossible to allocate temp buffer"); /* Shift all the spart ahead of the current empty position by 1 */ @@ -3771,32 +3779,33 @@ void cell_add_spart(const struct engine *e, struct cell *const c) { free(temp); /* Update the gpart->spart links (shift by 1) */ - for(int i = 0; i < n_copy; ++i) { - if(c->stars.parts[i+1].gpart != NULL) { - c->stars.parts[i+1].gpart->id_or_neg_offset--; + for (int i = 0; i < n_copy; ++i) { + if (c->stars.parts[i + 1].gpart != NULL) { + c->stars.parts[i + 1].gpart->id_or_neg_offset--; } else { #ifdef SWIFT_DEBUG_CHECKS - // MATTHIEU - if(c->stars.parts[i+1].time_bin != time_bin_not_created) - error("Incorrectly linked spart!"); + // MATTHIEU + if (c->stars.parts[i + 1].time_bin != time_bin_not_created) + error("Incorrectly linked spart!"); #endif } - } + } } - /* 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); - + /* 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]; sp->time_bin = time_bin_not_created; - + #ifdef SWIFT_DEBUG_CHECKS sp->ti_drift = e->ti_current; #endif @@ -3901,8 +3910,8 @@ void cell_remove_spart(const struct engine *e, struct cell *c, * * @return Pointer to the gpart the part has become. */ -struct gpart* 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) @@ -3943,8 +3952,8 @@ struct gpart* cell_convert_part_to_gpart(const struct engine *e, struct cell *c, * * @return Pointer to the gpart the spart has become. */ -struct gpart* 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) diff --git a/src/cell.h b/src/cell.h index a78aeb5bec..4177082ee0 100644 --- a/src/cell.h +++ b/src/cell.h @@ -672,10 +672,10 @@ 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); -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 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); 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); diff --git a/src/runner.c b/src/runner.c index 7a79146609..5f4738bb57 100644 --- a/src/runner.c +++ b/src/runner.c @@ -523,7 +523,8 @@ 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.7e7 && e->step > 2) { - message("c->cellID=%d Removing particle id=%lld rho=%e", c->cellID, p->id, rho); + message("c->cellID=%d Removing particle id=%lld rho=%e", c->cellID, + p->id, rho); cell_convert_part_to_gpart(e, c, p, xp); cell_add_spart(e, c); } -- GitLab From f71f7ba0221d0dea0da396729d4d505f45d4b110 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Fri, 16 Nov 2018 17:31:36 +0100 Subject: [PATCH 32/52] Correctly link the newly created star to the ddestroyed gpart. Can survive time integration for a few steps. Crashes at next rebuild still. --- src/cell.c | 31 ++++++++++--------------------- src/cell.h | 4 +++- src/runner.c | 45 +++++++++++++++++++++++++++++++++++++++++++-- 3 files changed, 56 insertions(+), 24 deletions(-) diff --git a/src/cell.c b/src/cell.c index 6839fc765e..eb9f0f5b01 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3725,12 +3725,12 @@ void cell_recursively_shift_sparts(struct cell *c, c->stars.parts++; } -void cell_add_spart(const struct engine *e, struct cell *const c) { +struct spart *cell_add_spart(const struct engine *e, struct cell *const c) { if (c->nodeID != engine_rank) error("Adding spart on a foreign node"); /* Progeny number at each level */ - int progeny[space_cell_maxdepth] = {0}; + int progeny[space_cell_maxdepth]; #ifdef SWIFT_DEBUG_CHECKS for (int i = 0; i < space_cell_maxdepth; ++i) progeny[i] = -1; #endif @@ -3747,15 +3747,6 @@ void cell_add_spart(const struct engine *e, struct cell *const c) { top = top->parent; } - message("top->cellID=%d", top->cellID); - message("depth=%d progenies=[%d %d %d %d %d %d %d %d]", c->depth, progeny[0], - progeny[1], progeny[2], progeny[3], progeny[4], progeny[5], - progeny[6], progeny[7]); - -#ifdef SWIFT_DEBUG_CHECKS - if (top->depth != 0) error("Cell-linking issue"); -#endif - /* Are there any extra particles left? */ if (top->stars.count == top->stars.count_total) error("We ran out of star particles!"); @@ -3763,11 +3754,11 @@ void cell_add_spart(const struct engine *e, struct cell *const c) { /* Number of particles to shift in order to get a free space. */ const int n_copy = &top->stars.parts[top->stars.count] - c->stars.parts; - message("top->count=%d c->count=%d n_copy=%d", top->stars.count, - c->stars.count, n_copy); - if (n_copy > 0) { + // MATTHIEU: This can be improved. We don't need to copy everything, just + // need to swap a few particles. + struct spart *temp = NULL; if (posix_memalign((void **)&temp, spart_align, n_copy * sizeof(struct spart)) != 0) @@ -3784,9 +3775,7 @@ void cell_add_spart(const struct engine *e, struct cell *const c) { c->stars.parts[i + 1].gpart->id_or_neg_offset--; } else { #ifdef SWIFT_DEBUG_CHECKS - // MATTHIEU - if (c->stars.parts[i + 1].time_bin != time_bin_not_created) - error("Incorrectly linked spart!"); + error("Incorrectly linked spart!"); #endif } } @@ -3804,15 +3793,15 @@ void cell_add_spart(const struct engine *e, struct cell *const c) { 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]; - sp->time_bin = time_bin_not_created; + + /* Set it to the current time-bin */ + sp->time_bin = e->min_active_bin; #ifdef SWIFT_DEBUG_CHECKS sp->ti_drift = e->ti_current; #endif -#ifdef SWIFT_DEBUG_CHECKS - cell_check_spart_pos(top, e->s->sparts); -#endif + return sp; } /** diff --git a/src/cell.h b/src/cell.h index 4177082ee0..9fd4158a6e 100644 --- a/src/cell.h +++ b/src/cell.h @@ -664,8 +664,10 @@ void cell_activate_drift_gpart(struct cell *c, struct scheduler *s); void cell_activate_sorts(struct cell *c, int sid, struct scheduler *s); void cell_clear_drift_flags(struct cell *c, void *data); void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data); +void cell_check_spart_pos(const struct cell *c, + const struct spart *global_sparts); int cell_has_tasks(struct cell *c); -void cell_add_spart(const struct engine *e, struct cell *c); +struct spart *cell_add_spart(const struct engine *e, struct cell *c); void cell_remove_part(const struct engine *e, struct cell *c, struct part *p, struct xpart *xp); void cell_remove_gpart(const struct engine *e, struct cell *c, diff --git a/src/runner.c b/src/runner.c index 5f4738bb57..04403c95a4 100644 --- a/src/runner.c +++ b/src/runner.c @@ -525,13 +525,54 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { if (rho > 1.7e7 && e->step > 2) { message("c->cellID=%d Removing particle id=%lld rho=%e", c->cellID, p->id, rho); - cell_convert_part_to_gpart(e, c, p, xp); - cell_add_spart(e, c); + + /* Destroy the gas particle and get it's gpart friend */ + struct gpart *gp = cell_convert_part_to_gpart(e, c, p, xp); + + /* Create a fresh (empty) spart */ + struct spart *sp = cell_add_spart(e, c); + + /* 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; +#endif + + /* Set a smoothing length */ + sp->h = max(c->stars.h_max, c->hydro.h_max); + + /* Set everything to a valid state */ + stars_init_spart(sp); } } } } + if (c->super == c) { +#ifdef SWIFT_DEBUG_CHECKS + /* Check that everything went OK */ + cell_check_spart_pos(c, e->s->sparts); +#endif + } + if (timer) TIMER_TOC(timer_do_star_formation); } -- GitLab From a72b92dde10d71147b0e6b4f035f9208f2503708 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Fri, 16 Nov 2018 19:03:42 +0100 Subject: [PATCH 33/52] Update the counter of extra particles when creating one. --- src/cell.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/cell.c b/src/cell.c index eb9f0f5b01..4e55906eb3 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3801,6 +3801,8 @@ struct spart *cell_add_spart(const struct engine *e, struct cell *const c) { sp->ti_drift = e->ti_current; #endif + e->s->nr_extra_sparts--; + return sp; } -- GitLab From 21f48d7e7c3a7b6447c053ab672657f76307ec00 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <matthieu.schaller@durham.ac.uk> Date: Sat, 17 Nov 2018 14:46:08 +0000 Subject: [PATCH 34/52] Correctly update the counters when adding a spart to the global array. Handle the case where we run out of free slots. --- src/cell.c | 13 +++++++++---- src/cell.h | 2 +- src/runner.c | 18 +++++++----------- 3 files changed, 17 insertions(+), 16 deletions(-) diff --git a/src/cell.c b/src/cell.c index 4e55906eb3..33fad78ccc 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3725,7 +3725,7 @@ void cell_recursively_shift_sparts(struct cell *c, c->stars.parts++; } -struct spart *cell_add_spart(const struct engine *e, struct cell *const c) { +struct spart *cell_add_spart(struct engine *e, struct cell *const c) { if (c->nodeID != engine_rank) error("Adding spart on a foreign node"); @@ -3748,8 +3748,11 @@ struct spart *cell_add_spart(const struct engine *e, struct cell *const c) { } /* Are there any extra particles left? */ - if (top->stars.count == top->stars.count_total) - error("We ran out of star particles!"); + if (top->stars.count == top->stars.count_total) { + 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 int n_copy = &top->stars.parts[top->stars.count] - c->stars.parts; @@ -3798,10 +3801,12 @@ struct spart *cell_add_spart(const struct engine *e, struct cell *const c) { sp->time_bin = e->min_active_bin; #ifdef SWIFT_DEBUG_CHECKS + /* Specify it was drifted to this point */ sp->ti_drift = e->ti_current; #endif - e->s->nr_extra_sparts--; + /* Register that we used one of the free slots. */ + atomic_dec(&e->s->nr_extra_sparts); return sp; } diff --git a/src/cell.h b/src/cell.h index 9fd4158a6e..11699a2ce5 100644 --- a/src/cell.h +++ b/src/cell.h @@ -667,7 +667,7 @@ 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); -struct spart *cell_add_spart(const struct engine *e, struct cell *c); +struct spart *cell_add_spart(struct engine *e, struct cell *c); void cell_remove_part(const struct engine *e, struct cell *c, struct part *p, struct xpart *xp); void cell_remove_gpart(const struct engine *e, struct cell *c, diff --git a/src/runner.c b/src/runner.c index 04403c95a4..6edb69d06e 100644 --- a/src/runner.c +++ b/src/runner.c @@ -492,7 +492,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; @@ -526,12 +526,15 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { message("c->cellID=%d Removing particle id=%lld rho=%e", c->cellID, p->id, rho); - /* Destroy the gas particle and get it's gpart friend */ - struct gpart *gp = cell_convert_part_to_gpart(e, c, p, xp); - /* Create a fresh (empty) spart */ struct spart *sp = cell_add_spart(e, c); + /* Did we run out of free spart slots? */ + if (sp == NULL) continue; + + /* 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; @@ -566,13 +569,6 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { } } - if (c->super == c) { -#ifdef SWIFT_DEBUG_CHECKS - /* Check that everything went OK */ - cell_check_spart_pos(c, e->s->sparts); -#endif - } - if (timer) TIMER_TOC(timer_do_star_formation); } -- GitLab From cf71e5b9c7ec69d6842891593cf250a9fa2f30e5 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Sat, 17 Nov 2018 16:45:29 +0100 Subject: [PATCH 35/52] Re-arrange the debugginc code to be more efficient. --- src/cell.c | 7 +++---- src/runner.c | 4 ++-- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/cell.c b/src/cell.c index 33fad78ccc..1b62d81b99 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3774,13 +3774,12 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) { /* Update the gpart->spart links (shift by 1) */ for (int i = 0; i < n_copy; ++i) { - if (c->stars.parts[i + 1].gpart != NULL) { - c->stars.parts[i + 1].gpart->id_or_neg_offset--; - } else { #ifdef SWIFT_DEBUG_CHECKS + if (c->stars.parts[i + 1].gpart == NULL) { error("Incorrectly linked spart!"); -#endif } +#endif + c->stars.parts[i + 1].gpart->id_or_neg_offset--; } } diff --git a/src/runner.c b/src/runner.c index 6edb69d06e..ab74d8e19c 100644 --- a/src/runner.c +++ b/src/runner.c @@ -529,8 +529,8 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { /* Create a fresh (empty) spart */ struct spart *sp = cell_add_spart(e, c); - /* Did we run out of free spart slots? */ - if (sp == NULL) continue; + /* Did we run out of free spart slots? */ + if (sp == NULL) continue; /* Destroy the gas particle and get it's gpart friend */ struct gpart *gp = cell_convert_part_to_gpart(e, c, p, xp); -- GitLab From 2d6c012394c07525fc1cbaa0adac72e0ba1a1f68 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Fri, 23 Nov 2018 10:11:56 +0100 Subject: [PATCH 36/52] Some improvements to the addition of stars on-the-fly. --- src/cell.c | 29 +++++++++-------------------- src/runner.c | 3 +-- 2 files changed, 10 insertions(+), 22 deletions(-) diff --git a/src/cell.c b/src/cell.c index 1b62d81b99..76a189cc5f 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3749,31 +3749,22 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) { /* Are there any extra particles left? */ if (top->stars.count == top->stars.count_total) { - message("We ran out of star particles!"); + //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 int n_copy = &top->stars.parts[top->stars.count] - c->stars.parts; + const size_t n_copy = &top->stars.parts[top->stars.count] - c->stars.parts; if (n_copy > 0) { // MATTHIEU: This can be improved. We don't need to copy everything, just // need to swap a few particles. - - struct spart *temp = NULL; - if (posix_memalign((void **)&temp, spart_align, - n_copy * sizeof(struct spart)) != 0) - error("Impossible to allocate temp buffer"); - - /* Shift all the spart ahead of the current empty position by 1 */ - memcpy(temp, c->stars.parts, n_copy * sizeof(struct spart)); - memcpy(c->stars.parts + 1, temp, n_copy * sizeof(struct spart)); - free(temp); + memmove(&c->stars.parts[1], &c->stars.parts[0], n_copy * sizeof(struct spart)); /* Update the gpart->spart links (shift by 1) */ - for (int i = 0; i < n_copy; ++i) { + 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!"); @@ -3805,7 +3796,8 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) { #endif /* Register that we used one of the free slots. */ - atomic_dec(&e->s->nr_extra_sparts); + const size_t one = 1; + atomic_sub(&e->s->nr_extra_sparts, one); return sp; } @@ -3990,7 +3982,6 @@ 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; - const int count_total = count_real + space_extra_parts; if (c->depth != 0 || c->nodeID != engine_rank) error("This function should only be called on local top-level cells!"); @@ -4008,7 +3999,7 @@ void cell_reorder_extra_parts(struct cell *c, const ptrdiff_t parts_offset) { } #ifdef SWIFT_DEBUG_CHECKS - if (first_not_extra >= count_total) + if (first_not_extra >= count_real + space_extra_parts) error("Looking for extra particles beyond this cell's range!"); #endif @@ -4033,7 +4024,6 @@ 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; - const int count_total = count_real + space_extra_sparts; if (c->depth != 0 || c->nodeID != engine_rank) error("This function should only be called on local top-level cells!"); @@ -4051,7 +4041,7 @@ void cell_reorder_extra_sparts(struct cell *c, const ptrdiff_t sparts_offset) { } #ifdef SWIFT_DEBUG_CHECKS - if (first_not_extra >= count_total) + if (first_not_extra >= count_real + space_extra_sparts) error("Looking for extra particles beyond this cell's range!"); #endif @@ -4076,7 +4066,6 @@ void cell_reorder_extra_gparts(struct cell *c, struct part *parts, struct gpart *gparts = c->grav.parts; const int count_real = c->grav.count; - const int count_total = count_real + space_extra_gparts; if (c->depth != 0 || c->nodeID != engine_rank) error("This function should only be called on local top-level cells!"); @@ -4094,7 +4083,7 @@ void cell_reorder_extra_gparts(struct cell *c, struct part *parts, } #ifdef SWIFT_DEBUG_CHECKS - if (first_not_extra >= count_total) + if (first_not_extra >= count_real + space_extra_gparts) error("Looking for extra particles beyond this cell's range!"); #endif diff --git a/src/runner.c b/src/runner.c index ab74d8e19c..dddeaad7fd 100644 --- a/src/runner.c +++ b/src/runner.c @@ -523,8 +523,7 @@ 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.7e7 && e->step > 2) { - message("c->cellID=%d Removing particle id=%lld rho=%e", c->cellID, - p->id, rho); + message("Removing particle id=%lld rho=%e", p->id, rho); /* Create a fresh (empty) spart */ struct spart *sp = cell_add_spart(e, c); -- GitLab From e80945ec54515b11d9d6e7f41906eced6263e552 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Fri, 23 Nov 2018 10:49:33 +0100 Subject: [PATCH 37/52] Corrected condition to create more spares at rebuild time. --- src/space.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/space.c b/src/space.c index dce5fff8fd..b924c4360e 100644 --- a/src/space.c +++ b/src/space.c @@ -862,7 +862,7 @@ void space_allocate_extras(struct space *s, int verbose) { /* Do we have enough space for the extra sparts (i.e. we haven't used up any) * ? */ - if (nr_sparts + expected_num_extra_sparts > size_sparts) { + if (nr_actual_sparts + expected_num_extra_sparts > nr_sparts) { /* Ok... need to put some more in the game */ -- GitLab From 7591ce2e1638d3947fd57de54310713e93e578e6 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Fri, 23 Nov 2018 12:10:23 +0100 Subject: [PATCH 38/52] Additional debugging checks when re-ordering sparts in top-level cells. --- src/cell.c | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/cell.c b/src/cell.c index 76a189cc5f..a74d44a47c 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3749,7 +3749,7 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) { /* Are there any extra particles left? */ if (top->stars.count == top->stars.count_total) { - //message("We ran out of star particles!"); + // message("We ran out of star particles!"); atomic_inc(&e->forcerebuild); return NULL; } @@ -3761,7 +3761,8 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) { // 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)); + 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) { @@ -4049,6 +4050,11 @@ void cell_reorder_extra_sparts(struct cell *c, const ptrdiff_t sparts_offset) { 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 } } } -- GitLab From 1a1184030ef0fcd7174fb536f7e4c81fef87d87c Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Fri, 23 Nov 2018 12:10:47 +0100 Subject: [PATCH 39/52] When creating a star from gas, also synchronise the drifting points. --- src/runner.c | 1 + 1 file changed, 1 insertion(+) diff --git a/src/runner.c b/src/runner.c index dddeaad7fd..ba13d06261 100644 --- a/src/runner.c +++ b/src/runner.c @@ -556,6 +556,7 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { #ifdef SWIFT_DEBUG_CHECKS sp->ti_kick = gp->ti_kick; + gp->ti_drift = sp->drift; #endif /* Set a smoothing length */ -- GitLab From daa89f8f01932490947e02357fa21beec9be1cbe Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <matthieu.schaller@durham.ac.uk> Date: Mon, 26 Nov 2018 12:33:15 +0000 Subject: [PATCH 40/52] Correctly activate the cell's hierarchy when creating a new spart. --- src/cell.c | 22 ++++++++++++++++++---- src/runner.c | 1 + 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/cell.c b/src/cell.c index 76a189cc5f..7274e9faba 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3681,7 +3681,7 @@ void cell_check_spart_pos(const struct cell *c, error("Unlinked spart!"); if (&global_sparts[-gp->id_or_neg_offset] != sp) - error("Incorrectly linked spart!"); + error("Incorrectly linked spart! sp->created=%d", sp->created); } } @@ -3727,7 +3727,9 @@ void cell_recursively_shift_sparts(struct cell *c, 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!"); /* Progeny number at each level */ int progeny[space_cell_maxdepth]; @@ -3748,8 +3750,8 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) { } /* Are there any extra particles left? */ - if (top->stars.count == top->stars.count_total) { - //message("We ran out of star particles!"); + if (top->stars.count == top->stars.count_total - 1) { + message("We ran out of star particles!"); atomic_inc(&e->forcerebuild); return NULL; } @@ -3757,6 +3759,11 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) { /* 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 @@ -3767,7 +3774,7 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) { 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!"); + message("Incorrectly linked spart!"); } #endif c->stars.parts[i + 1].gpart->id_or_neg_offset--; @@ -3790,6 +3797,13 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) { /* 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; diff --git a/src/runner.c b/src/runner.c index dddeaad7fd..2e73937e06 100644 --- a/src/runner.c +++ b/src/runner.c @@ -556,6 +556,7 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { #ifdef SWIFT_DEBUG_CHECKS sp->ti_kick = gp->ti_kick; + gp->ti_drift = sp->ti_drift; #endif /* Set a smoothing length */ -- GitLab From 1a20336d188a9c359b9eb25d4df548645881012a Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <matthieu.schaller@durham.ac.uk> Date: Mon, 26 Nov 2018 12:41:26 +0000 Subject: [PATCH 41/52] Restore error message when the moved spart are incorrectly linked. --- src/cell.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cell.c b/src/cell.c index 7274e9faba..47e7168849 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3774,7 +3774,7 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) { for (size_t i = 0; i < n_copy; ++i) { #ifdef SWIFT_DEBUG_CHECKS if (c->stars.parts[i + 1].gpart == NULL) { - message("Incorrectly linked spart!"); + error("Incorrectly linked spart!"); } #endif c->stars.parts[i + 1].gpart->id_or_neg_offset--; -- GitLab From d9ef48ce3187634ba868563b1e0530d43260bc10 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 26 Nov 2018 15:45:19 +0100 Subject: [PATCH 42/52] Post-merge fixes. --- src/cell.h | 34 +++++++++++++++++++++++++++++++++- src/runner.c | 2 +- src/space.c | 2 +- 3 files changed, 35 insertions(+), 3 deletions(-) diff --git a/src/cell.h b/src/cell.h index 73f4949835..d71df6fe4e 100644 --- a/src/cell.h +++ b/src/cell.h @@ -348,7 +348,7 @@ struct cell { /*! Do any of this cell's sub-cells need to be sorted? */ char do_sub_sort; - #ifdef SWIFT_DEBUG_CHECKS +#ifdef SWIFT_DEBUG_CHECKS /*! Last (integer) time the cell's sort arrays were updated. */ integertime_t ti_sort; @@ -469,6 +469,9 @@ struct cell { /*! Linked list of the tasks computing this cell's star density. */ struct link *density; + /*! The task computing this cell's sorts. */ + struct task *sorts; + /*! Max smoothing length in this cell. */ double h_max; @@ -490,6 +493,30 @@ struct cell { /*! Values of dx_max before the drifts, used for sub-cell tasks. */ float dx_max_part_old; + /*! Maximum particle movement in this cell since the last sort. */ + float dx_max_sort; + + /*! Values of dx_max_sort before the drifts, used for sub-cell tasks. */ + float dx_max_sort_old; + + /*! Bit mask of sort directions that will be needed in the next timestep. */ + unsigned int requires_sorts; + + /*! Pointer for the sorted indices. */ + struct entry *sort[13]; + + /*! Bit-mask indicating the sorted directions */ + unsigned int sorted; + + /*! Bit mask of sorts that need to be computed for this cell. */ + unsigned int do_sort; + + /*! Do any of this cell's sub-cells need to be sorted? */ + char do_sub_sort; + + /*! Maximum end of (integer) time step in this cell for gravity tasks. */ + integertime_t ti_end_min; + /*! Number of #spart updated in this cell. */ int updated; @@ -499,6 +526,11 @@ struct cell { /*! Is the #spart data of this cell being used in a sub-cell? */ int hold; +#ifdef SWIFT_DEBUG_CHECKS + /*! Last (integer) time the cell's sort arrays were updated. */ + integertime_t ti_sort; +#endif + } stars; #ifdef WITH_MPI diff --git a/src/runner.c b/src/runner.c index 130335c714..8b0ac61c2d 100644 --- a/src/runner.c +++ b/src/runner.c @@ -556,7 +556,7 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { #ifdef SWIFT_DEBUG_CHECKS sp->ti_kick = gp->ti_kick; - gp->ti_drift = sp->ti_drift; + gp->ti_drift = sp->ti_drift; #endif /* Set a smoothing length */ diff --git a/src/space.c b/src/space.c index 2776e54929..66cea7a7c5 100644 --- a/src/space.c +++ b/src/space.c @@ -252,7 +252,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements, #endif if (s->with_self_gravity) bzero(c->grav.multipole, sizeof(struct gravity_tensors)); - for (int i = 0; i < 13; i++) + for (int i = 0; i < 13; i++) { if (c->hydro.sort[i] != NULL) { free(c->hydro.sort[i]); c->hydro.sort[i] = NULL; -- GitLab From 92fe81e9389358fd7cd684c9fab7ed69566098f6 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 28 Nov 2018 12:19:29 +0100 Subject: [PATCH 43/52] Add runtime flag for star formation. Remove the now defunct sourceterm flag. --- README | 1 + examples/main.c | 22 ++++++++++++---------- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/README b/README index 747fe6da47..8462c9aead 100644 --- a/README +++ b/README @@ -34,6 +34,7 @@ Valid options are: -r Continue using restart files. -s Run with hydrodynamics. -S Run with stars. + -F Run with star formation. -b Run with stars feedback. -t {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified. -T Print timers every time-step. diff --git a/examples/main.c b/examples/main.c index 7340f244e3..6b77b9ffaa 100644 --- a/examples/main.c +++ b/examples/main.c @@ -98,6 +98,7 @@ void print_help_message(void) { printf(" %2s %14s %s\n", "-r", "", "Continue using restart files."); printf(" %2s %14s %s\n", "-s", "", "Run with hydrodynamics."); printf(" %2s %14s %s\n", "-S", "", "Run with stars."); + printf(" %2s %14s %s\n", "-F", "", "Run with star formation."); printf(" %2s %14s %s\n", "-b", "", "Run with stars feedback."); printf(" %2s %14s %s\n", "-t", "{int}", "The number of threads to use on each MPI rank. Defaults to 1 if not " @@ -189,11 +190,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; @@ -250,7 +251,7 @@ int main(int argc, char *argv[]) { } break; case 'F': - with_sourceterms = 1; + with_star_formation = 1; break; case 'g': with_external_gravity = 1; @@ -398,6 +399,14 @@ int main(int argc, char *argv[]) { if (myrank == 0) print_help_message(); return 1; } + if (!with_stars && with_star_formation) { + if (myrank == 0) + printf( + "Error: Cannot process star formation without stars, -S must be " + "chosen.\n"); + if (myrank == 0) print_help_message(); + return 1; + } /* Let's pin the main thread, now we know if affinity will be used. */ #if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE) @@ -916,10 +925,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; @@ -931,15 +936,12 @@ 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], -- GitLab From 868f5f2643fb39a6b98338233afa9401014742a2 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 28 Nov 2018 12:25:55 +0100 Subject: [PATCH 44/52] Completely remove the source term task. --- src/engine.c | 7 +- src/engine.h | 11 +- src/engine_maketasks.c | 7 - src/runner.c | 40 ---- src/sourceterms.c | 85 -------- src/sourceterms.h | 78 ------- src/sourceterms/sn_feedback/sn_feedback.h | 192 ------------------ .../sn_feedback/sn_feedback_struct.h | 45 ---- src/sourceterms_struct.h | 26 --- src/task.c | 2 - src/task.h | 1 - src/timers.c | 1 - src/timers.h | 1 - tools/task_plots/analyse_tasks.py | 1 - tools/task_plots/plot_tasks.py | 1 - 15 files changed, 6 insertions(+), 492 deletions(-) delete mode 100644 src/sourceterms.c delete mode 100644 src/sourceterms.h delete mode 100644 src/sourceterms/sn_feedback/sn_feedback.h delete mode 100644 src/sourceterms/sn_feedback/sn_feedback_struct.h delete mode 100644 src/sourceterms_struct.h diff --git a/src/engine.c b/src/engine.c index 472800b499..0827f10c63 100644 --- a/src/engine.c +++ b/src/engine.c @@ -110,7 +110,6 @@ const char *engine_policy_names[] = {"none", "drift everything", "reconstruct multi-poles", "cooling", - "sourceterms", "stars", "structure finding", "star formation", @@ -1923,9 +1922,6 @@ int engine_estimate_nr_tasks(struct engine *e) { if (e->policy & engine_policy_star_formation) { n1 += 1; } - if (e->policy & engine_policy_sourceterms) { - n1 += 2; - } if (e->policy & engine_policy_stars) { /* 1 self, 1 sort, 26/2 pairs */ n1 += 15; @@ -2579,8 +2575,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; } diff --git a/src/engine.h b/src/engine.h index eb73dc32d0..01c4aacb3a 100644 --- a/src/engine.h +++ b/src/engine.h @@ -71,13 +71,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]; /** diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c index 68841aa599..2c9fe3cb70 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/runner.c b/src/runner.c index 8b0ac61c2d..c6a31021fe 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" @@ -100,42 +99,6 @@ /* Import the stars loop functions. */ #include "runner_doiact_stars.h" -/** - * @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. @@ -2945,9 +2908,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/sourceterms.c b/src/sourceterms.c deleted file mode 100644 index 993045e615..0000000000 --- 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 407d2f1936..0000000000 --- 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 411673c37e..0000000000 --- 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 dd1842a671..0000000000 --- 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 b3c38986db..0000000000 --- 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/task.c b/src/task.c index a90482b269..f68d10f7cf 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; diff --git a/src/task.h b/src/task.h index 994b2b14c0..5bce55d6a2 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/timers.c b/src/timers.c index ccec0a9657..bd1a7e6274 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 48ca1e2763..d5d3529bd9 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 ca41970c68..f79a0090b0 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 1fe7bcbd11..1ff722a607 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", -- GitLab From 9f17e0cac97d2f4131d1fd54ca1fe6c457277f8e Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 28 Nov 2018 12:28:43 +0100 Subject: [PATCH 45/52] Remove the i/o debugging code for spare particles. --- src/gravity/Default/gravity_io.h | 13 +------------ src/hydro/Gadget2/hydro_io.h | 15 +-------------- src/stars/Default/stars_io.h | 16 ++-------------- 3 files changed, 4 insertions(+), 40 deletions(-) diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h index 0c3c643102..1ba3899e7e 100644 --- a/src/gravity/Default/gravity_io.h +++ b/src/gravity/Default/gravity_io.h @@ -92,15 +92,6 @@ INLINE static void darkmatter_read_particles(struct gpart* gparts, UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset); } -INLINE static void test_time_bin_gpart(const struct engine* e, - const struct gpart* gp, float* ret) { - - if (gp->time_bin >= time_bin_inhibited) - error("Writing inhibited or extra particle time_bin=%d", gp->time_bin); - - *ret = gp->time_bin; -} - /** * @brief Specifies which g-particle fields to write to a dataset * @@ -113,7 +104,7 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts, int* num_fields) { /* Say how much we want to write */ - *num_fields = 5; + *num_fields = 4; /* List what we want to write */ list[0] = io_make_output_field_convert_gpart( @@ -124,8 +115,6 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts, io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, gparts, mass); list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset); - list[4] = io_make_output_field_convert_gpart( - "TimeBin", FLOAT, 1, UNIT_CONV_NO_UNITS, gparts, test_time_bin_gpart); } #endif /* SWIFT_DEFAULT_GRAVITY_IO_H */ diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h index d776951767..ec7d34f7ad 100644 --- a/src/hydro/Gadget2/hydro_io.h +++ b/src/hydro/Gadget2/hydro_io.h @@ -128,16 +128,6 @@ INLINE static void convert_part_potential(const struct engine* e, ret[0] = 0.f; } -INLINE static void test_time_bin_part(const struct engine* e, - const struct part* p, - const struct xpart* xp, float* ret) { - - if (p->time_bin >= time_bin_inhibited) - error("Writing inhibited or extra particle time_bin=%d", p->time_bin); - - *ret = p->time_bin; -} - /** * @brief Specifies which particle fields to write to a dataset * @@ -151,7 +141,7 @@ INLINE static void hydro_write_particles(const struct part* parts, struct io_props* list, int* num_fields) { - *num_fields = 11; + *num_fields = 10; /* List what we want to write */ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3, @@ -178,9 +168,6 @@ INLINE static void hydro_write_particles(const struct part* parts, list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1, UNIT_CONV_POTENTIAL, parts, xparts, convert_part_potential); - list[10] = - io_make_output_field_convert_part("TimeBin", FLOAT, 1, UNIT_CONV_NO_UNITS, - parts, xparts, test_time_bin_part); #ifdef DEBUG_INTERACTIONS_SPH diff --git a/src/stars/Default/stars_io.h b/src/stars/Default/stars_io.h index 82d06d292f..a6c2768f71 100644 --- a/src/stars/Default/stars_io.h +++ b/src/stars/Default/stars_io.h @@ -22,15 +22,6 @@ #include "io_properties.h" #include "stars_part.h" -INLINE static void test_time_bin_spart(const struct engine *e, - const struct spart *sp, float *ret) { - - if (sp->time_bin >= time_bin_inhibited) - error("Writing inhibited or extra particle time_bin=%d", sp->time_bin); - - *ret = sp->time_bin; -} - /** * @brief Specifies which s-particle fields to read from a dataset * @@ -69,8 +60,8 @@ INLINE static void stars_write_particles(const struct spart *sparts, struct io_props *list, int *num_fields) { - /* Say how much we want to write */ - *num_fields = 6; + /* Say how much we want to read */ + *num_fields = 5; /* List what we want to read */ list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, @@ -83,9 +74,6 @@ INLINE static void stars_write_particles(const struct spart *sparts, sparts, id); list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH, sparts, h); - - list[5] = io_make_output_field_convert_spart( - "TimeBin", FLOAT, 1, UNIT_CONV_NO_UNITS, sparts, test_time_bin_spart); } /** -- GitLab From fbc6c104c5e44eb769b7c5e90afd27fb2e05d44e Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 28 Nov 2018 15:20:39 +0100 Subject: [PATCH 46/52] Also remove the sourceterm files from the Makefile. --- src/Makefile.am | 5 ++--- src/engine.c | 11 +---------- src/engine.h | 7 +------ src/swift.h | 1 - 4 files changed, 4 insertions(+), 20 deletions(-) diff --git a/src/Makefile.am b/src/Makefile.am index 9b0610667b..78846115d0 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 \ @@ -57,7 +57,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 \ @@ -75,7 +75,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/engine.c b/src/engine.c index 0827f10c63..59fb2faf98 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" @@ -4015,8 +4014,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params, struct pm_mesh *mesh, const struct external_potential *potential, const 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)); @@ -4079,7 +4077,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 @@ -5175,7 +5172,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); @@ -5272,11 +5268,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 01c4aacb3a..ffcdb14e74 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" @@ -359,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; @@ -416,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, const 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/swift.h b/src/swift.h index 153c4ae0d4..7279fc58b3 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" -- GitLab From 8dfc2f12ac46b1a49b59f28e504c6c8e98a4a1f6 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 28 Nov 2018 15:33:11 +0100 Subject: [PATCH 47/52] Only activate the creation of extra spart when we are running the star formation policy. --- examples/main.c | 5 ++--- src/space.c | 18 ++++++++++++------ src/space.h | 7 +++++-- 3 files changed, 19 insertions(+), 11 deletions(-) diff --git a/examples/main.c b/examples/main.c index 6b77b9ffaa..596344d327 100644 --- a/examples/main.c +++ b/examples/main.c @@ -140,7 +140,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; @@ -844,7 +843,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); @@ -947,7 +946,7 @@ int main(int argc, char *argv[]) { 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/src/space.c b/src/space.c index 66cea7a7c5..7cfdaa9a71 100644 --- a/src/space.c +++ b/src/space.c @@ -977,7 +977,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { space_regrid(s, verbose); /* Allocate extra space for particles that will be created */ - space_allocate_extras(s, verbose); + 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; @@ -1027,9 +1027,9 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { /* Initialise the counters, including buffer space for future particles */ for (int i = 0; i < s->nr_cells; ++i) { - cell_part_counts[i] = 0; // space_extra_parts; - cell_gpart_counts[i] = 0; // space_extra_gparts; - cell_spart_counts[i] = 0; // space_extra_sparts; + cell_part_counts[i] = 0; + cell_gpart_counts[i] = 0; + cell_spart_counts[i] = 0; } /* Run through the particles and get their cell index. */ @@ -1567,7 +1567,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) { /* Re-order the extra particles such that they are at the end of their cell's memory pool. */ - space_reorder_extras(s, verbose); + 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. */ @@ -3630,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. * @@ -3643,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)); @@ -3655,6 +3657,7 @@ void space_init(struct space *s, struct swift_params *params, s->periodic = periodic; s->with_self_gravity = self_gravity; s->with_hydro = hydro; + s->with_star_formation = star_formation; s->nr_parts = Npart; s->nr_gparts = Ngpart; s->nr_sparts = Nspart; @@ -3857,6 +3860,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); } diff --git a/src/space.h b/src/space.h index c422bce97b..a1280945d2 100644 --- a/src/space.h +++ b/src/space.h @@ -46,7 +46,7 @@ struct cosmology; #define space_maxsize_default 8000000 #define space_extra_parts_default 0 #define space_extra_gparts_default 0 -#define space_extra_sparts_default 40 +#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 @@ -96,6 +96,9 @@ struct space { /*! Are we doing gravity? */ int with_self_gravity; + /*! Are we doing star formation? */ + int with_star_formation; + /*! Width of the top-level cells. */ double width[3]; @@ -266,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); -- GitLab From 9cf92b55ae3f29edbe8055af11f8c0fb75cceb99 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 28 Nov 2018 15:46:36 +0100 Subject: [PATCH 48/52] Read the number of extra part/gpart/spart from the yaml file. --- examples/parameter_example.yml | 3 +++ src/space.c | 6 ++++++ 2 files changed, 9 insertions(+) diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 9ae5fef665..78e46e9fff 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/space.c b/src/space.c index 7cfdaa9a71..fa34146d80 100644 --- a/src/space.c +++ b/src/space.c @@ -3754,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, -- GitLab From b44437b1cf2323cda1752d3aa7c40e30ddc38e7a Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 28 Nov 2018 15:51:27 +0100 Subject: [PATCH 49/52] Prevent users from running with star formation when running the MPI version --- examples/main.c | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/examples/main.c b/examples/main.c index 596344d327..719a6a459d 100644 --- a/examples/main.c +++ b/examples/main.c @@ -513,10 +513,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) -- GitLab From 865819781b9cc5a7453b3b1e1a43ca7f08356993 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 28 Nov 2018 17:43:17 +0100 Subject: [PATCH 50/52] Turned the full part --> spart conversion into a single function. --- src/cell.c | 95 ++++++++++++++++++++++++++++++++++++++++++++++++++-- src/cell.h | 4 ++- src/runner.c | 36 ++------------------ 3 files changed, 98 insertions(+), 37 deletions(-) diff --git a/src/cell.c b/src/cell.c index 91a1eef7b0..98d51b09fc 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3893,11 +3893,26 @@ void cell_recursively_shift_sparts(struct cell *c, c->stars.parts++; } +/** + * @breif "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]; @@ -4071,6 +4086,9 @@ 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. @@ -4078,7 +4096,8 @@ void cell_remove_spart(const struct engine *e, struct cell *c, * @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. + * @return Pointer to the #gpart the #part has become. It carries the + * ID of the #part and has a dark matter type. */ struct gpart *cell_convert_part_to_gpart(const struct engine *e, struct cell *c, struct part *p, struct xpart *xp) { @@ -4114,13 +4133,17 @@ struct gpart *cell_convert_part_to_gpart(const struct engine *e, struct cell *c, * @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. + * @return Pointer to the #gpart the #spart has become. It carries the + * ID of the #spart and has a dark matter type. */ struct gpart *cell_convert_spart_to_gpart(const struct engine *e, struct cell *c, struct spart *sp) { @@ -4152,6 +4175,74 @@ struct gpart *cell_convert_spart_to_gpart(const struct engine *e, return gp; } +/** + * @param "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. diff --git a/src/cell.h b/src/cell.h index d71df6fe4e..7178698c42 100644 --- a/src/cell.h +++ b/src/cell.h @@ -707,17 +707,19 @@ 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); -struct spart *cell_add_spart(struct engine *e, struct cell *c); void cell_remove_part(const struct engine *e, struct cell *c, struct part *p, struct xpart *xp); void cell_remove_gpart(const struct engine *e, struct cell *c, struct gpart *gp); void cell_remove_spart(const struct engine *e, struct cell *c, struct spart *sp); +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); diff --git a/src/runner.c b/src/runner.c index c6a31021fe..dd0e627ed1 100644 --- a/src/runner.c +++ b/src/runner.c @@ -488,43 +488,11 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) { if (rho > 1.7e7 && e->step > 2) { message("Removing particle id=%lld rho=%e", p->id, rho); - /* Create a fresh (empty) spart */ - struct spart *sp = cell_add_spart(e, c); + struct spart *sp = cell_convert_part_to_spart(e, c, p, xp); - /* Did we run out of free spart slots? */ + /* Did we run out of fresh particles? */ if (sp == NULL) continue; - /* 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); - /* Set everything to a valid state */ stars_init_spart(sp); } -- GitLab From c74e0ea52f33e7ed9b9d7414e84490df26093345 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Thu, 29 Nov 2018 10:25:04 +0100 Subject: [PATCH 51/52] Doxygen and RTD documentation fixes --- README | 2 +- doc/RTD/source/CommandLineOptions/index.rst | 4 ++++ examples/main.c | 2 +- src/cell.c | 4 ++-- src/engine.c | 1 - src/runner.c | 3 --- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/README b/README index 8462c9aead..bccc45e543 100644 --- a/README +++ b/README @@ -31,7 +31,7 @@ Valid options are: -n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop. -o {str} Generate a default output parameter file. -P {sec:par:val} Set parameter value and overwrites values read from the parameters file. Can be used more than once. - -r Continue using restart files. + -r Resume a run using restart files. -s Run with hydrodynamics. -S Run with stars. -F Run with star formation. diff --git a/doc/RTD/source/CommandLineOptions/index.rst b/doc/RTD/source/CommandLineOptions/index.rst index a8ce7af082..e23eccd5bf 100644 --- a/doc/RTD/source/CommandLineOptions/index.rst +++ b/doc/RTD/source/CommandLineOptions/index.rst @@ -31,13 +31,17 @@ can be found by typing ``./swift -h``. + ``-o``: {str} Generate a default output parameter file. + ``-P``: {sec:par:val} Set parameter value and overwrites values read from the parameters file. Can be used more than once. ++ ``-r``: Resume a run using restart files. + ``-s``: Run with hydrodynamics. + ``-S``: Run with stars. ++ ``-F``: Run with star formation. ++ ``-b``: Run with stars feedback. + ``-t``: {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified. + ``-T``: Print timers every time-step. + ``-v``: [12] Increase the level of verbosity: 1, MPI-rank 0 writes, 2, All MPI-ranks write. ++ ``-x``: Run with structure finding. + ``-y``: {int} Time-step frequency at which task graphs are dumped. + ``-Y``: {int} Time-step frequency at which threadpool tasks are dumped. + ``-h``: Print a help message and exit. diff --git a/examples/main.c b/examples/main.c index 719a6a459d..4e0cc86f1d 100644 --- a/examples/main.c +++ b/examples/main.c @@ -95,7 +95,7 @@ void print_help_message(void) { printf(" %2s %14s %s\n", "-P", "{sec:par:val}", "Set parameter value and overwrites values read from the parameters " "file. Can be used more than once."); - printf(" %2s %14s %s\n", "-r", "", "Continue using restart files."); + printf(" %2s %14s %s\n", "-r", "", "Resume a run using restart files."); printf(" %2s %14s %s\n", "-s", "", "Run with hydrodynamics."); printf(" %2s %14s %s\n", "-S", "", "Run with stars."); printf(" %2s %14s %s\n", "-F", "", "Run with star formation."); diff --git a/src/cell.c b/src/cell.c index 98d51b09fc..ed54163a00 100644 --- a/src/cell.c +++ b/src/cell.c @@ -3894,7 +3894,7 @@ void cell_recursively_shift_sparts(struct cell *c, } /** - * @breif "Add" a #spart in a given #cell. + * @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 @@ -4176,7 +4176,7 @@ struct gpart *cell_convert_spart_to_gpart(const struct engine *e, } /** - * @param "Remove" a #part from a #cell and replace it with a #spart + * @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 diff --git a/src/engine.c b/src/engine.c index 59fb2faf98..09010fdc60 100644 --- a/src/engine.c +++ b/src/engine.c @@ -4002,7 +4002,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, diff --git a/src/runner.c b/src/runner.c index dd0e627ed1..dacefa498c 100644 --- a/src/runner.c +++ b/src/runner.c @@ -586,9 +586,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) { \ -- GitLab From 187062ab1f72ad4807089c1c3b932223713ace60 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Sat, 1 Dec 2018 15:16:55 +0100 Subject: [PATCH 52/52] Change the -F flag from sourceterms to star formation. --- README.md | 2 +- examples/main.c | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 94e95776cd..92925017d9 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/examples/main.c b/examples/main.c index 02289bfae7..63e8511b26 100644 --- a/examples/main.c +++ b/examples/main.c @@ -185,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, -- GitLab