Commit e6632efb by Matthieu Schaller

### Sanitize top-level cells at startups to avoid situations with ridiculously large h values.

parent 6fae7e8e
 ... ... @@ -679,6 +679,55 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset) { part_relink_parts(gparts, gcount, parts - parts_offset); } /** * @brief Sanitizes the smoothing length values of cells by setting large * outliers to more sensible values. * * We compute the mean and standard deviation of the smoothing lengths in * logarithmic space and limit values to mean + 4 sigma. * * @param c The cell. */ void cell_sanitize(struct cell *c) { const int count = c->count; struct part *parts = c->parts; /* First collect some statistics */ float h_mean = 0.f, h_mean2 = 0.f; float h_min = FLT_MAX, h_max = 0.f; for (int i = 0; i < count; ++i) { const float h = log10f(parts[i].h); h_mean += h; h_mean2 += h * h; h_max = max(h_max, h); h_min = min(h_min, h); } h_mean /= count; h_mean2 /= count; const float h_var = h_mean2 - h_mean * h_mean; const float h_std = sqrtf(h_var); /* Choose a cut */ const float h_limit = pow(10.f, h_mean + 4.f * h_std); /* Be verbose this is not innocuous */ message("Cell properties: h_min= %f h_max= %f geometric mean= %f", powf(10.f, h_min), powf(10.f, h_max), powf(10.f, h_mean)); if (c->h_max > h_limit) { message("Smoothing lengths will be limited to (mean + 4sigma)= %f", h_limit); /* Apply the cut */ for (int i = 0; i < count; ++i) parts->h = min(parts[i].h, h_limit); c->h_max = h_limit; } } /** * @brief Converts hydro quantities to a valid state after the initial density * calculation ... ...
 ... ... @@ -218,6 +218,7 @@ struct cell { /* Function prototypes. */ void cell_split(struct cell *c, ptrdiff_t parts_offset); void cell_sanitize(struct cell *c); int cell_locktree(struct cell *c); void cell_unlocktree(struct cell *c); int cell_glocktree(struct cell *c); ... ...
 ... ... @@ -2243,6 +2243,8 @@ void engine_rebuild(struct engine *e) { /* Re-build the space. */ space_rebuild(e->s, 0.0, e->verbose); if (e->ti_current == 0) space_sanitize(e->s); /* If in parallel, exchange the cell structure. */ #ifdef WITH_MPI engine_exchange_cells(e); ... ...
 ... ... @@ -58,6 +58,7 @@ int space_splitsize = space_splitsize_default; int space_subsize = space_subsize_default; int space_maxsize = space_maxsize_default; int space_maxcount = space_maxcount_default; /* Map shift vector to sortlist. */ const int sortlistID[27] = { ... ... @@ -698,7 +699,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) { /* At this point, we have the upper-level cells, old or new. Now make sure that the parts in each cell are ok. */ space_split(s, cells_top, verbose); space_split(s, cells_top, s->nr_cells, verbose); if (verbose) message("took %.3f %s.", clocks_from_ticks(getticks() - tic), ... ... @@ -711,14 +712,16 @@ void space_rebuild(struct space *s, double cell_max, int verbose) { * This is done in parallel using threads in the #threadpool. * * @param s The #space. * @param cells The cell hierarchy * @param cells The cell hierarchy. * @param nr_cells The number of cells. * @param verbose Are we talkative ? */ void space_split(struct space *s, struct cell *cells, int verbose) { void space_split(struct space *s, struct cell *cells, int nr_cells, int verbose) { const ticks tic = getticks(); threadpool_map(&s->e->threadpool, space_split_mapper, cells, s->nr_cells, threadpool_map(&s->e->threadpool, space_split_mapper, cells, nr_cells, sizeof(struct cell), 1, s); if (verbose) ... ... @@ -726,6 +729,29 @@ void space_split(struct space *s, struct cell *cells, int verbose) { clocks_getunit()); } /** * @brief Runs through the top-level cells and checks whether tasks associated * with them can be split. If not, try to sanitize the cells. * * @param s The #space to act upon. */ void space_sanitize(struct space *s) { for (int k = 0; k < s->nr_cells; k++) { struct cell *c = &s->cells_top[k]; const double min_width = c->dmin; /* Do we have a problem ? */ if (c->h_max * kernel_gamma * space_stretch > min_width * 0.5 && c->count > space_maxcount) { /* Ok, clean-up the mess */ cell_sanitize(c); } } } /** * @brief Sort the particles and condensed particles according to the given * indices. ... ... @@ -1247,17 +1273,17 @@ void space_map_cells_pre(struct space *s, int full, * too many particles. * * @param map_data Pointer towards the top-cells. * @param num_elements The number of cells to treat. * @param num_cells The number of cells to treat. * @param extra_data Pointers to the #space. */ void space_split_mapper(void *map_data, int num_elements, void *extra_data) { void space_split_mapper(void *map_data, int num_cells, void *extra_data) { /* Unpack the inputs. */ struct space *s = (struct space *)extra_data; struct cell *restrict cells_top = (struct cell *)map_data; struct engine *e = s->e; for (int ind = 0; ind < num_elements; ind++) { for (int ind = 0; ind < num_cells; ind++) { struct cell *c = &cells_top[ind]; ... ... @@ -1578,6 +1604,8 @@ void space_init(struct space *s, const struct swift_params *params, space_subsize_default); space_splitsize = parser_get_opt_param_int( params, "Scheduler:cell_split_size", space_splitsize_default); space_maxcount = parser_get_opt_param_int(params, "Scheduler:cell_max_count", space_maxcount_default); if (verbose) message("max_size set to %d, sub_size set to %d, split_size set to %d", space_maxsize, space_subsize, space_splitsize); ... ...
 ... ... @@ -41,6 +41,7 @@ #define space_splitsize_default 400 #define space_maxsize_default 8000000 #define space_subsize_default 64000000 #define space_maxcount_default 100000 #define space_stretch 1.10f #define space_maxreldx 0.25f ... ... @@ -51,6 +52,7 @@ extern int space_splitsize; extern int space_maxsize; extern int space_subsize; extern int space_maxcount; /* Map shift vector to sortlist. */ extern const int sortlistID[27]; ... ... @@ -148,6 +150,7 @@ void space_init(struct space *s, const struct swift_params *params, double dim[3], struct part *parts, struct gpart *gparts, size_t Npart, size_t Ngpart, int periodic, int gravity, 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); void space_map_parts(struct space *s, ... ... @@ -164,7 +167,8 @@ void space_gparts_sort_mapper(void *map_data, int num_elements, void *extra_data); void space_rebuild(struct space *s, double h_max, int verbose); void space_recycle(struct space *s, struct cell *c); void space_split(struct space *s, struct cell *cells, int verbose); void space_split(struct space *s, struct cell *cells, int nr_cells, int verbose); void space_split_mapper(void *map_data, int num_elements, void *extra_data); void space_do_parts_sort(); void space_do_gparts_sort(); ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!