From cf636e161d52b73282dbaecad472e402aa6abbe6 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Wed, 6 Apr 2022 16:25:39 +0200 Subject: [PATCH] Applied code formatting script --- examples/main.c | 24 +- src/cell.h | 60 +- src/cell_split.c | 4 +- src/cell_unskip.c | 6 +- src/common_io.h | 25 +- src/common_io_cells.c | 43 +- src/distributed_io.c | 74 +- src/engine_maketasks.c | 75 +- src/engine_marktasks.c | 6 +- src/engine_proxy.c | 2 +- src/engine_redistribute.c | 36 +- src/mesh_gravity.c | 48 +- src/mesh_gravity.h | 8 +- src/parallel_io.c | 36 +- src/partition.c | 31 +- src/runner_doiact_grav.c | 3 +- src/runner_others.c | 15 +- src/runner_time_integration.c | 13 +- src/serial_io.c | 36 +- src/single_io.c | 38 +- src/space.c | 45 +- src/space.h | 25 +- src/space_cell_index.c | 62 +- src/space_extras.c | 27 +- src/space_rebuild.c | 61 +- src/space_regrid.c | 424 +++++----- src/space_split.c | 30 +- src/task.c | 3 +- src/zoom_partition.c | 169 ++-- src/zoom_region.c | 1431 +++++++++++++++++---------------- src/zoom_region.h | 31 +- src/zoom_regrid.c | 477 +++++------ tests/testFFT.c | 4 +- 33 files changed, 1807 insertions(+), 1565 deletions(-) diff --git a/examples/main.c b/examples/main.c index 337af20999..bc4ad63c5d 100644 --- a/examples/main.c +++ b/examples/main.c @@ -1312,22 +1312,22 @@ int main(int argc, char *argv[]) { neutrino_props_init(&neutrino_properties, &prog_const, &us, params, &cosmo, with_neutrinos); - /* Initialise the gravity properties */ - bzero(&gravity_properties, sizeof(struct gravity_props)); - if (with_self_gravity) - gravity_props_init( - &gravity_properties, params, &prog_const, &cosmo, with_cosmology, - with_external_gravity, with_baryon_particles, with_DM_particles, - with_neutrinos, with_DM_background_particles, periodic, dim); + /* Initialise the gravity properties */ + bzero(&gravity_properties, sizeof(struct gravity_props)); + if (with_self_gravity) + gravity_props_init( + &gravity_properties, params, &prog_const, &cosmo, with_cosmology, + with_external_gravity, with_baryon_particles, with_DM_particles, + with_neutrinos, with_DM_background_particles, periodic, dim); /* Initialize the space with these data. */ if (myrank == 0) clocks_gettime(&tic); space_init(&s, params, &cosmo, dim, &hydro_properties, &gravity_properties, - parts, gparts, sinks, sparts, bparts, Ngas, Ngpart, Nsink, Nspart, - Nbpart, Nnupart, periodic, replicate, remap_ids, generate_gas_in_ics, - with_hydro, with_self_gravity, with_star_formation, with_sink, - with_DM_background_particles, with_neutrinos, talking, dry_run, - nr_nodes); + parts, gparts, sinks, sparts, bparts, Ngas, Ngpart, Nsink, + Nspart, Nbpart, Nnupart, periodic, replicate, remap_ids, + generate_gas_in_ics, with_hydro, with_self_gravity, + with_star_formation, with_sink, with_DM_background_particles, + with_neutrinos, talking, dry_run, nr_nodes); /* Initialise the line of sight properties. */ if (with_line_of_sight) los_init(s.dim, &los_properties, params); diff --git a/src/cell.h b/src/cell.h index 06ce875cf3..2443acdc35 100644 --- a/src/cell.h +++ b/src/cell.h @@ -326,9 +326,9 @@ enum cell_flags { * @brief What kind of top level cell is this cell? * * 0 = A background top level cell. - * 1 = A background top level cell that is within the gravity criterion of the zoom region. - * 2 = An ignored background top level cell (as it contains the zoom region). - * 3 = A top level zoom cell. + * 1 = A background top level cell that is within the gravity criterion of the + * zoom region. 2 = An ignored background top level cell (as it contains the + * zoom region). 3 = A top level zoom cell. */ enum tl_cell_types { tl_cell, tl_cell_neighbour, void_tl_cell, zoom_tl_cell }; #endif @@ -666,7 +666,8 @@ int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj, * @param x, y, z Coordinates of particle/cell. */ __attribute__((always_inline)) INLINE int cell_getid_pos(const struct space *s, - const double x, const double y, + const double x, + const double y, const double z) { /* Define variable to output */ int cell_id; @@ -684,7 +685,6 @@ __attribute__((always_inline)) INLINE int cell_getid_pos(const struct space *s, const int j = y * s->iwidth[1]; const int k = z * s->iwidth[2]; cell_id = cell_getid(s->cdim, i, j, k); - } #else /* Not compiled with zoom regions so we can use the simple version */ @@ -696,7 +696,7 @@ __attribute__((always_inline)) INLINE int cell_getid_pos(const struct space *s, return cell_id; } - /** +/** * @brief Does a #cell contain no particle at all. * * @param c The #cell. @@ -722,12 +722,24 @@ __attribute__((always_inline)) INLINE static double cell_min_dist2_same_size( const int periodic, const double dim[3]) { #ifdef SWIFT_DEBUG_CHECKS - if (ci->width[0] != cj->width[0]) error("x cells of different size! (ci->width=[%f %f %f] cj->width=[%f %f %f])", - ci->width[0], ci->width[1], ci->width[2], cj->width[0], cj->width[1], cj->width[2]); - if (ci->width[1] != cj->width[1]) error("y cells of different size! (ci->width=[%f %f %f] cj->width=[%f %f %f])", - ci->width[0], ci->width[1], ci->width[2], cj->width[0], cj->width[1], cj->width[2]); - if (ci->width[2] != cj->width[2]) error("z cells of different size! (ci->width=[%f %f %f] cj->width=[%f %f %f])", - ci->width[0], ci->width[1], ci->width[2], cj->width[0], cj->width[1], cj->width[2]); + if (ci->width[0] != cj->width[0]) + error( + "x cells of different size! (ci->width=[%f %f %f] cj->width=[%f %f " + "%f])", + ci->width[0], ci->width[1], ci->width[2], cj->width[0], cj->width[1], + cj->width[2]); + if (ci->width[1] != cj->width[1]) + error( + "y cells of different size! (ci->width=[%f %f %f] cj->width=[%f %f " + "%f])", + ci->width[0], ci->width[1], ci->width[2], cj->width[0], cj->width[1], + cj->width[2]); + if (ci->width[2] != cj->width[2]) + error( + "z cells of different size! (ci->width=[%f %f %f] cj->width=[%f %f " + "%f])", + ci->width[0], ci->width[1], ci->width[2], cj->width[0], cj->width[1], + cj->width[2]); #endif const double cix_min = ci->loc[0]; @@ -1403,36 +1415,40 @@ __attribute__((always_inline)) INLINE void cell_assign_top_level_cell_index( #ifdef WITH_ZOOM_REGION - /* Get some information from the space */ - const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; + /* Get some information from the space */ + const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; - /* Get some zoom information from the space */ - const int zoom_cdim[3] = {s->zoom_props->cdim[0], s->zoom_props->cdim[1], s->zoom_props->cdim[2]}; + /* Get some zoom information from the space */ + const int zoom_cdim[3] = {s->zoom_props->cdim[0], s->zoom_props->cdim[1], + s->zoom_props->cdim[2]}; - if (((cdim[0] * cdim[1] * cdim[2]) + (zoom_cdim[0] * zoom_cdim[1] * zoom_cdim[2])) > 32 * 32 * 32) { - /* print warning only once */ + if (((cdim[0] * cdim[1] * cdim[2]) + + (zoom_cdim[0] * zoom_cdim[1] * zoom_cdim[2])) > 32 * 32 * 32) { + /* print warning only once */ if (last_cell_id == 1ULL) { message( "WARNING: Got (%d x %d x %d + %d x %d x %d) top level cells. " "Cell IDs are only guaranteed to be " "reproduceably unique if count is < 32^3", - cdim[0], cdim[1], cdim[2], zoom_cdim[0], zoom_cdim[1], zoom_cdim[2]); + cdim[0], cdim[1], cdim[2], zoom_cdim[0], zoom_cdim[1], + zoom_cdim[2]); } /* Do this in same line. Otherwise, bad things happen. */ c->cellID = atomic_inc(&last_cell_id); } else { - c->cellID = (unsigned long long)(cell_getid_pos(s, c->loc[0], c->loc[1], c->loc[2])); + c->cellID = (unsigned long long)(cell_getid_pos(s, c->loc[0], c->loc[1], + c->loc[2])); } /* in both cases, keep track of first prodigy index */ atomic_inc(&last_leaf_cell_id); #else /* Get some information from the space */ - const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; + const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; const double iwidth[3] = {s->iwidth[0], s->iwidth[1], s->iwidth[2]}; if (cdim[0] * cdim[1] * cdim[2] > 32 * 32 * 32) { - /* print warning only once */ + /* print warning only once */ if (last_cell_id == 1ULL) { message( "WARNING: Got > %d x %d x %d top level cells. " diff --git a/src/cell_split.c b/src/cell_split.c index 74b486b0b1..e942a1e1a0 100644 --- a/src/cell_split.c +++ b/src/cell_split.c @@ -76,12 +76,12 @@ void cell_split(struct cell *c, const ptrdiff_t parts_offset, int bucket_offset[9]; #ifdef SWIFT_DEBUG_CHECKS - /* Check that the buffs are OK and particles are within the cells prior to split. */ + /* Check that the buffs are OK and particles are within the cells prior to + * split. */ for (int k = 0; k < count; k++) { if (buff[k].x[0] != parts[k].x[0] || buff[k].x[1] != parts[k].x[1] || buff[k].x[2] != parts[k].x[2]) error("Inconsistent buff contents."); - } for (int k = 0; k < gcount; k++) { if (gbuff[k].x[0] != gparts[k].x[0] || gbuff[k].x[1] != gparts[k].x[1] || diff --git a/src/cell_unskip.c b/src/cell_unskip.c index f47e1c728f..38d6100256 100644 --- a/src/cell_unskip.c +++ b/src/cell_unskip.c @@ -1524,8 +1524,10 @@ int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s) { #ifdef SWIFT_DEBUG_CHECKS /* Ensure we are not rebuilding on a zoom and natural cell pair */ if (ci->tl_cell_type <= 2 || cj->tl_cell_type <= 2) - error("We just decided to rebuild based on a hydro zoom and natural cell pair. " - "This should never happen!"); + error( + "We just decided to rebuild based on a hydro zoom and natural " + "cell pair. " + "This should never happen!"); #endif } diff --git a/src/common_io.h b/src/common_io.h index 58c41cd5a3..bd6ef258f7 100644 --- a/src/common_io.h +++ b/src/common_io.h @@ -107,19 +107,18 @@ void io_write_code_description(hid_t h_file); void io_write_engine_policy(hid_t h_file, const struct engine* e); void io_write_part_type_names(hid_t h_grp); -void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const int zoom_cdim[3], const double dim[3], - const struct cell* cells_top, const int nr_cells, - const int nr_zoomcells, const int nr_bkgcells, - const double width[3], const double zoom_width[3], const int nodeID, - const int distributed, - const int subsample[swift_type_count], - const float subsample_fraction[swift_type_count], - const int snap_num, - const long long global_counts[swift_type_count], - const long long global_offsets[swift_type_count], - const int num_fields[swift_type_count], - const struct unit_system* internal_units, - const struct unit_system* snapshot_units, const int with_zoom); +void io_write_cell_offsets( + hid_t h_grp, const int cdim[3], const int zoom_cdim[3], const double dim[3], + const struct cell* cells_top, const int nr_cells, const int nr_zoomcells, + const int nr_bkgcells, const double width[3], const double zoom_width[3], + const int nodeID, const int distributed, + const int subsample[swift_type_count], + const float subsample_fraction[swift_type_count], const int snap_num, + const long long global_counts[swift_type_count], + const long long global_offsets[swift_type_count], + const int num_fields[swift_type_count], + const struct unit_system* internal_units, + const struct unit_system* snapshot_units, const int with_zoom); void io_read_unit_system(hid_t h_file, struct unit_system* ic_units, const struct unit_system* internal_units, diff --git a/src/common_io_cells.c b/src/common_io_cells.c index e48a2a5ef1..ee689da146 100644 --- a/src/common_io_cells.c +++ b/src/common_io_cells.c @@ -449,14 +449,18 @@ void io_write_array(hid_t h_grp, const int n, const void* array, * * @param h_grp the hdf5 group to write to. * @param cdim The number of top-level cells along each axis. - * @param zoom_cdim The number of top-level zoom cells along each axis. (only used with zoom region) + * @param zoom_cdim The number of top-level zoom cells along each axis. (only + * used with zoom region) * @param dim The box size. * @param cells_top The top-level cells. * @param nr_cells The number of top-level cells. - * @param nr_zoomcells The number of top-level zoom cells. (only used with zoom region) - * @param nr_bkgcells The number of top-level background cells. (only used with zoom region) + * @param nr_zoomcells The number of top-level zoom cells. (only used with zoom + * region) + * @param nr_bkgcells The number of top-level background cells. (only used with + * zoom region) * @param width The width of a top-level cell. - * @param zoom_width The width of a top-level zoom cell. (only used with zoom region) + * @param zoom_width The width of a top-level zoom cell. (only used with zoom + * region) * @param nodeID This node's MPI rank. * @param distributed Is this a distributed snapshot? * @param subsample Are we subsampling the different particle types? @@ -468,21 +472,21 @@ void io_write_array(hid_t h_grp, const int n, const void* array, * @param numFields The number of fields to write for each particle type. * @param internal_units The internal unit system. * @param snapshot_units The snapshot unit system. - * @param with_zoom Flag for whether we're running with a zoom region. (only used with zoom region) + * @param with_zoom Flag for whether we're running with a zoom region. (only + * used with zoom region) */ -void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const int zoom_cdim[3], const double dim[3], - const struct cell* cells_top, const int nr_cells, - const int nr_zoomcells, const int nr_bkgcells, - const double width[3], const double zoom_width[3], const int nodeID, - const int distributed, - const int subsample[swift_type_count], - const float subsample_fraction[swift_type_count], - const int snap_num, - const long long global_counts[swift_type_count], - const long long global_offsets[swift_type_count], - const int num_fields[swift_type_count], - const struct unit_system* internal_units, - const struct unit_system* snapshot_units, const int with_zoom) { +void io_write_cell_offsets( + hid_t h_grp, const int cdim[3], const int zoom_cdim[3], const double dim[3], + const struct cell* cells_top, const int nr_cells, const int nr_zoomcells, + const int nr_bkgcells, const double width[3], const double zoom_width[3], + const int nodeID, const int distributed, + const int subsample[swift_type_count], + const float subsample_fraction[swift_type_count], const int snap_num, + const long long global_counts[swift_type_count], + const long long global_offsets[swift_type_count], + const int num_fields[swift_type_count], + const struct unit_system* internal_units, + const struct unit_system* snapshot_units, const int with_zoom) { #ifdef SWIFT_DEBUG_CHECKS if (distributed) { @@ -503,8 +507,7 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const int zoom_cdim[3 #ifdef WITH_ZOOM_REGION if (with_zoom) - for (int ijk = 0; ijk < 3; ijk++) - zoom_cell_width[ijk] = zoom_width[ijk]; + for (int ijk = 0; ijk < 3; ijk++) zoom_cell_width[ijk] = zoom_width[ijk]; #endif /* Temporary memory for the cell-by-cell information */ diff --git a/src/distributed_io.c b/src/distributed_io.c index ea1b5f41e3..866560b4bb 100644 --- a/src/distributed_io.c +++ b/src/distributed_io.c @@ -1074,27 +1074,33 @@ void write_output_distributed(struct engine* e, h_grp = H5Gcreate(h_file, "/Cells", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); if (h_grp < 0) error("Error while creating cells group"); - /* Write the location of the particles in the arrays */ + /* Write the location of the particles in the arrays */ #ifdef WITH_ZOOM_REGION if (e->s->with_zoom_region) { - io_write_cell_offsets(h_grp, e->s->cdim, e->s->zoom_props->cdim, e->s->dim, e->s->cells_top, - e->s->nr_cells, e->s->zoom_props->nr_zoom_cells, e->s->zoom_props->nr_bkg_cells, e->s->width, - e->s->zoom_props->width, mpi_rank, /*distributed=*/1, subsample, subsample_fraction, - e->snapshot_output_count, N_total, global_offsets, - numFields, internal_units, snapshot_units, /*with_zoom=*/1); + io_write_cell_offsets( + h_grp, e->s->cdim, e->s->zoom_props->cdim, e->s->dim, e->s->cells_top, + e->s->nr_cells, e->s->zoom_props->nr_zoom_cells, + e->s->zoom_props->nr_bkg_cells, e->s->width, e->s->zoom_props->width, + mpi_rank, /*distributed=*/1, subsample, subsample_fraction, + e->snapshot_output_count, N_total, global_offsets, numFields, + internal_units, snapshot_units, /*with_zoom=*/1); } else { - io_write_cell_offsets(h_grp, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, - e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, e->s->width, - /*zoom_width=*/NULL, mpi_rank, /*distributed=*/1, subsample, subsample_fraction, - e->snapshot_output_count, N_total, global_offsets, - numFields, internal_units, snapshot_units, /*with_zoom=*/0); + io_write_cell_offsets( + h_grp, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, + e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, + e->s->width, + /*zoom_width=*/NULL, mpi_rank, /*distributed=*/1, subsample, + subsample_fraction, e->snapshot_output_count, N_total, global_offsets, + numFields, internal_units, snapshot_units, /*with_zoom=*/0); } #else - io_write_cell_offsets(h_grp, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, - e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, e->s->width, - /*zoom_width=*/NULL, mpi_rank, /*distributed=*/1, subsample, subsample_fraction, - e->snapshot_output_count, N_total, global_offsets, - numFields, internal_units, snapshot_units, /*with_zoom=*/0); + io_write_cell_offsets(h_grp, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, + e->s->cells_top, e->s->nr_cells, /*nr_zoomcells=*/0, + /*nr_bkgcells=*/e->s->nr_cells, e->s->width, + /*zoom_width=*/NULL, mpi_rank, /*distributed=*/1, + subsample, subsample_fraction, e->snapshot_output_count, + N_total, global_offsets, numFields, internal_units, + snapshot_units, /*with_zoom=*/0); #endif H5Gclose(h_grp); @@ -1484,24 +1490,30 @@ void write_output_distributed(struct engine* e, /* Write the location of the particles in the arrays */ #ifdef WITH_ZOOM_REGION if (e->s->with_zoom_region) { - io_write_cell_offsets(h_grp_cells, e->s->cdim, e->s->zoom_props->cdim, e->s->dim, e->s->cells_top, - e->s->nr_cells, e->s->zoom_props->nr_zoom_cells, e->s->zoom_props->nr_bkg_cells, e->s->width, - e->s->zoom_props->width, mpi_rank, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, global_offsets, - numFields, internal_units, snapshot_units, /*with_zoom=*/1); + io_write_cell_offsets( + h_grp_cells, e->s->cdim, e->s->zoom_props->cdim, e->s->dim, + e->s->cells_top, e->s->nr_cells, e->s->zoom_props->nr_zoom_cells, + e->s->zoom_props->nr_bkg_cells, e->s->width, e->s->zoom_props->width, + mpi_rank, /*distributed=*/0, subsample, subsample_fraction, + e->snapshot_output_count, N_total, global_offsets, numFields, + internal_units, snapshot_units, /*with_zoom=*/1); } else { - io_write_cell_offsets(h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, - e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, e->s->width, - /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, global_offsets, - numFields, internal_units, snapshot_units, /*with_zoom=*/0); + io_write_cell_offsets( + h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, + e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, + e->s->width, + /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, subsample, + subsample_fraction, e->snapshot_output_count, N_total, global_offsets, + numFields, internal_units, snapshot_units, /*with_zoom=*/0); } #else - io_write_cell_offsets(h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, - e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, e->s->width, - /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, global_offsets, - numFields, internal_units, snapshot_units, /*with_zoom=*/0); + io_write_cell_offsets(h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, + e->s->cells_top, e->s->nr_cells, /*nr_zoomcells=*/0, + /*nr_bkgcells=*/e->s->nr_cells, e->s->width, + /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, + subsample, subsample_fraction, e->snapshot_output_count, + N_total, global_offsets, numFields, internal_units, + snapshot_units, /*with_zoom=*/0); #endif /* Close everything */ diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c index f861bef659..8296c9cded 100644 --- a/src/engine_maketasks.c +++ b/src/engine_maketasks.c @@ -1280,10 +1280,11 @@ void engine_make_hierarchical_tasks_gravity(struct engine *e, struct cell *c) { /* Gravity non-neighbouring pm calculations */ if (c->top->tl_cell_type == 3) { c->grav.long_range = scheduler_addtask( - s, task_type_grav_long_range, task_subtype_none, 0, 0, c, NULL); + s, task_type_grav_long_range, task_subtype_none, 0, 0, c, NULL); } else { - c->grav.long_range = scheduler_addtask( - s, task_type_grav_long_range_bkg, task_subtype_none, 0, 0, c, NULL); + c->grav.long_range = + scheduler_addtask(s, task_type_grav_long_range_bkg, + task_subtype_none, 0, 0, c, NULL); } /* Gravity recursive down-pass */ @@ -1861,7 +1862,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements, ci, cj); #ifdef SWIFT_DEBUG_CHECKS - #ifdef WITH_MPI +#ifdef WITH_MPI /* Let's cross-check that we had a proxy for that cell */ if (ci->nodeID == nodeID && cj->nodeID != engine_rank) { @@ -1907,14 +1908,13 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements, } #endif /* WITH_MPI */ #endif /* SWIFT_DEBUG_CHECKS */ - } - } - } - } - } + } + } + } + } + } } - /** * @brief Constructs the top-level tasks for the external gravity. * @@ -4259,16 +4259,17 @@ void engine_make_fof_tasks(struct engine *e) { /* Construct a FOF loop over neighbours */ if (e->policy & engine_policy_fof) { #ifdef WITH_ZOOM_REGION - if (s->with_zoom_region) { - threadpool_map(&e->threadpool, engine_make_fofloop_tasks_mapper_with_zoom, NULL, - s->zoom_props->tl_cell_offset, 1, threadpool_auto_chunk_size, e); - } else { - threadpool_map(&e->threadpool, engine_make_fofloop_tasks_mapper, NULL, - s->nr_cells, 1, threadpool_auto_chunk_size, e); - } + if (s->with_zoom_region) { + threadpool_map(&e->threadpool, engine_make_fofloop_tasks_mapper_with_zoom, + NULL, s->zoom_props->tl_cell_offset, 1, + threadpool_auto_chunk_size, e); + } else { + threadpool_map(&e->threadpool, engine_make_fofloop_tasks_mapper, NULL, + s->nr_cells, 1, threadpool_auto_chunk_size, e); + } #else - threadpool_map(&e->threadpool, engine_make_fofloop_tasks_mapper, NULL, - s->nr_cells, 1, threadpool_auto_chunk_size, e); + threadpool_map(&e->threadpool, engine_make_fofloop_tasks_mapper, NULL, + s->nr_cells, 1, threadpool_auto_chunk_size, e); #endif } @@ -4327,16 +4328,17 @@ void engine_maketasks(struct engine *e) { /* Construct the first hydro loop over neighbours */ if (e->policy & engine_policy_hydro) { #ifdef WITH_ZOOM_REGION - if (s->with_zoom_region) { - threadpool_map(&e->threadpool, engine_make_hydroloop_tasks_mapper_with_zoom, NULL, - s->zoom_props->nr_zoom_cells, 1, threadpool_auto_chunk_size, e); + if (s->with_zoom_region) { + threadpool_map( + &e->threadpool, engine_make_hydroloop_tasks_mapper_with_zoom, NULL, + s->zoom_props->nr_zoom_cells, 1, threadpool_auto_chunk_size, e); } else { - threadpool_map(&e->threadpool, engine_make_hydroloop_tasks_mapper, NULL, - s->nr_cells, 1, threadpool_auto_chunk_size, e); + threadpool_map(&e->threadpool, engine_make_hydroloop_tasks_mapper, NULL, + s->nr_cells, 1, threadpool_auto_chunk_size, e); } #else - threadpool_map(&e->threadpool, engine_make_hydroloop_tasks_mapper, NULL, - s->nr_cells, 1, threadpool_auto_chunk_size, e); + threadpool_map(&e->threadpool, engine_make_hydroloop_tasks_mapper, NULL, + s->nr_cells, 1, threadpool_auto_chunk_size, e); #endif } @@ -4350,19 +4352,22 @@ void engine_maketasks(struct engine *e) { if (e->policy & engine_policy_self_gravity) { #ifdef WITH_ZOOM_REGION if (s->with_zoom_region) { - threadpool_map(&e->threadpool, engine_make_self_gravity_tasks_mapper_zoom_cells, NULL, - s->zoom_props->nr_zoom_cells, 1, threadpool_auto_chunk_size, e); - threadpool_map(&e->threadpool, engine_make_self_gravity_tasks_mapper_natural_cells, NULL, - s->zoom_props->nr_bkg_cells, 1, threadpool_auto_chunk_size, e); - threadpool_map(&e->threadpool, engine_make_self_gravity_tasks_mapper_with_zoom_diffsize, NULL, - s->nr_cells, 1, threadpool_auto_chunk_size, e); + threadpool_map( + &e->threadpool, engine_make_self_gravity_tasks_mapper_zoom_cells, + NULL, s->zoom_props->nr_zoom_cells, 1, threadpool_auto_chunk_size, e); + threadpool_map( + &e->threadpool, engine_make_self_gravity_tasks_mapper_natural_cells, + NULL, s->zoom_props->nr_bkg_cells, 1, threadpool_auto_chunk_size, e); + threadpool_map(&e->threadpool, + engine_make_self_gravity_tasks_mapper_with_zoom_diffsize, + NULL, s->nr_cells, 1, threadpool_auto_chunk_size, e); } else { - threadpool_map(&e->threadpool, engine_make_self_gravity_tasks_mapper, NULL, - s->nr_cells, 1, threadpool_auto_chunk_size, e); + threadpool_map(&e->threadpool, engine_make_self_gravity_tasks_mapper, + NULL, s->nr_cells, 1, threadpool_auto_chunk_size, e); } #else threadpool_map(&e->threadpool, engine_make_self_gravity_tasks_mapper, NULL, - s->nr_cells, 1, threadpool_auto_chunk_size, e); + s->nr_cells, 1, threadpool_auto_chunk_size, e); #endif } diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c index 531580bdc3..d169eea782 100644 --- a/src/engine_marktasks.c +++ b/src/engine_marktasks.c @@ -805,8 +805,10 @@ void engine_marktasks_mapper(void *map_data, int num_elements, #ifdef SWIFT_DEBUG_CHECKS /* Ensure we are not rebuilding on a zoom and natural cell pair */ - if (ci->tl_cell_type <= 2 || cj->tl_cell_type <= 2) - error("We just decided to rebuild based on a hydro zoom and natural cell pair. " + if (ci->tl_cell_type <= 2 || cj->tl_cell_type <= 2) + error( + "We just decided to rebuild based on a hydro zoom and natural " + "cell pair. " "This should never happen!"); #endif } diff --git a/src/engine_proxy.c b/src/engine_proxy.c index 3787747611..434632940e 100644 --- a/src/engine_proxy.c +++ b/src/engine_proxy.c @@ -42,7 +42,7 @@ void engine_makeproxies(struct engine *e) { #ifdef WITH_MPI #ifdef WITH_ZOOM_REGION - engine_makeproxies_with_zoom_region(e); + engine_makeproxies_with_zoom_region(e); #else /* Let's time this */ diff --git a/src/engine_redistribute.c b/src/engine_redistribute.c index cfc81ae590..5993696c1f 100644 --- a/src/engine_redistribute.c +++ b/src/engine_redistribute.c @@ -275,8 +275,8 @@ struct redist_mapper_data { else if (parts[k].x[j] >= s->dim[j]) \ parts[k].x[j] -= s->dim[j]; \ } \ - const int cid = cell_getid_pos(s, parts[k].x[0], \ - parts[k].x[1], parts[k].x[2]); \ + const int cid = \ + cell_getid_pos(s, parts[k].x[0], parts[k].x[1], parts[k].x[2]); \ dest[k] = s->cells_top[cid].nodeID; \ size_t ind = mydata->nodeID * mydata->nr_nodes + dest[k]; \ lcounts[ind] += 1; \ @@ -754,7 +754,7 @@ void engine_redistribute(struct engine *e) { if (sp->time_bin == time_bin_not_created) error("Inhibited particle found after sorting!"); - + /* New cell index */ const int new_cid = cell_getid_pos(s, sp->x[0], sp->x[1], sp->x[2]); @@ -818,7 +818,7 @@ void engine_redistribute(struct engine *e) { if (bp->time_bin == time_bin_not_created) error("Inhibited particle found after sorting!"); - + /* New cell index */ const int new_cid = cell_getid_pos(s, bp->x[0], bp->x[1], bp->x[2]); @@ -891,13 +891,19 @@ void engine_redistribute(struct engine *e) { const int new_node = c->nodeID; if (g_dest[k] != new_node) - error("gpart's new node index not matching sorted index (%d != %d). cid=%d pos=%f %f %f cwidth=%f", - g_dest[k], new_node, new_cid, gp->x[0], gp->x[1], gp->x[2], c->width[0]); + error( + "gpart's new node index not matching sorted index (%d != %d). cid=%d " + "pos=%f %f %f cwidth=%f", + g_dest[k], new_node, new_cid, gp->x[0], gp->x[1], gp->x[2], + c->width[0]); if (gp->x[0] < c->loc[0] || gp->x[0] > c->loc[0] + c->width[0] || gp->x[1] < c->loc[1] || gp->x[1] > c->loc[1] + c->width[1] || gp->x[2] < c->loc[2] || gp->x[2] > c->loc[2] + c->width[2]) - error("gpart not sorted into the right top-level cell! (cellid: %d, cell_type: %d)", new_cid, c->tl_cell_type); + error( + "gpart not sorted into the right top-level cell! (cellid: %d, " + "cell_type: %d)", + new_cid, c->tl_cell_type); } #endif @@ -1147,8 +1153,8 @@ void engine_redistribute(struct engine *e) { /* Verify that all parts are in the right place. */ for (size_t k = 0; k < nr_parts_new; k++) { - const int cid = cell_getid_pos(s, s->parts[k].x[0], - s->parts[k].x[1], s->parts[k].x[2]); + const int cid = + cell_getid_pos(s, s->parts[k].x[0], s->parts[k].x[1], s->parts[k].x[2]); if (cells[cid].nodeID != nodeID) error("Received particle (%zu) that does not belong here (nodeID=%i).", k, @@ -1156,8 +1162,8 @@ void engine_redistribute(struct engine *e) { } for (size_t k = 0; k < nr_gparts_new; k++) { - const int cid = cell_getid_pos(s, s->gparts[k].x[0], - s->gparts[k].x[1], s->gparts[k].x[2]); + const int cid = cell_getid_pos(s, s->gparts[k].x[0], s->gparts[k].x[1], + s->gparts[k].x[2]); if (cells[cid].nodeID != nodeID) error("Received g-particle (%zu) that does not belong here (nodeID=%i).", @@ -1165,8 +1171,8 @@ void engine_redistribute(struct engine *e) { } for (size_t k = 0; k < nr_sparts_new; k++) { - const int cid = cell_getid_pos(s, s->sparts[k].x[0], - s->sparts[k].x[1], s->sparts[k].x[2]); + const int cid = cell_getid_pos(s, s->sparts[k].x[0], s->sparts[k].x[1], + s->sparts[k].x[2]); if (cells[cid].nodeID != nodeID) error("Received s-particle (%zu) that does not belong here (nodeID=%i).", @@ -1174,8 +1180,8 @@ void engine_redistribute(struct engine *e) { } for (size_t k = 0; k < nr_bparts_new; k++) { - const int cid = cell_getid_pos(s, s->bparts[k].x[0], - s->bparts[k].x[1], s->bparts[k].x[2]); + const int cid = cell_getid_pos(s, s->bparts[k].x[0], s->bparts[k].x[1], + s->bparts[k].x[2]); if (cells[cid].nodeID != nodeID) error("Received b-particle (%zu) that does not belong here (nodeID=%i).", diff --git a/src/mesh_gravity.c b/src/mesh_gravity.c index 8a202b5488..5ad2916df2 100644 --- a/src/mesh_gravity.c +++ b/src/mesh_gravity.c @@ -281,7 +281,8 @@ void cell_gpart_to_mesh_CIC_mapper(void* map_data, int num, void* extra) { #ifdef WITH_ZOOM_REGION /** - * @brief A wrapper for the threadpool mapper function for the mesh CIC assignment of a background cell. + * @brief A wrapper for the threadpool mapper function for the mesh CIC + * assignment of a background cell. * * @param map_data A chunk of the list of local cells. * @param num The number of cells in the chunk. @@ -292,7 +293,8 @@ void bkg_cell_gpart_to_mesh_CIC_mapper(void* map_data, int num, void* extra) { } /** - * @brief A wrapper for the threadpool mapper function for the mesh CIC assignment of a zoom cell. + * @brief A wrapper for the threadpool mapper function for the mesh CIC + * assignment of a zoom cell. * * @param map_data A chunk of the list of local cells. * @param num The number of cells in the chunk. @@ -508,7 +510,8 @@ void cell_mesh_to_gpart_CIC_mapper(void* map_data, int num, void* extra) { #ifdef WITH_ZOOM_REGION /** - * @brief A wrapper for the threadpool mapper function for the mesh CIC assignment of a background cell. + * @brief A wrapper for the threadpool mapper function for the mesh CIC + * assignment of a background cell. * * @param map_data A chunk of the list of local cells. * @param num The number of cells in the chunk. @@ -519,7 +522,8 @@ void bkg_cell_mesh_to_gpart_CIC_mapper(void* map_data, int num, void* extra) { } /** - * @brief A wrapper for the threadpool mapper function for the mesh CIC assignment of a zoom cell. + * @brief A wrapper for the threadpool mapper function for the mesh CIC + * assignment of a zoom cell. * * @param map_data A chunk of the list of local cells. * @param num The number of cells in the chunk. @@ -955,16 +959,18 @@ void compute_potential_global(struct pm_mesh* mesh, const struct space* s, * the local top-level cells */ #ifdef WITH_ZOOM_REGION if (s->with_zoom_region) { - threadpool_map(tp, bkg_cell_gpart_to_mesh_CIC_mapper, (void*)s->zoom_props->local_bkg_cells_with_particles_top, - s->zoom_props->nr_local_bkg_cells_with_particles, sizeof(int), threadpool_auto_chunk_size, - (void*)&data); - threadpool_map(tp, zoom_cell_gpart_to_mesh_CIC_mapper, (void*)s->zoom_props->local_zoom_cells_with_particles_top, - s->zoom_props->nr_local_zoom_cells_with_particles, sizeof(int), threadpool_auto_chunk_size, - (void*)&data); + threadpool_map(tp, bkg_cell_gpart_to_mesh_CIC_mapper, + (void*)s->zoom_props->local_bkg_cells_with_particles_top, + s->zoom_props->nr_local_bkg_cells_with_particles, + sizeof(int), threadpool_auto_chunk_size, (void*)&data); + threadpool_map(tp, zoom_cell_gpart_to_mesh_CIC_mapper, + (void*)s->zoom_props->local_zoom_cells_with_particles_top, + s->zoom_props->nr_local_zoom_cells_with_particles, + sizeof(int), threadpool_auto_chunk_size, (void*)&data); } else { threadpool_map(tp, cell_gpart_to_mesh_CIC_mapper, (void*)local_cells, - nr_local_cells, sizeof(int), threadpool_auto_chunk_size, - (void*)&data); + nr_local_cells, sizeof(int), threadpool_auto_chunk_size, + (void*)&data); } #else threadpool_map(tp, cell_gpart_to_mesh_CIC_mapper, (void*)local_cells, @@ -1073,16 +1079,18 @@ void compute_potential_global(struct pm_mesh* mesh, const struct space* s, the local top-level cells */ #ifdef WITH_ZOOM_REGION if (s->with_zoom_region) { - threadpool_map(tp, bkg_cell_mesh_to_gpart_CIC_mapper, (void*)s->zoom_props->local_bkg_cells_with_particles_top, - s->zoom_props->nr_local_bkg_cells_with_particles, sizeof(int), threadpool_auto_chunk_size, - (void*)&data); - threadpool_map(tp, zoom_cell_mesh_to_gpart_CIC_mapper, (void*)s->zoom_props->local_zoom_cells_with_particles_top, - s->zoom_props->nr_local_zoom_cells_with_particles, sizeof(int), threadpool_auto_chunk_size, - (void*)&data); + threadpool_map(tp, bkg_cell_mesh_to_gpart_CIC_mapper, + (void*)s->zoom_props->local_bkg_cells_with_particles_top, + s->zoom_props->nr_local_bkg_cells_with_particles, + sizeof(int), threadpool_auto_chunk_size, (void*)&data); + threadpool_map(tp, zoom_cell_mesh_to_gpart_CIC_mapper, + (void*)s->zoom_props->local_zoom_cells_with_particles_top, + s->zoom_props->nr_local_zoom_cells_with_particles, + sizeof(int), threadpool_auto_chunk_size, (void*)&data); } else { threadpool_map(tp, cell_mesh_to_gpart_CIC_mapper, (void*)local_cells, - nr_local_cells, sizeof(int), threadpool_auto_chunk_size, - (void*)&data); + nr_local_cells, sizeof(int), threadpool_auto_chunk_size, + (void*)&data); } #else threadpool_map(tp, cell_mesh_to_gpart_CIC_mapper, (void*)local_cells, diff --git a/src/mesh_gravity.h b/src/mesh_gravity.h index b9a1ddeaeb..01bb448a06 100644 --- a/src/mesh_gravity.h +++ b/src/mesh_gravity.h @@ -97,9 +97,9 @@ void pm_mesh_free(struct pm_mesh *mesh); /* Dump/restore. */ void pm_mesh_struct_dump(const struct pm_mesh *p, FILE *stream); void pm_mesh_struct_restore(struct pm_mesh *p, FILE *stream); -void bkg_cell_mesh_to_gpart_CIC_mapper(void* map_data, int num, void* extra); -void zoom_cell_mesh_to_gpart_CIC_mapper(void* map_data, int num, void* extra); -void bkg_cell_gpart_to_mesh_CIC_mapper(void* map_data, int num, void* extra); -void zoom_cell_gpart_to_mesh_CIC_mapper(void* map_data, int num, void* extra); +void bkg_cell_mesh_to_gpart_CIC_mapper(void *map_data, int num, void *extra); +void zoom_cell_mesh_to_gpart_CIC_mapper(void *map_data, int num, void *extra); +void bkg_cell_gpart_to_mesh_CIC_mapper(void *map_data, int num, void *extra); +void zoom_cell_gpart_to_mesh_CIC_mapper(void *map_data, int num, void *extra); #endif /* SWIFT_MESH_GRAVITY_H */ diff --git a/src/parallel_io.c b/src/parallel_io.c index 4ec4fb2b2a..ca455a08d0 100644 --- a/src/parallel_io.c +++ b/src/parallel_io.c @@ -1626,24 +1626,30 @@ void write_output_parallel(struct engine* e, /* Write the location of the particles in the arrays */ #ifdef WITH_ZOOM_REGION if (e->s->with_zoom_region) { - io_write_cell_offsets(h_grp_cells, e->s->cdim, e->s->zoom_props->cdim, e->s->dim, e->s->cells_top, - e->s->nr_cells, e->s->zoom_props->nr_zoom_cells, e->s->zoom_props->nr_bkg_cells, e->s->width, - e->s->zoom_props->width, mpi_rank, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, offset, - numFields, internal_units, snapshot_units, /*with_zoom=*/1); + io_write_cell_offsets( + h_grp_cells, e->s->cdim, e->s->zoom_props->cdim, e->s->dim, + e->s->cells_top, e->s->nr_cells, e->s->zoom_props->nr_zoom_cells, + e->s->zoom_props->nr_bkg_cells, e->s->width, e->s->zoom_props->width, + mpi_rank, /*distributed=*/0, subsample, subsample_fraction, + e->snapshot_output_count, N_total, offset, numFields, internal_units, + snapshot_units, /*with_zoom=*/1); } else { - io_write_cell_offsets(h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, - e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, e->s->width, - /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, offset, - numFields, internal_units, snapshot_units, /*with_zoom=*/0); + io_write_cell_offsets( + h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, + e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, + e->s->width, + /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, subsample, + subsample_fraction, e->snapshot_output_count, N_total, offset, + numFields, internal_units, snapshot_units, /*with_zoom=*/0); } #else - io_write_cell_offsets(h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, - e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, e->s->width, - /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, offset, - numFields, internal_units, snapshot_units, /*with_zoom=*/0); + io_write_cell_offsets(h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, + e->s->cells_top, e->s->nr_cells, /*nr_zoomcells=*/0, + /*nr_bkgcells=*/e->s->nr_cells, e->s->width, + /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, + subsample, subsample_fraction, e->snapshot_output_count, + N_total, offset, numFields, internal_units, + snapshot_units, /*with_zoom=*/0); #endif /* Close everything */ diff --git a/src/partition.c b/src/partition.c index 063525659c..a8aad5af4c 100644 --- a/src/partition.c +++ b/src/partition.c @@ -106,7 +106,8 @@ static int repart_init_fixed_costs(void); * @param nregions the number of regions * @param samplecells the list of sample cell positions, size of 3*nregions */ -static void pick_vector(struct space *s, int *cdim, int nregions, int *samplecells) { +static void pick_vector(struct space *s, int *cdim, int nregions, + int *samplecells) { /* Get length of space and divide up. */ int length = cdim[0] * cdim[1] * cdim[2]; @@ -128,8 +129,8 @@ static void pick_vector(struct space *s, int *cdim, int nregions, int *samplecel samplecells[m++] = k; l++; } - n++; - if (n == step) n = 0; + n++; + if (n == step) n = 0; } } } @@ -143,7 +144,8 @@ static void pick_vector(struct space *s, int *cdim, int nregions, int *samplecel * Using the sample positions as seeds pick cells that are geometrically * closest and apply the partition to the space. */ -static void split_vector(struct space *s, int *cdim, int nregions, int *samplecells, int offset) { +static void split_vector(struct space *s, int *cdim, int nregions, + int *samplecells, int offset) { int n = 0; for (int i = 0; i < cdim[0]; i++) { @@ -356,8 +358,8 @@ struct counts_mapper_data { else if (parts[k].x[j] >= dim[j]) \ parts[k].x[j] -= dim[j]; \ } \ - const int cid = cell_getid_pos(s, parts[k].x[0], \ - parts[k].x[1], parts[k].x[2]); \ + const int cid = \ + cell_getid_pos(s, parts[k].x[0], parts[k].x[1], parts[k].x[2]); \ if (cid > ucid) ucid = cid; \ if (cid < lcid) lcid = cid; \ } \ @@ -365,8 +367,8 @@ struct counts_mapper_data { if ((lcounts = (double *)calloc(sizeof(double), nused)) == NULL) \ error("Failed to allocate counts thread-specific buffer"); \ for (int k = 0; k < num_elements; k++) { \ - const int cid = cell_getid_pos(s, parts[k].x[0], \ - parts[k].x[1], parts[k].x[2]); \ + const int cid = \ + cell_getid_pos(s, parts[k].x[0], parts[k].x[1], parts[k].x[2]); \ lcounts[cid - lcid] += size; \ } \ for (int k = 0; k < nused; k++) \ @@ -1925,8 +1927,8 @@ void partition_initial_partition(struct partition *initial_partition, for (j = 0; j < 3; j++) { if (s->with_zoom_region) { if (c->tl_cell_type == 3) { - ind[j] = (c->loc[j] - s->zoom_props->region_bounds[2 * j]) - / s->zoom_props->dim[j] * initial_partition->grid[j]; + ind[j] = (c->loc[j] - s->zoom_props->region_bounds[2 * j]) / + s->zoom_props->dim[j] * initial_partition->grid[j]; } else { ind[j] = c->loc[j] / s->dim[j] * initial_partition->grid[j]; } @@ -2052,11 +2054,13 @@ void partition_initial_partition(struct partition *initial_partition, if (s->with_zoom_region) { /* With a zoom region we must apply the background offset */ - split_vector(s, s->cdim, nr_nodes, samplecells, s->zoom_props->tl_cell_offset); + split_vector(s, s->cdim, nr_nodes, samplecells, + s->zoom_props->tl_cell_offset); free(samplecells); int *zoom_samplecells = NULL; - if ((zoom_samplecells = (int *)malloc(sizeof(int) * nr_nodes * 3)) == NULL) + if ((zoom_samplecells = (int *)malloc(sizeof(int) * nr_nodes * 3)) == + NULL) error("Failed to allocate zoom_samplecells"); if (nodeID == 0) { @@ -2064,7 +2068,8 @@ void partition_initial_partition(struct partition *initial_partition, } /* Share the zoom_samplecells around all the nodes. */ - res = MPI_Bcast(zoom_samplecells, nr_nodes * 3, MPI_INT, 0, MPI_COMM_WORLD); + res = + MPI_Bcast(zoom_samplecells, nr_nodes * 3, MPI_INT, 0, MPI_COMM_WORLD); if (res != MPI_SUCCESS) mpi_error(res, "Failed to bcast the partition sample cells."); diff --git a/src/runner_doiact_grav.c b/src/runner_doiact_grav.c index ea45f9a872..5e9b48236b 100644 --- a/src/runner_doiact_grav.c +++ b/src/runner_doiact_grav.c @@ -2489,7 +2489,8 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, #ifdef WITH_ZOOM_REGION const double min_radius2 = cell_min_dist2(top, cj, periodic, dim); #else - const double min_radius2 = cell_min_dist2_same_size(top, cj, periodic, dim); + const double min_radius2 = + cell_min_dist2_same_size(top, cj, periodic, dim); #endif /* Are we beyond the distance where the truncated forces are 0 ?*/ diff --git a/src/runner_others.c b/src/runner_others.c index d93379785b..d5b108f895 100644 --- a/src/runner_others.c +++ b/src/runner_others.c @@ -799,13 +799,18 @@ void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer) { my_id = gp->id_or_neg_offset; error( - "g-particle (id=%lld, type=%s, pos=[%f, %f, %f]) did not interact " + "g-particle (id=%lld, type=%s, pos=[%f, %f, %f]) did not " + "interact " "gravitationally with all other gparts " "gp->num_interacted=%lld, total_gparts=%lld (local " - "num_gparts=%zd inhibited_gparts=%lld cell_type=%d cell_loc=[%f %f %f] cell_width=[%f %f %f] rmax=%f grav_count=%d)", - my_id, part_type_names[gp->type], gp->x[0], gp->x[1], gp->x[2], gp->num_interacted, - e->total_nr_gparts, e->s->nr_gparts, e->count_inhibited_gparts, c->tl_cell_type, - c->loc[0], c->loc[1], c->loc[2], c->width[0], c->width[1], c->width[2], c->grav.multipole->r_max, c->grav.count); + "num_gparts=%zd inhibited_gparts=%lld cell_type=%d " + "cell_loc=[%f %f %f] cell_width=[%f %f %f] rmax=%f " + "grav_count=%d)", + my_id, part_type_names[gp->type], gp->x[0], gp->x[1], gp->x[2], + gp->num_interacted, e->total_nr_gparts, e->s->nr_gparts, + e->count_inhibited_gparts, c->tl_cell_type, c->loc[0], + c->loc[1], c->loc[2], c->width[0], c->width[1], c->width[2], + c->grav.multipole->r_max, c->grav.count); } } #endif diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c index 47b6407918..f340c5afb2 100644 --- a/src/runner_time_integration.c +++ b/src/runner_time_integration.c @@ -455,11 +455,14 @@ void runner_do_kick2(struct runner *r, struct cell *c, const int timer) { #ifdef SWIFT_DEBUG_CHECKS /* Check that kick and the drift are synchronized */ - if (p->ti_drift != p->ti_kick) error("Error integrating part in time. " - "(p->ti_drift=%lld, p->ti_kick=%lld, e->ti_current=%lld, p->x=[%f %f %f]" - "c->tl_cell_type=%d, c->hydro.count=%d)", - p->ti_drift, p->ti_kick, e->ti_current, p->x[0], p->x[1], p->x[2], - c->tl_cell_type, c->hydro.count); + if (p->ti_drift != p->ti_kick) + error( + "Error integrating part in time. " + "(p->ti_drift=%lld, p->ti_kick=%lld, e->ti_current=%lld, " + "p->x=[%f %f %f]" + "c->tl_cell_type=%d, c->hydro.count=%d)", + p->ti_drift, p->ti_kick, e->ti_current, p->x[0], p->x[1], p->x[2], + c->tl_cell_type, c->hydro.count); #endif /* Prepare the values to be drifted */ diff --git a/src/serial_io.c b/src/serial_io.c index bf346324c9..9759880f99 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -1289,24 +1289,30 @@ void write_output_serial(struct engine* e, /* Write the location of the particles in the arrays */ #ifdef WITH_ZOOM_REGION if (e->s->with_zoom_region) { - io_write_cell_offsets(h_grp_cells, e->s->cdim, e->s->zoom_props->cdim, e->s->dim, e->s->cells_top, - e->s->nr_cells, e->s->zoom_props->nr_zoom_cells, e->s->zoom_props->nr_bkg_cells, e->s->width, - e->s->zoom_props->width, mpi_rank, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, offset, - numFields, internal_units, snapshot_units, /*with_zoom=*/1); + io_write_cell_offsets( + h_grp_cells, e->s->cdim, e->s->zoom_props->cdim, e->s->dim, + e->s->cells_top, e->s->nr_cells, e->s->zoom_props->nr_zoom_cells, + e->s->zoom_props->nr_bkg_cells, e->s->width, e->s->zoom_props->width, + mpi_rank, /*distributed=*/0, subsample, subsample_fraction, + e->snapshot_output_count, N_total, offset, numFields, internal_units, + snapshot_units, /*with_zoom=*/1); } else { - io_write_cell_offsets(h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, - e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, e->s->width, - /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, offset, - numFields, internal_units, snapshot_units, /*with_zoom=*/0); + io_write_cell_offsets( + h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, + e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, + e->s->width, + /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, subsample, + subsample_fraction, e->snapshot_output_count, N_total, offset, + numFields, internal_units, snapshot_units, /*with_zoom=*/0); } #else - io_write_cell_offsets(h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, - e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, e->s->width, - /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, offset, - numFields, internal_units, snapshot_units, /*with_zoom=*/0); + io_write_cell_offsets(h_grp_cells, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, + e->s->cells_top, e->s->nr_cells, /*nr_zoomcells=*/0, + /*nr_bkgcells=*/e->s->nr_cells, e->s->width, + /*zoom_width=*/NULL, mpi_rank, /*distributed=*/0, + subsample, subsample_fraction, e->snapshot_output_count, + N_total, offset, numFields, internal_units, + snapshot_units, /*with_zoom=*/0); #endif /* Close everything */ diff --git a/src/single_io.c b/src/single_io.c index 1f6e9ea54f..1900ce57c9 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -1076,27 +1076,33 @@ void write_output_single(struct engine* e, h_grp = H5Gcreate(h_file, "/Cells", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); if (h_grp < 0) error("Error while creating cells group"); - /* Write the location of the particles in the arrays */ + /* Write the location of the particles in the arrays */ #ifdef WITH_ZOOM_REGION if (e->s->with_zoom_region) { - io_write_cell_offsets(h_grp, e->s->cdim, e->s->zoom_props->cdim, e->s->dim, e->s->cells_top, - e->s->nr_cells, e->s->zoom_props->nr_zoom_cells, e->s->zoom_props->nr_bkg_cells, e->s->width, - e->s->zoom_props->width, e->nodeID, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, global_offsets, - numFields, internal_units, snapshot_units, /*with_zoom=*/1); + io_write_cell_offsets( + h_grp, e->s->cdim, e->s->zoom_props->cdim, e->s->dim, e->s->cells_top, + e->s->nr_cells, e->s->zoom_props->nr_zoom_cells, + e->s->zoom_props->nr_bkg_cells, e->s->width, e->s->zoom_props->width, + e->nodeID, /*distributed=*/0, subsample, subsample_fraction, + e->snapshot_output_count, N_total, global_offsets, numFields, + internal_units, snapshot_units, /*with_zoom=*/1); } else { - io_write_cell_offsets(h_grp, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, - e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, e->s->width, - /*zoom_width=*/NULL, e->nodeID, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, global_offsets, - numFields, internal_units, snapshot_units, /*with_zoom=*/0); + io_write_cell_offsets( + h_grp, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, + e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, + e->s->width, + /*zoom_width=*/NULL, e->nodeID, /*distributed=*/0, subsample, + subsample_fraction, e->snapshot_output_count, N_total, global_offsets, + numFields, internal_units, snapshot_units, /*with_zoom=*/0); } #else - io_write_cell_offsets(h_grp, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, e->s->cells_top, - e->s->nr_cells, /*nr_zoomcells=*/0, /*nr_bkgcells=*/e->s->nr_cells, e->s->width, - /*zoom_width=*/NULL, e->nodeID, /*distributed=*/0, subsample, subsample_fraction, - e->snapshot_output_count, N_total, global_offsets, - numFields, internal_units, snapshot_units, /*with_zoom=*/0); + io_write_cell_offsets(h_grp, e->s->cdim, /*zoom_cdim=*/NULL, e->s->dim, + e->s->cells_top, e->s->nr_cells, /*nr_zoomcells=*/0, + /*nr_bkgcells=*/e->s->nr_cells, e->s->width, + /*zoom_width=*/NULL, e->nodeID, /*distributed=*/0, + subsample, subsample_fraction, e->snapshot_output_count, + N_total, global_offsets, numFields, internal_units, + snapshot_units, /*with_zoom=*/0); #endif H5Gclose(h_grp); diff --git a/src/space.c b/src/space.c index 8d46563be9..ae46f2325f 100644 --- a/src/space.c +++ b/src/space.c @@ -44,9 +44,9 @@ #include "atomic.h" #include "const.h" #include "cooling.h" -#include "gravity_properties.h" #include "engine.h" #include "error.h" +#include "gravity_properties.h" #include "kernel_hydro.h" #include "lock.h" #include "minmax.h" @@ -224,11 +224,13 @@ void space_reorder_extras(struct space *s, int verbose) { /* Re-order the gas particles */ if (space_extra_parts) { - /* In the zoom case we need to limit our loop to the cells containing parts (zoom cells). */ + /* In the zoom case we need to limit our loop to the cells containing parts + * (zoom cells). */ if (s->with_zoom_region) { threadpool_map(&s->e->threadpool, space_reorder_extra_parts_mapper, - s->zoom_props->local_zoom_cells_top, s->zoom_props->nr_local_zoom_cells, - sizeof(int), threadpool_auto_chunk_size, s); + s->zoom_props->local_zoom_cells_top, + s->zoom_props->nr_local_zoom_cells, sizeof(int), + threadpool_auto_chunk_size, s); } else { threadpool_map(&s->e->threadpool, space_reorder_extra_parts_mapper, s->local_cells_top, s->nr_local_cells, sizeof(int), @@ -245,11 +247,13 @@ void space_reorder_extras(struct space *s, int verbose) { /* Re-order the star particles */ if (space_extra_sparts) { - /* In the zoom case we need to limit our loop to the cells containing sparts (zoom cells). */ + /* In the zoom case we need to limit our loop to the cells containing sparts + * (zoom cells). */ if (s->with_zoom_region) { threadpool_map(&s->e->threadpool, space_reorder_extra_sparts_mapper, - s->zoom_props->local_zoom_cells_top, s->zoom_props->nr_local_zoom_cells, - sizeof(int), threadpool_auto_chunk_size, s); + s->zoom_props->local_zoom_cells_top, + s->zoom_props->nr_local_zoom_cells, sizeof(int), + threadpool_auto_chunk_size, s); } else { threadpool_map(&s->e->threadpool, space_reorder_extra_sparts_mapper, s->local_cells_top, s->nr_local_cells, sizeof(int), @@ -264,11 +268,13 @@ void space_reorder_extras(struct space *s, int verbose) { /* Re-order the sink particles */ if (space_extra_sinks) { - /* In the zoom case we need to limit our loop to the cells containing sinks (zoom cells). */ + /* In the zoom case we need to limit our loop to the cells containing sinks + * (zoom cells). */ if (s->with_zoom_region) { threadpool_map(&s->e->threadpool, space_reorder_extra_sinks_mapper, - s->zoom_props->local_zoom_cells_top, s->zoom_props->nr_local_zoom_cells, - sizeof(int), threadpool_auto_chunk_size, s); + s->zoom_props->local_zoom_cells_top, + s->zoom_props->nr_local_zoom_cells, sizeof(int), + threadpool_auto_chunk_size, s); } else { threadpool_map(&s->e->threadpool, space_reorder_extra_sinks_mapper, s->local_cells_top, s->nr_local_cells, sizeof(int), @@ -1139,8 +1145,9 @@ void space_collect_mean_masses(struct space *s, int verbose) { */ void space_init(struct space *s, struct swift_params *params, const struct cosmology *cosmo, double dim[3], - const struct hydro_props *hydro_properties, struct gravity_props *gravity_properties, - struct part *parts, struct gpart *gparts, struct sink *sinks, struct spart *sparts, + const struct hydro_props *hydro_properties, + struct gravity_props *gravity_properties, struct part *parts, + struct gpart *gparts, struct sink *sinks, struct spart *sparts, struct bpart *bparts, size_t Npart, size_t Ngpart, size_t Nsink, size_t Nspart, size_t Nbpart, size_t Nnupart, int periodic, int replicate, int remap_ids, int generate_gas_in_ics, @@ -1534,13 +1541,13 @@ void space_init(struct space *s, struct swift_params *params, #ifdef WITH_ZOOM_REGION if (!dry_run) { if (s->with_zoom_region) { - space_regrid_zoom(s, gravity_properties, verbose); - } else { - space_regrid(s, verbose); - } + space_regrid_zoom(s, gravity_properties, verbose); + } else { + space_regrid(s, verbose); + } } #else - if (!dry_run) space_regrid(s, verbose); + if (!dry_run) space_regrid(s, verbose); #endif /* Compute the max id for the generation of unique id. */ @@ -2730,9 +2737,9 @@ void space_write_cell(const struct space *s, FILE *f, const struct cell *c) { fprintf(f, "%g,%g,%i,%i", c->hydro.h_max, c->stars.h_max, c->depth, c->maxdepth); #ifdef WITH_ZOOM_REGION - fprintf(f, ",%i\n", c->tl_cell_type); + fprintf(f, ",%i\n", c->tl_cell_type); #else - fprintf(f, "\n"); + fprintf(f, "\n"); #endif /* Write children */ diff --git a/src/space.h b/src/space.h index d8276485d7..f7fe9594ce 100644 --- a/src/space.h +++ b/src/space.h @@ -380,15 +380,15 @@ struct zoom_region_properties { /*! The ijk integer coordinates of the zoom region's background cell. */ int zoom_cell_ijk[3]; - /*! The number of zoom cells along an axis in a natural top level cell, - * used to define zoom_cdim */ - int nr_zoom_per_bkg_cells; + /*! The number of zoom cells along an axis in a natural top level cell, + * used to define zoom_cdim */ + int nr_zoom_per_bkg_cells; - /*! Number of zoom cells */ - int nr_zoom_cells; + /*! Number of zoom cells */ + int nr_zoom_cells; - /*! Number of natural/background cells */ - int nr_bkg_cells; + /*! Number of natural/background cells */ + int nr_bkg_cells; /*! Number of *local* top-level zoom cells */ int nr_local_zoom_cells; @@ -414,7 +414,8 @@ struct zoom_region_properties { /*! The indices of the *local* top-level background cells */ int *local_bkg_cells_with_particles_top; - /*! Number of particles that have left the zoom region and been converted to dark matter */ + /*! Number of particles that have left the zoom region and been converted to + * dark matter */ size_t nr_wanderers; }; @@ -435,8 +436,9 @@ void space_sinks_sort(struct sink *sinks, int *ind, int *counts, int num_bins, void space_getcells(struct space *s, int nr_cells, struct cell **cells); void space_init(struct space *s, struct swift_params *params, const struct cosmology *cosmo, double dim[3], - const struct hydro_props *hydro_properties, struct gravity_props *gravity_properties, - struct part *parts, struct gpart *gparts, struct sink *sinks, struct spart *sparts, + const struct hydro_props *hydro_properties, + struct gravity_props *gravity_properties, struct part *parts, + struct gpart *gparts, struct sink *sinks, struct spart *sparts, struct bpart *bparts, size_t Npart, size_t Ngpart, size_t Nsink, size_t Nspart, size_t Nbpart, size_t Nnupart, int periodic, int replicate, int remap_ids, int generate_gas_in_ics, @@ -454,7 +456,8 @@ void space_map_parts_xparts(struct space *s, struct cell *c)); void space_map_cells_post(struct space *s, int full, void (*fun)(struct cell *c, void *data), void *data); -void space_rebuild(struct space *s, int repartitioned, struct gravity_props *gravity_properties, int verbose); +void space_rebuild(struct space *s, int repartitioned, + struct gravity_props *gravity_properties, int verbose); void space_recycle(struct space *s, struct cell *c); void space_recycle_list(struct space *s, struct cell *cell_list_begin, struct cell *cell_list_end, diff --git a/src/space_cell_index.c b/src/space_cell_index.c index 01f787a6c5..df9e442e60 100644 --- a/src/space_cell_index.c +++ b/src/space_cell_index.c @@ -124,24 +124,24 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts, if (pos_z == dim_z) pos_z = 0.0; /* Get its cell index */ - const int index = cell_getid_pos(s, pos_x, pos_y, pos_z); + const int index = cell_getid_pos(s, pos_x, pos_y, pos_z); #ifdef WITH_ZOOM_REGION - /* If this part has wandered into a background cell we need to - * convert it into a dark matter particle. */ - if (index >= s->zoom_props->tl_cell_offset) { - /* For this we need to get the xpart and the cell from the space. */ - struct xpart *restrict xp = &xparts[k]; - struct cell *c = &s->cells_top[index]; - cell_convert_part_to_gpart(s->e, c, p, xp); - - /* Increment the number of wanderers */ - s->zoom_props->nr_wanderers++; - } + /* If this part has wandered into a background cell we need to + * convert it into a dark matter particle. */ + if (index >= s->zoom_props->tl_cell_offset) { + /* For this we need to get the xpart and the cell from the space. */ + struct xpart *restrict xp = &xparts[k]; + struct cell *c = &s->cells_top[index]; + cell_convert_part_to_gpart(s->e, c, p, xp); + + /* Increment the number of wanderers */ + s->zoom_props->nr_wanderers++; + } #endif #if defined(SWIFT_DEBUG_CHECKS) && !defined(WITH_ZOOM_REGION) - if (index < 0 || index >= s->cdim[0] * s->cdim[1] * s->cdim[2]) - error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, s->cdim[0], - s->cdim[1], s->cdim[2], pos_x, pos_y, pos_z); + if (index < 0 || index >= s->cdim[0] * s->cdim[1] * s->cdim[2]) + error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, + s->cdim[0], s->cdim[1], s->cdim[2], pos_x, pos_y, pos_z); #endif #ifdef SWIFT_DEBUG_CHECKS if (pos_x >= dim_x || pos_y >= dim_y || pos_z >= dim_z || pos_x < 0. || @@ -261,8 +261,8 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts, #if defined(SWIFT_DEBUG_CHECKS) && !defined(WITH_ZOOM_REGION) if (index < 0 || index >= s->cdim[0] * s->cdim[1] * s->cdim[2]) - error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, s->cdim[0], - s->cdim[1], s->cdim[2], pos_x, pos_y, pos_z); + error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, + s->cdim[0], s->cdim[1], s->cdim[2], pos_x, pos_y, pos_z); #endif #ifdef SWIFT_DEBUG_CHECKS if (pos_x >= dim_x || pos_y >= dim_y || pos_z >= dim_z || pos_x < 0. || @@ -270,7 +270,7 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts, error("Particle outside of simulation box. p->x=[%e %e %e]", pos_x, pos_y, pos_z); #endif - + if (gp->time_bin == time_bin_inhibited) { /* Is this particle to be removed? */ ind[k] = -1; @@ -384,13 +384,13 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts, if (pos_y == dim_y) pos_y = 0.0; if (pos_z == dim_z) pos_z = 0.0; - /* Get its cell index */ - const int index = cell_getid_pos(s, pos_x, pos_y, pos_z); + /* Get its cell index */ + const int index = cell_getid_pos(s, pos_x, pos_y, pos_z); #if defined(SWIFT_DEBUG_CHECKS) && !defined(WITH_ZOOM_REGION) if (index < 0 || index >= s->cdim[0] * s->cdim[1] * s->cdim[2]) - error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, s->cdim[0], - s->cdim[1], s->cdim[2], pos_x, pos_y, pos_z); + error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, + s->cdim[0], s->cdim[1], s->cdim[2], pos_x, pos_y, pos_z); #endif #ifdef SWIFT_DEBUG_CHECKS if (pos_x >= dim_x || pos_y >= dim_y || pos_z >= dim_z || pos_x < 0. || @@ -508,13 +508,13 @@ void space_bparts_get_cell_index_mapper(void *map_data, int nr_bparts, if (pos_y == dim_y) pos_y = 0.0; if (pos_z == dim_z) pos_z = 0.0; - /* Get its cell index */ - const int index = cell_getid_pos(s, pos_x, pos_y, pos_z); + /* Get its cell index */ + const int index = cell_getid_pos(s, pos_x, pos_y, pos_z); #if defined(SWIFT_DEBUG_CHECKS) && !defined(WITH_ZOOM_REGION) if (index < 0 || index >= s->cdim[0] * s->cdim[1] * s->cdim[2]) - error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, s->cdim[0], - s->cdim[1], s->cdim[2], pos_x, pos_y, pos_z); + error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, + s->cdim[0], s->cdim[1], s->cdim[2], pos_x, pos_y, pos_z); #endif #ifdef SWIFT_DEBUG_CHECKS if (pos_x >= dim_x || pos_y >= dim_y || pos_z >= dim_z || pos_x < 0. || @@ -632,13 +632,13 @@ void space_sinks_get_cell_index_mapper(void *map_data, int nr_sinks, if (pos_y == dim_y) pos_y = 0.0; if (pos_z == dim_z) pos_z = 0.0; - /* Get its cell index */ - const int index = cell_getid_pos(s, pos_x, pos_y, pos_z); + /* Get its cell index */ + const int index = cell_getid_pos(s, pos_x, pos_y, pos_z); #if defined(SWIFT_DEBUG_CHECKS) && !defined(WITH_ZOOM_REGION) - if (index < 0 || index >= s->cdim[0] * s->cdim[1] * s->cdim[2]) - error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, s->cdim[0], - s->cdim[1], s->cdim[2], pos_x, pos_y, pos_z); + if (index < 0 || index >= s->cdim[0] * s->cdim[1] * s->cdim[2]) + error("Invalid index=%d cdim=[%d %d %d] p->x=[%e %e %e]", index, + s->cdim[0], s->cdim[1], s->cdim[2], pos_x, pos_y, pos_z); #endif #ifdef SWIFT_DEBUG_CHECKS if (pos_x >= dim_x || pos_y >= dim_y || pos_z >= dim_z || pos_x < 0. || diff --git a/src/space_extras.c b/src/space_extras.c index 59f1ce7304..45fd58955d 100644 --- a/src/space_extras.c +++ b/src/space_extras.c @@ -93,12 +93,13 @@ void space_allocate_extras(struct space *s, int verbose) { ++nr_local_cells; } } - + /* Define a number of cells for non-gparts (baryons). */ size_t nr_baryon_cells = nr_local_cells; #ifdef WITH_ZOOM_REGION - int *local_zoom_cells = (int *)malloc(sizeof(int) * s->zoom_props->nr_zoom_cells); + int *local_zoom_cells = + (int *)malloc(sizeof(int) * s->zoom_props->nr_zoom_cells); if (local_zoom_cells == NULL) error("Failed to allocate list of local top-level cells"); @@ -214,9 +215,12 @@ void space_allocate_extras(struct space *s, int verbose) { 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] + (0.5 * cells[current_cell].width[0]); - s->gparts[i].x[1] = cells[current_cell].loc[1] + (0.5 * cells[current_cell].width[1]); - s->gparts[i].x[2] = cells[current_cell].loc[2] + (0.5 * cells[current_cell].width[2]); + s->gparts[i].x[0] = + cells[current_cell].loc[0] + (0.5 * cells[current_cell].width[0]); + s->gparts[i].x[1] = + cells[current_cell].loc[1] + (0.5 * cells[current_cell].width[1]); + s->gparts[i].x[2] = + cells[current_cell].loc[2] + (0.5 * cells[current_cell].width[2]); ++count_in_cell; count_extra_gparts++; } @@ -230,7 +234,6 @@ void space_allocate_extras(struct space *s, int verbose) { current_cell = local_cells[local_cell_id]; count_in_cell = 0; - } } @@ -320,7 +323,8 @@ void space_allocate_extras(struct space *s, int verbose) { if (count_in_cell == space_extra_parts) { ++local_cell_id; - /* When running with a zoom region we only ever need to consider the zoom cells. */ + /* When running with a zoom region we only ever need to consider the + * zoom cells. */ if (s->with_zoom_region) { if (local_cell_id == nr_baryon_cells) break; @@ -410,7 +414,8 @@ void space_allocate_extras(struct space *s, int verbose) { if (count_in_cell == space_extra_sinks) { ++local_cell_id; - /* When running with a zoom region we only ever need to consider the zoom cells. */ + /* When running with a zoom region we only ever need to consider the + * zoom cells. */ if (s->with_zoom_region) { if (local_cell_id == nr_baryon_cells) break; @@ -501,7 +506,8 @@ void space_allocate_extras(struct space *s, int verbose) { if (count_in_cell == space_extra_sparts) { ++local_cell_id; - /* When running with a zoom region we only ever need to consider the zoom cells. */ + /* When running with a zoom region we only ever need to consider the + * zoom cells. */ if (s->with_zoom_region) { if (local_cell_id == nr_baryon_cells) break; @@ -592,7 +598,8 @@ void space_allocate_extras(struct space *s, int verbose) { if (count_in_cell == space_extra_bparts) { ++local_cell_id; - /* When running with a zoom region we only ever need to consider the zoom cells. */ + /* When running with a zoom region we only ever need to consider the + * zoom cells. */ if (s->with_zoom_region) { if (local_cell_id == nr_baryon_cells) break; diff --git a/src/space_rebuild.c b/src/space_rebuild.c index f01e5d0cf5..3d47c776a9 100644 --- a/src/space_rebuild.c +++ b/src/space_rebuild.c @@ -48,7 +48,8 @@ extern unsigned long long last_leaf_cell_id; * @param repartitioned Did we just repartition? * @param verbose Print messages to stdout or not */ -void space_rebuild(struct space *s, int repartitioned, struct gravity_props *gravity_properties, int verbose) { +void space_rebuild(struct space *s, int repartitioned, + struct gravity_props *gravity_properties, int verbose) { const ticks tic = getticks(); @@ -65,11 +66,11 @@ void space_rebuild(struct space *s, int repartitioned, struct gravity_props *gra /* Re-grid if necessary, or just re-set the cell data. */ #ifdef WITH_ZOOM_REGION - if (s->with_zoom_region) { + if (s->with_zoom_region) { space_regrid_zoom(s, gravity_properties, verbose); - } else { - space_regrid(s, verbose); - } + } else { + space_regrid(s, verbose); + } #else space_regrid(s, verbose); #endif @@ -867,8 +868,8 @@ void space_rebuild(struct space *s, int repartitioned, struct gravity_props *gra if (gp->x[0] < c->loc[0] || gp->x[0] > c->loc[0] + c->width[0] || gp->x[1] < c->loc[1] || gp->x[1] > c->loc[1] + c->width[1] || gp->x[2] < c->loc[2] || gp->x[2] > c->loc[2] + c->width[2]) { - error("gpart not sorted into the right top-level cell!"); - } + error("gpart not sorted into the right top-level cell!"); + } } #endif /* SWIFT_DEBUG_CHECKS */ @@ -897,11 +898,11 @@ void space_rebuild(struct space *s, int repartitioned, struct gravity_props *gra verbose); #endif - /* Define variables to count particles in cell types */ - int bkg_cell_particles = 0; - int zoom_cell_particles = 0; + /* Define variables to count particles in cell types */ + int bkg_cell_particles = 0; + int zoom_cell_particles = 0; - /* Hook the cells up to the parts. Make list of local and non-empty cells */ + /* Hook the cells up to the parts. Make list of local and non-empty cells */ const ticks tic3 = getticks(); struct part *finger = s->parts; struct xpart *xfinger = s->xparts; @@ -959,9 +960,13 @@ void space_rebuild(struct space *s, int repartitioned, struct gravity_props *gra if (s->with_zoom_region) { /* Add the number of particles to the cell counter */ if (c->tl_cell_type <= 2) { - bkg_cell_particles += (c->hydro.count + c->grav.count + c->stars.count + c->sinks.count + c->black_holes.count); + bkg_cell_particles += + (c->hydro.count + c->grav.count + c->stars.count + + c->sinks.count + c->black_holes.count); } else { - zoom_cell_particles += (c->hydro.count + c->grav.count + c->stars.count + c->sinks.count + c->black_holes.count); + zoom_cell_particles += + (c->hydro.count + c->grav.count + c->stars.count + + c->sinks.count + c->black_holes.count); } } #endif @@ -982,14 +987,16 @@ void space_rebuild(struct space *s, int repartitioned, struct gravity_props *gra /* Add this cell to the appropriate list of local cells */ if (c->tl_cell_type == zoom_tl_cell) { - s->zoom_props->local_zoom_cells_top[s->zoom_props->nr_local_zoom_cells] = k; + s->zoom_props + ->local_zoom_cells_top[s->zoom_props->nr_local_zoom_cells] = k; s->zoom_props->nr_local_zoom_cells++; - } else if (c->tl_cell_type == tl_cell || c->tl_cell_type == tl_cell_neighbour) { + } else if (c->tl_cell_type == tl_cell || + c->tl_cell_type == tl_cell_neighbour) { - s->zoom_props->local_bkg_cells_top[s->zoom_props->nr_local_bkg_cells] = k; + s->zoom_props + ->local_bkg_cells_top[s->zoom_props->nr_local_bkg_cells] = k; s->zoom_props->nr_local_bkg_cells++; - } } #endif @@ -1006,14 +1013,16 @@ void space_rebuild(struct space *s, int repartitioned, struct gravity_props *gra /* Add this cell to the appropriate list of local cells */ if (c->tl_cell_type == zoom_tl_cell) { - s->zoom_props->local_zoom_cells_with_particles_top[s->zoom_props->nr_local_zoom_cells_with_particles] = k; + s->zoom_props->local_zoom_cells_with_particles_top + [s->zoom_props->nr_local_zoom_cells_with_particles] = k; s->zoom_props->nr_local_zoom_cells_with_particles++; - } else if (c->tl_cell_type == tl_cell || c->tl_cell_type == tl_cell_neighbour) { + } else if (c->tl_cell_type == tl_cell || + c->tl_cell_type == tl_cell_neighbour) { - s->zoom_props->local_bkg_cells_with_particles_top[s->zoom_props->nr_local_bkg_cells_with_particles] = k; + s->zoom_props->local_bkg_cells_with_particles_top + [s->zoom_props->nr_local_bkg_cells_with_particles] = k; s->zoom_props->nr_local_bkg_cells_with_particles++; - } } #endif @@ -1030,14 +1039,14 @@ void space_rebuild(struct space *s, int repartitioned, struct gravity_props *gra if (s->with_zoom_region) { message("Have %d local particles in background cells", bkg_cell_particles); - message("Have %d local particles in zoom cells", - zoom_cell_particles); + message("Have %d local particles in zoom cells", zoom_cell_particles); } #endif - /* Lets report how many wanderers have been removed */ - if (s->with_hydro && s->zoom_props->nr_wanderers > 0) - message("Converted %zu wandering particles to dark matter thus far", s->zoom_props->nr_wanderers); + /* Lets report how many wanderers have been removed */ + if (s->with_hydro && s->zoom_props->nr_wanderers > 0) + message("Converted %zu wandering particles to dark matter thus far", + s->zoom_props->nr_wanderers); } /* Re-order the extra particles such that they are at the end of their cell's diff --git a/src/space_regrid.c b/src/space_regrid.c index e24c1af135..9c828b82c4 100644 --- a/src/space_regrid.c +++ b/src/space_regrid.c @@ -38,76 +38,76 @@ */ void space_regrid(struct space *s, int verbose) { - const size_t nr_parts = s->nr_parts; - const size_t nr_sparts = s->nr_sparts; - const size_t nr_bparts = s->nr_bparts; - const size_t nr_sinks = s->nr_sinks; - const ticks tic = getticks(); - const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0; - - /* Run through the cells and get the current h_max. */ - // tic = getticks(); - float h_max = s->cell_min / kernel_gamma / space_stretch; - if (nr_parts > 0) { - - /* Can we use the list of local non-empty top-level cells? */ - if (s->local_cells_with_particles_top != NULL) { - for (int k = 0; k < s->nr_local_cells_with_particles; ++k) { - const struct cell *c = - &s->cells_top[s->local_cells_with_particles_top[k]]; - if (c->hydro.h_max > h_max) { - h_max = c->hydro.h_max; - } - if (c->stars.h_max > h_max) { - h_max = c->stars.h_max; - } - if (c->black_holes.h_max > h_max) { - h_max = c->black_holes.h_max; - } - if (c->sinks.r_cut_max > h_max) { - h_max = c->sinks.r_cut_max / kernel_gamma; - } - } - - /* Can we instead use all the top-level cells? */ - } else if (s->cells_top != NULL) { - for (int k = 0; k < s->nr_cells; k++) { - const struct cell *c = &s->cells_top[k]; - if (c->nodeID == engine_rank && c->hydro.h_max > h_max) { - h_max = c->hydro.h_max; - } - if (c->nodeID == engine_rank && c->stars.h_max > h_max) { - h_max = c->stars.h_max; - } - if (c->nodeID == engine_rank && c->black_holes.h_max > h_max) { - h_max = c->black_holes.h_max; - } - if (c->nodeID == engine_rank && c->sinks.r_cut_max > h_max) { - h_max = c->sinks.r_cut_max / kernel_gamma; - } - } - - /* Last option: run through the particles */ - } else { - for (size_t k = 0; k < nr_parts; k++) { - if (s->parts[k].h > h_max) h_max = s->parts[k].h; - } - for (size_t k = 0; k < nr_sparts; k++) { - if (s->sparts[k].h > h_max) h_max = s->sparts[k].h; - } - for (size_t k = 0; k < nr_bparts; k++) { - if (s->bparts[k].h > h_max) h_max = s->bparts[k].h; - } - for (size_t k = 0; k < nr_sinks; k++) { - if (s->sinks[k].r_cut > h_max) h_max = s->sinks[k].r_cut / kernel_gamma; - } - } - } + const size_t nr_parts = s->nr_parts; + const size_t nr_sparts = s->nr_sparts; + const size_t nr_bparts = s->nr_bparts; + const size_t nr_sinks = s->nr_sinks; + const ticks tic = getticks(); + const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0; + + /* Run through the cells and get the current h_max. */ + // tic = getticks(); + float h_max = s->cell_min / kernel_gamma / space_stretch; + if (nr_parts > 0) { + + /* Can we use the list of local non-empty top-level cells? */ + if (s->local_cells_with_particles_top != NULL) { + for (int k = 0; k < s->nr_local_cells_with_particles; ++k) { + const struct cell *c = + &s->cells_top[s->local_cells_with_particles_top[k]]; + if (c->hydro.h_max > h_max) { + h_max = c->hydro.h_max; + } + if (c->stars.h_max > h_max) { + h_max = c->stars.h_max; + } + if (c->black_holes.h_max > h_max) { + h_max = c->black_holes.h_max; + } + if (c->sinks.r_cut_max > h_max) { + h_max = c->sinks.r_cut_max / kernel_gamma; + } + } + + /* Can we instead use all the top-level cells? */ + } else if (s->cells_top != NULL) { + for (int k = 0; k < s->nr_cells; k++) { + const struct cell *c = &s->cells_top[k]; + if (c->nodeID == engine_rank && c->hydro.h_max > h_max) { + h_max = c->hydro.h_max; + } + if (c->nodeID == engine_rank && c->stars.h_max > h_max) { + h_max = c->stars.h_max; + } + if (c->nodeID == engine_rank && c->black_holes.h_max > h_max) { + h_max = c->black_holes.h_max; + } + if (c->nodeID == engine_rank && c->sinks.r_cut_max > h_max) { + h_max = c->sinks.r_cut_max / kernel_gamma; + } + } + + /* Last option: run through the particles */ + } else { + for (size_t k = 0; k < nr_parts; k++) { + if (s->parts[k].h > h_max) h_max = s->parts[k].h; + } + for (size_t k = 0; k < nr_sparts; k++) { + if (s->sparts[k].h > h_max) h_max = s->sparts[k].h; + } + for (size_t k = 0; k < nr_bparts; k++) { + if (s->bparts[k].h > h_max) h_max = s->bparts[k].h; + } + for (size_t k = 0; k < nr_sinks; k++) { + if (s->sinks[k].r_cut > h_max) h_max = s->sinks[k].r_cut / kernel_gamma; + } + } + } /* If we are running in parallel, make sure everybody agrees on how large the largest cell should be. */ #ifdef WITH_MPI - { + { float buff; if (MPI_Allreduce(&h_max, &buff, 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD) != MPI_SUCCESS) @@ -155,7 +155,7 @@ void space_regrid(struct space *s, int verbose) { * global partition is recomputed and the particles redistributed. * Be prepared to do that. */ #ifdef WITH_MPI - double oldwidth[3] = {0., 0., 0.}; + double oldwidth[3] = {0., 0., 0.}; double oldcdim[3] = {0., 0., 0.}; int *oldnodeIDs = NULL; if (cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] || cdim[2] < s->cdim[2]) { @@ -188,157 +188,157 @@ void space_regrid(struct space *s, int verbose) { const int no_regrid = (s->cells_top == NULL && oldnodeIDs == NULL); #endif - /* Do we need to re-build the upper-level cells? */ - // tic = getticks(); - if (s->cells_top == NULL || cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] || - cdim[2] < s->cdim[2]) { + /* Do we need to re-build the upper-level cells? */ + // tic = getticks(); + if (s->cells_top == NULL || cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] || + cdim[2] < s->cdim[2]) { /* Be verbose about this. */ #ifdef SWIFT_DEBUG_CHECKS - message("(re)griding space cdim=(%d %d %d)", cdim[0], cdim[1], cdim[2]); + message("(re)griding space cdim=(%d %d %d)", cdim[0], cdim[1], cdim[2]); fflush(stdout); #endif - /* Free the old cells, if they were allocated. */ - if (s->cells_top != NULL) { - space_free_cells(s); - swift_free("local_cells_with_tasks_top", s->local_cells_with_tasks_top); - swift_free("local_cells_top", s->local_cells_top); - swift_free("cells_with_particles_top", s->cells_with_particles_top); - swift_free("local_cells_with_particles_top", - s->local_cells_with_particles_top); - swift_free("cells_top", s->cells_top); - swift_free("multipoles_top", s->multipoles_top); - } - - /* Also free the task arrays, these will be regenerated and we can use the - * memory while copying the particle arrays. */ - if (s->e != NULL) scheduler_free_tasks(&s->e->sched); - - /* Set the new cell dimensions only if smaller. */ - for (int k = 0; k < 3; k++) { - s->cdim[k] = cdim[k]; - s->width[k] = s->dim[k] / cdim[k]; - s->iwidth[k] = 1.0 / s->width[k]; - } - const float dmin = min3(s->width[0], s->width[1], s->width[2]); - - /* Allocate the highest level of cells. */ - s->tot_cells = s->nr_cells = cdim[0] * cdim[1] * cdim[2]; - - if (swift_memalign("cells_top", (void **)&s->cells_top, cell_align, - s->nr_cells * sizeof(struct cell)) != 0) - error("Failed to allocate top-level cells."); - bzero(s->cells_top, s->nr_cells * sizeof(struct cell)); - - /* Allocate the multipoles for the top-level cells. */ - if (s->with_self_gravity) { - if (swift_memalign("multipoles_top", (void **)&s->multipoles_top, - multipole_align, - s->nr_cells * sizeof(struct gravity_tensors)) != 0) - error("Failed to allocate top-level multipoles."); - bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors)); - } - - /* Allocate the indices of local cells */ - if (swift_memalign("local_cells_top", (void **)&s->local_cells_top, - SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) - error("Failed to allocate indices of local top-level cells."); - bzero(s->local_cells_top, s->nr_cells * sizeof(int)); - - /* Allocate the indices of local cells with tasks */ - if (swift_memalign("local_cells_with_tasks_top", - (void **)&s->local_cells_with_tasks_top, - SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) - error("Failed to allocate indices of local top-level cells with tasks."); - bzero(s->local_cells_with_tasks_top, s->nr_cells * sizeof(int)); - - /* Allocate the indices of cells with particles */ - if (swift_memalign("cells_with_particles_top", - (void **)&s->cells_with_particles_top, - SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) - error("Failed to allocate indices of top-level cells with particles."); - bzero(s->cells_with_particles_top, s->nr_cells * sizeof(int)); - - /* Allocate the indices of local cells with particles */ - if (swift_memalign("local_cells_with_particles_top", - (void **)&s->local_cells_with_particles_top, - SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) - error( - "Failed to allocate indices of local top-level cells with " - "particles."); - bzero(s->local_cells_with_particles_top, s->nr_cells * sizeof(int)); - - /* Set the cells' locks */ - for (int k = 0; k < s->nr_cells; k++) { - if (lock_init(&s->cells_top[k].hydro.lock) != 0) - error("Failed to init spinlock for hydro."); - if (lock_init(&s->cells_top[k].grav.plock) != 0) - error("Failed to init spinlock for gravity."); - if (lock_init(&s->cells_top[k].grav.mlock) != 0) - error("Failed to init spinlock for multipoles."); - if (lock_init(&s->cells_top[k].grav.star_formation_lock) != 0) - error("Failed to init spinlock for star formation (gpart)."); - if (lock_init(&s->cells_top[k].stars.lock) != 0) - error("Failed to init spinlock for stars."); - if (lock_init(&s->cells_top[k].sinks.lock) != 0) - error("Failed to init spinlock for sinks."); - if (lock_init(&s->cells_top[k].sinks.sink_formation_lock) != 0) - error("Failed to init spinlock for sink formation."); - if (lock_init(&s->cells_top[k].black_holes.lock) != 0) - error("Failed to init spinlock for black holes."); - if (lock_init(&s->cells_top[k].stars.star_formation_lock) != 0) - error("Failed to init spinlock for star formation (spart)."); - } - - /* Set the cell location and sizes. */ - for (int i = 0; i < cdim[0]; i++) - for (int j = 0; j < cdim[1]; j++) - for (int k = 0; k < cdim[2]; k++) { - const size_t cid = cell_getid(cdim, i, j, k); - struct cell *restrict c = &s->cells_top[cid]; - c->loc[0] = i * s->width[0]; - c->loc[1] = j * s->width[1]; - c->loc[2] = k * s->width[2]; - c->width[0] = s->width[0]; - c->width[1] = s->width[1]; - c->width[2] = s->width[2]; - c->dmin = dmin; + /* Free the old cells, if they were allocated. */ + if (s->cells_top != NULL) { + space_free_cells(s); + swift_free("local_cells_with_tasks_top", s->local_cells_with_tasks_top); + swift_free("local_cells_top", s->local_cells_top); + swift_free("cells_with_particles_top", s->cells_with_particles_top); + swift_free("local_cells_with_particles_top", + s->local_cells_with_particles_top); + swift_free("cells_top", s->cells_top); + swift_free("multipoles_top", s->multipoles_top); + } + + /* Also free the task arrays, these will be regenerated and we can use the + * memory while copying the particle arrays. */ + if (s->e != NULL) scheduler_free_tasks(&s->e->sched); + + /* Set the new cell dimensions only if smaller. */ + for (int k = 0; k < 3; k++) { + s->cdim[k] = cdim[k]; + s->width[k] = s->dim[k] / cdim[k]; + s->iwidth[k] = 1.0 / s->width[k]; + } + const float dmin = min3(s->width[0], s->width[1], s->width[2]); + + /* Allocate the highest level of cells. */ + s->tot_cells = s->nr_cells = cdim[0] * cdim[1] * cdim[2]; + + if (swift_memalign("cells_top", (void **)&s->cells_top, cell_align, + s->nr_cells * sizeof(struct cell)) != 0) + error("Failed to allocate top-level cells."); + bzero(s->cells_top, s->nr_cells * sizeof(struct cell)); + + /* Allocate the multipoles for the top-level cells. */ + if (s->with_self_gravity) { + if (swift_memalign("multipoles_top", (void **)&s->multipoles_top, + multipole_align, + s->nr_cells * sizeof(struct gravity_tensors)) != 0) + error("Failed to allocate top-level multipoles."); + bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors)); + } + + /* Allocate the indices of local cells */ + if (swift_memalign("local_cells_top", (void **)&s->local_cells_top, + SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) + error("Failed to allocate indices of local top-level cells."); + bzero(s->local_cells_top, s->nr_cells * sizeof(int)); + + /* Allocate the indices of local cells with tasks */ + if (swift_memalign("local_cells_with_tasks_top", + (void **)&s->local_cells_with_tasks_top, + SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) + error("Failed to allocate indices of local top-level cells with tasks."); + bzero(s->local_cells_with_tasks_top, s->nr_cells * sizeof(int)); + + /* Allocate the indices of cells with particles */ + if (swift_memalign("cells_with_particles_top", + (void **)&s->cells_with_particles_top, + SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) + error("Failed to allocate indices of top-level cells with particles."); + bzero(s->cells_with_particles_top, s->nr_cells * sizeof(int)); + + /* Allocate the indices of local cells with particles */ + if (swift_memalign("local_cells_with_particles_top", + (void **)&s->local_cells_with_particles_top, + SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) + error( + "Failed to allocate indices of local top-level cells with " + "particles."); + bzero(s->local_cells_with_particles_top, s->nr_cells * sizeof(int)); + + /* Set the cells' locks */ + for (int k = 0; k < s->nr_cells; k++) { + if (lock_init(&s->cells_top[k].hydro.lock) != 0) + error("Failed to init spinlock for hydro."); + if (lock_init(&s->cells_top[k].grav.plock) != 0) + error("Failed to init spinlock for gravity."); + if (lock_init(&s->cells_top[k].grav.mlock) != 0) + error("Failed to init spinlock for multipoles."); + if (lock_init(&s->cells_top[k].grav.star_formation_lock) != 0) + error("Failed to init spinlock for star formation (gpart)."); + if (lock_init(&s->cells_top[k].stars.lock) != 0) + error("Failed to init spinlock for stars."); + if (lock_init(&s->cells_top[k].sinks.lock) != 0) + error("Failed to init spinlock for sinks."); + if (lock_init(&s->cells_top[k].sinks.sink_formation_lock) != 0) + error("Failed to init spinlock for sink formation."); + if (lock_init(&s->cells_top[k].black_holes.lock) != 0) + error("Failed to init spinlock for black holes."); + if (lock_init(&s->cells_top[k].stars.star_formation_lock) != 0) + error("Failed to init spinlock for star formation (spart)."); + } + + /* Set the cell location and sizes. */ + for (int i = 0; i < cdim[0]; i++) + for (int j = 0; j < cdim[1]; j++) + for (int k = 0; k < cdim[2]; k++) { + const size_t cid = cell_getid(cdim, i, j, k); + struct cell *restrict c = &s->cells_top[cid]; + c->loc[0] = i * s->width[0]; + c->loc[1] = j * s->width[1]; + c->loc[2] = k * s->width[2]; + c->width[0] = s->width[0]; + c->width[1] = s->width[1]; + c->width[2] = s->width[2]; + c->dmin = dmin; c->tl_cell_type = tl_cell; - c->depth = 0; - c->split = 0; - c->hydro.count = 0; - c->grav.count = 0; - c->stars.count = 0; - c->sinks.count = 0; - c->top = c; - c->super = c; - c->hydro.super = c; - c->grav.super = c; - c->hydro.ti_old_part = ti_current; - c->grav.ti_old_part = ti_current; - c->stars.ti_old_part = ti_current; - c->sinks.ti_old_part = ti_current; - c->black_holes.ti_old_part = ti_current; - c->grav.ti_old_multipole = ti_current; + c->depth = 0; + c->split = 0; + c->hydro.count = 0; + c->grav.count = 0; + c->stars.count = 0; + c->sinks.count = 0; + c->top = c; + c->super = c; + c->hydro.super = c; + c->grav.super = c; + c->hydro.ti_old_part = ti_current; + c->grav.ti_old_part = ti_current; + c->stars.ti_old_part = ti_current; + c->sinks.ti_old_part = ti_current; + c->black_holes.ti_old_part = ti_current; + c->grav.ti_old_multipole = ti_current; #ifdef WITH_MPI - c->mpi.tag = -1; + c->mpi.tag = -1; c->mpi.recv = NULL; c->mpi.send = NULL; #endif // WITH_MPI - if (s->with_self_gravity) c->grav.multipole = &s->multipoles_top[cid]; + if (s->with_self_gravity) c->grav.multipole = &s->multipoles_top[cid]; #if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_CELL_GRAPH) - cell_assign_top_level_cell_index(c, s); + cell_assign_top_level_cell_index(c, s); #endif - } + } - /* Be verbose about the change. */ - if (verbose) - message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1], - cdim[2]); + /* Be verbose about the change. */ + if (verbose) + message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1], + cdim[2]); #ifdef WITH_MPI - if (oldnodeIDs != NULL) { + if (oldnodeIDs != NULL) { /* We have changed the top-level cell dimension, so need to redistribute * cells around the nodes. We repartition using the old space node * positions as a grid to resample. */ @@ -387,17 +387,17 @@ void space_regrid(struct space *s, int verbose) { } #endif /* WITH_MPI */ - // message( "rebuilding upper-level cells took %.3f %s." , - // clocks_from_ticks(double)(getticks() - tic), clocks_getunit()); + // message( "rebuilding upper-level cells took %.3f %s." , + // clocks_from_ticks(double)(getticks() - tic), clocks_getunit()); - } /* re-build upper-level cells? */ - else { /* Otherwise, just clean up the cells. */ + } /* re-build upper-level cells? */ + else { /* Otherwise, just clean up the cells. */ - /* Free the old cells, if they were allocated. */ - space_free_cells(s); - } + /* Free the old cells, if they were allocated. */ + space_free_cells(s); + } - if (verbose) - message("took %.3f %s.", clocks_from_ticks(getticks() - tic), - clocks_getunit()); + if (verbose) + message("took %.3f %s.", clocks_from_ticks(getticks() - tic), + clocks_getunit()); } diff --git a/src/space_split.c b/src/space_split.c index 08d72afe44..f57e126697 100644 --- a/src/space_split.c +++ b/src/space_split.c @@ -720,8 +720,8 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) { #ifdef WITH_ZOOM_REGION /** - * @brief A wrapper for #threadpool mapper function to split background cells if they contain - * too many particles. + * @brief A wrapper for #threadpool mapper function to split background cells if + * they contain too many particles. * * @param map_data Pointer towards the top-cells. * @param num_cells The number of cells to treat. @@ -732,8 +732,8 @@ void bkg_space_split_mapper(void *map_data, int num_cells, void *extra_data) { } /** - * @brief A wrapper for #threadpool mapper function to split zoom cells if they contain - * too many particles. + * @brief A wrapper for #threadpool mapper function to split zoom cells if they + * contain too many particles. * * @param map_data Pointer towards the top-cells. * @param num_cells The number of cells to treat. @@ -760,20 +760,20 @@ void space_split(struct space *s, int verbose) { #ifdef WITH_ZOOM_REGION if (s->with_zoom_region) { threadpool_map(&s->e->threadpool, bkg_space_split_mapper, - s->zoom_props->local_bkg_cells_with_particles_top, - s->zoom_props->nr_local_bkg_cells_with_particles, sizeof(int), - threadpool_auto_chunk_size, s); + s->zoom_props->local_bkg_cells_with_particles_top, + s->zoom_props->nr_local_bkg_cells_with_particles, + sizeof(int), threadpool_auto_chunk_size, s); threadpool_map(&s->e->threadpool, zoom_space_split_mapper, - s->zoom_props->local_zoom_cells_with_particles_top, - s->zoom_props->nr_local_zoom_cells_with_particles, sizeof(int), - threadpool_auto_chunk_size, s); - } else{ + s->zoom_props->local_zoom_cells_with_particles_top, + s->zoom_props->nr_local_zoom_cells_with_particles, + sizeof(int), threadpool_auto_chunk_size, s); + } else { threadpool_map(&s->e->threadpool, space_split_mapper, - s->local_cells_with_particles_top, - s->nr_local_cells_with_particles, sizeof(int), - threadpool_auto_chunk_size, s); + s->local_cells_with_particles_top, + s->nr_local_cells_with_particles, sizeof(int), + threadpool_auto_chunk_size, s); } -# else +#else threadpool_map(&s->e->threadpool, space_split_mapper, s->local_cells_with_particles_top, s->nr_local_cells_with_particles, sizeof(int), diff --git a/src/task.c b/src/task.c index 486315255e..070ee2069f 100644 --- a/src/task.c +++ b/src/task.c @@ -1174,7 +1174,8 @@ void task_print(const struct task *t) { */ void task_get_group_name(int type, int subtype, char *cluster) { - if (type == task_type_grav_long_range || type == task_type_grav_long_range_bkg || type == task_type_grav_mm) { + if (type == task_type_grav_long_range || + type == task_type_grav_long_range_bkg || type == task_type_grav_mm) { strcpy(cluster, "Gravity"); return; diff --git a/src/zoom_partition.c b/src/zoom_partition.c index 5277ef7001..778d9b8dda 100644 --- a/src/zoom_partition.c +++ b/src/zoom_partition.c @@ -23,17 +23,16 @@ /* Config parameters. */ #include "../config.h" - -#include <float.h> - #include "cell.h" -#include "gravity_properties.h" #include "engine.h" -#include "proxy.h" +#include "gravity_properties.h" #include "partition.h" +#include "proxy.h" #include "space.h" #include "zoom_region.h" +#include <float.h> + /* MPI headers. */ #ifdef WITH_MPI #include <mpi.h> @@ -50,27 +49,27 @@ */ static int check_complete(struct space *s, int verbose, int nregions) { - int *present = NULL; - if ((present = (int *)malloc(sizeof(int) * nregions)) == NULL) - error("Failed to allocate present array"); + int *present = NULL; + if ((present = (int *)malloc(sizeof(int) * nregions)) == NULL) + error("Failed to allocate present array"); - int failed = 0; - for (int i = 0; i < nregions; i++) present[i] = 0; - for (int i = 0; i < s->nr_cells; i++) { - if (s->cells_top[i].nodeID <= nregions) - present[s->cells_top[i].nodeID]++; - else - message("Bad nodeID: s->cells_top[%d].nodeID = %d", i, - s->cells_top[i].nodeID); - } - for (int i = 0; i < nregions; i++) { - if (!present[i]) { - failed = 1; - if (verbose) message("Region %d is not present in partition", i); - } - } - free(present); - return (!failed); + int failed = 0; + for (int i = 0; i < nregions; i++) present[i] = 0; + for (int i = 0; i < s->nr_cells; i++) { + if (s->cells_top[i].nodeID <= nregions) + present[s->cells_top[i].nodeID]++; + else + message("Bad nodeID: s->cells_top[%d].nodeID = %d", i, + s->cells_top[i].nodeID); + } + for (int i = 0; i < nregions; i++) { + if (!present[i]) { + failed = 1; + if (verbose) message("Region %d is not present in partition", i); + } + } + free(present); + return (!failed); } /** @@ -96,48 +95,52 @@ static int check_complete(struct space *s, int verbose, int nregions) { * * @return 1 if the new space contains nodeIDs from all nodes, 0 otherwise. */ -int partition_space_to_space_zoom(double *oldh, double *oldcdim, double *oldzoomh, - double *oldzoomcdim, int *oldnodeIDs, struct space *s) { +int partition_space_to_space_zoom(double *oldh, double *oldcdim, + double *oldzoomh, double *oldzoomcdim, + int *oldnodeIDs, struct space *s) { - /* Define the old tl_cell_offset */ - const int old_bkg_cell_offset = oldzoomcdim[0] * oldzoomcdim[1] * oldzoomcdim[2]; + /* Define the old tl_cell_offset */ + const int old_bkg_cell_offset = + oldzoomcdim[0] * oldzoomcdim[1] * oldzoomcdim[2]; - /* Loop over all the new zoom cells. */ - for (int i = 0; i < s->zoom_props->cdim[0]; i++) { - for (int j = 0; j < s->zoom_props->cdim[1]; j++) { - for (int k = 0; k < s->zoom_props->cdim[2]; k++) { + /* Loop over all the new zoom cells. */ + for (int i = 0; i < s->zoom_props->cdim[0]; i++) { + for (int j = 0; j < s->zoom_props->cdim[1]; j++) { + for (int k = 0; k < s->zoom_props->cdim[2]; k++) { - /* Scale indices to old cell space. */ - const int ii = rint(i * s->zoom_props->iwidth[0] * oldzoomh[0]); - const int jj = rint(j * s->zoom_props->iwidth[1] * oldzoomh[1]); - const int kk = rint(k * s->zoom_props->iwidth[2] * oldzoomh[2]); + /* Scale indices to old cell space. */ + const int ii = rint(i * s->zoom_props->iwidth[0] * oldzoomh[0]); + const int jj = rint(j * s->zoom_props->iwidth[1] * oldzoomh[1]); + const int kk = rint(k * s->zoom_props->iwidth[2] * oldzoomh[2]); - const int cid = cell_getid(s->zoom_props->cdim, i, j, k); - const int oldcid = cell_getid(oldzoomcdim, ii, jj, kk); - s->cells_top[cid].nodeID = oldnodeIDs[oldcid]; - } - } - } + const int cid = cell_getid(s->zoom_props->cdim, i, j, k); + const int oldcid = cell_getid(oldzoomcdim, ii, jj, kk); + s->cells_top[cid].nodeID = oldnodeIDs[oldcid]; + } + } + } - /* Loop over all the new cells. */ - for (int i = 0; i < s->cdim[0]; i++) { - for (int j = 0; j < s->cdim[1]; j++) { - for (int k = 0; k < s->cdim[2]; k++) { + /* Loop over all the new cells. */ + for (int i = 0; i < s->cdim[0]; i++) { + for (int j = 0; j < s->cdim[1]; j++) { + for (int k = 0; k < s->cdim[2]; k++) { - /* Scale indices to old cell space. */ - const int ii = rint(i * s->iwidth[0] * oldh[0]); - const int jj = rint(j * s->iwidth[1] * oldh[1]); - const int kk = rint(k * s->iwidth[2] * oldh[2]); + /* Scale indices to old cell space. */ + const int ii = rint(i * s->iwidth[0] * oldh[0]); + const int jj = rint(j * s->iwidth[1] * oldh[1]); + const int kk = rint(k * s->iwidth[2] * oldh[2]); - const int cid = cell_getid(s->cdim, i, j, k) + s->zoom_props->tl_cell_offset; - const int oldcid = cell_getid(oldcdim, ii, jj, kk) + old_bkg_cell_offset; - s->cells_top[cid].nodeID = oldnodeIDs[oldcid]; - } - } - } + const int cid = + cell_getid(s->cdim, i, j, k) + s->zoom_props->tl_cell_offset; + const int oldcid = + cell_getid(oldcdim, ii, jj, kk) + old_bkg_cell_offset; + s->cells_top[cid].nodeID = oldnodeIDs[oldcid]; + } + } + } - /* Check we have all nodeIDs present in the resample. */ - return check_complete(s, 1, s->e->nr_nodes); + /* Check we have all nodeIDs present in the resample. */ + return check_complete(s, 1, s->e->nr_nodes); } /* Vectorisation support */ @@ -159,28 +162,31 @@ int partition_space_to_space_zoom(double *oldh, double *oldcdim, double *oldzoom void pick_vector_zoom(struct space *s, int nregions, int *samplecells) { /* Get length of space and divide up. */ - int length = (s->cdim[0] * s->cdim[1] * s->cdim[2]) - + (s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * s->zoom_props->cdim[2]); + int length = (s->cdim[0] * s->cdim[1] * s->cdim[2]) + + (s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * + s->zoom_props->cdim[2]); if (nregions > length) { error("Too few cells (%d) for this number of regions (%d)", length, nregions); } - /* Set up variables */ - int step = (s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * s->zoom_props->cdim[2]) / nregions; + /* Set up variables */ + int step = (s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * + s->zoom_props->cdim[2]) / + nregions; int n = 0; int m = 0; int l = 0; /* Loop over zoom grid */ - for (int i = 0; i < s->zoom_props->cdim[0]; i++) { - for (int j = 0; j < s->zoom_props->cdim[1]; j++) { - for (int k = 0; k < s->zoom_props->cdim[2]; k++) { - if (n == 0 && l < nregions) { - samplecells[m++] = i; - samplecells[m++] = j; - samplecells[m++] = k; - l++; + for (int i = 0; i < s->zoom_props->cdim[0]; i++) { + for (int j = 0; j < s->zoom_props->cdim[1]; j++) { + for (int k = 0; k < s->zoom_props->cdim[2]; k++) { + if (n == 0 && l < nregions) { + samplecells[m++] = i; + samplecells[m++] = j; + samplecells[m++] = k; + l++; } n++; if (n == step) n = 0; @@ -188,19 +194,19 @@ void pick_vector_zoom(struct space *s, int nregions, int *samplecells) { } } - step = (s->cdim[0] * s->cdim[1] * s->cdim[2]) / nregions; + step = (s->cdim[0] * s->cdim[1] * s->cdim[2]) / nregions; n = 0; l = 0; /* Loop over natural grid */ - for (int i = 0; i < s->cdim[0]; i++) { - for (int j = 0; j < s->cdim[1]; j++) { - for (int k = 0; k < s->cdim[2]; k++) { - if (n == 0 && l < nregions) { - samplecells[m++] = i; - samplecells[m++] = j; - samplecells[m++] = k; - l++; + for (int i = 0; i < s->cdim[0]; i++) { + for (int j = 0; j < s->cdim[1]; j++) { + for (int k = 0; k < s->cdim[2]; k++) { + if (n == 0 && l < nregions) { + samplecells[m++] = i; + samplecells[m++] = j; + samplecells[m++] = k; + l++; } n++; if (n == step) n = 0; @@ -221,7 +227,7 @@ void pick_vector_zoom(struct space *s, int nregions, int *samplecells) { */ void split_vector_zoom(struct space *s, int nregions, int *samplecells) { - /* Define variables for selection */ + /* Define variables for selection */ int cid = 0; /* Loop over zoom cells*/ @@ -270,4 +276,3 @@ void split_vector_zoom(struct space *s, int nregions, int *samplecells) { } #endif #endif - diff --git a/src/zoom_region.c b/src/zoom_region.c index 87d97a86e2..22cb6beecf 100644 --- a/src/zoom_region.c +++ b/src/zoom_region.c @@ -22,16 +22,16 @@ ******************************************************************************/ /* Config parameters. */ -#include "../config.h" - -#include <float.h> +#include "zoom_region.h" +#include "../config.h" #include "cell.h" -#include "gravity_properties.h" #include "engine.h" +#include "gravity_properties.h" #include "proxy.h" #include "space.h" -#include "zoom_region.h" + +#include <float.h> /* MPI headers. */ #ifdef WITH_MPI @@ -39,15 +39,18 @@ #endif /** - * @brief Read parameter file for "ZoomRegion" properties, and initialize the zoom_region struct. + * @brief Read parameter file for "ZoomRegion" properties, and initialize the + * zoom_region struct. * * @param params Swift parameter structure. * @param s The space */ -void zoom_region_init(struct swift_params *params, struct space *s, int verbose) { +void zoom_region_init(struct swift_params *params, struct space *s, + int verbose) { #ifdef WITH_ZOOM_REGION /* Are we running with a zoom region? */ - s->with_zoom_region = parser_get_opt_param_int(params, "ZoomRegion:enable", 0); + s->with_zoom_region = + parser_get_opt_param_int(params, "ZoomRegion:enable", 0); /* If so... */ if (s->with_zoom_region) { @@ -59,18 +62,24 @@ void zoom_region_init(struct swift_params *params, struct space *s, int verbose) error("Error allocating memory for the zoom parameters."); /* Are refining the background cells? */ - s->zoom_props->refine_bkg = parser_get_opt_param_int(params, "ZoomRegion:enable_bkg_refinement", 1); + s->zoom_props->refine_bkg = + parser_get_opt_param_int(params, "ZoomRegion:enable_bkg_refinement", 1); /* Set the zoom cdim. */ - s->zoom_props->cdim[0] = parser_get_opt_param_int(params, "Scheduler:max_top_level_cells", - space_max_top_level_cells_default); - s->zoom_props->cdim[1] = parser_get_opt_param_int(params, "Scheduler:max_top_level_cells", - space_max_top_level_cells_default); - s->zoom_props->cdim[2] = parser_get_opt_param_int(params, "Scheduler:max_top_level_cells", - space_max_top_level_cells_default); - - /* Extract the zoom width boost factor (used to define the buffer around the zoom region). */ - s->zoom_props->zoom_boost_factor = parser_get_opt_param_float(params, "ZoomRegion:zoom_boost_factor", 1.1); + s->zoom_props->cdim[0] = + parser_get_opt_param_int(params, "Scheduler:max_top_level_cells", + space_max_top_level_cells_default); + s->zoom_props->cdim[1] = + parser_get_opt_param_int(params, "Scheduler:max_top_level_cells", + space_max_top_level_cells_default); + s->zoom_props->cdim[2] = + parser_get_opt_param_int(params, "Scheduler:max_top_level_cells", + space_max_top_level_cells_default); + + /* Extract the zoom width boost factor (used to define the buffer around the + * zoom region). */ + s->zoom_props->zoom_boost_factor = + parser_get_opt_param_float(params, "ZoomRegion:zoom_boost_factor", 1.1); /* Set the number of zoom cells in a natural cell. */ s->zoom_props->nr_zoom_per_bkg_cells = s->zoom_props->cdim[0]; @@ -78,7 +87,8 @@ void zoom_region_init(struct swift_params *params, struct space *s, int verbose) /* Initialise the number of wanders (unused if with_hydro == False)*/ s->zoom_props->nr_wanderers = 0; - /* Get an initial dimension for the zoom region and its geometric mid point. */ + /* Get an initial dimension for the zoom region and its geometric mid point. + */ double new_zoom_boundary[6] = {1e20, -1e20, 1e20, -1e20, 1e20, -1e20}; double midpoint[3] = {0.0, 0.0, 0.0}; const size_t nr_gparts = s->nr_gparts; @@ -88,9 +98,10 @@ void zoom_region_init(struct swift_params *params, struct space *s, int verbose) /* Get the shift from the ICs since this hasn't been applied yet. */ double shift[3] = {0.0, 0.0, 0.0}; parser_get_opt_param_double_array(params, "InitialConditions:shift", 3, - shift); + shift); - /* Find the min/max location in each dimension for each mask gravity particle, and their COM. */ + /* Find the min/max location in each dimension for each mask gravity + * particle, and their COM. */ for (size_t k = 0; k < nr_gparts; k++) { if (s->gparts[k].type != swift_type_dark_matter) continue; @@ -107,18 +118,12 @@ void zoom_region_init(struct swift_params *params, struct space *s, int verbose) } /* Ammend boundaries for this particle. */ - if (x < new_zoom_boundary[0]) - new_zoom_boundary[0] = x; - if (x > new_zoom_boundary[1]) - new_zoom_boundary[1] = x; - if (y < new_zoom_boundary[2]) - new_zoom_boundary[2] = y; - if (y > new_zoom_boundary[3]) - new_zoom_boundary[3] = y; - if (z < new_zoom_boundary[4]) - new_zoom_boundary[4] = z; - if (z > new_zoom_boundary[5]) - new_zoom_boundary[5] = z; + if (x < new_zoom_boundary[0]) new_zoom_boundary[0] = x; + if (x > new_zoom_boundary[1]) new_zoom_boundary[1] = x; + if (y < new_zoom_boundary[2]) new_zoom_boundary[2] = y; + if (y > new_zoom_boundary[3]) new_zoom_boundary[3] = y; + if (z < new_zoom_boundary[4]) new_zoom_boundary[4] = z; + if (z > new_zoom_boundary[5]) new_zoom_boundary[5] = z; /* Total up mass and position for COM. */ mtot += s->gparts[k].mass; @@ -158,21 +163,22 @@ void zoom_region_init(struct swift_params *params, struct space *s, int verbose) /* Store result. */ for (int ijk = 0; ijk < 3; ijk++) s->zoom_props->com[ijk] = com[ijk]; - if (verbose) - message("com: [%f %f %f]", com[0], com[1], com[2]); + if (verbose) message("com: [%f %f %f]", com[0], com[1], com[2]); if (verbose) - message("initial_dim: [%f %f %f] initial_zoom_boundary: [%f-%f %f-%f %f-%f]", - new_zoom_boundary[1] - new_zoom_boundary[0], - new_zoom_boundary[3] - new_zoom_boundary[2], - new_zoom_boundary[5] - new_zoom_boundary[4], - new_zoom_boundary[0], new_zoom_boundary[1], new_zoom_boundary[2], - new_zoom_boundary[3], new_zoom_boundary[4], new_zoom_boundary[5]); + message( + "initial_dim: [%f %f %f] initial_zoom_boundary: [%f-%f %f-%f %f-%f]", + new_zoom_boundary[1] - new_zoom_boundary[0], + new_zoom_boundary[3] - new_zoom_boundary[2], + new_zoom_boundary[5] - new_zoom_boundary[4], new_zoom_boundary[0], + new_zoom_boundary[1], new_zoom_boundary[2], new_zoom_boundary[3], + new_zoom_boundary[4], new_zoom_boundary[5]); /* Get the initial dimensions and midpoint. */ double ini_dim[3] = {0.0, 0.0, 0.0}; for (int ijk = 0; ijk < 3; ijk++) { - ini_dim[ijk] = (new_zoom_boundary[(ijk * 2) + 1] - new_zoom_boundary[ijk * 2]); + ini_dim[ijk] = + (new_zoom_boundary[(ijk * 2) + 1] - new_zoom_boundary[ijk * 2]); midpoint[ijk] = new_zoom_boundary[(ijk * 2) + 1] - (ini_dim[ijk] / 2); } @@ -181,44 +187,51 @@ void zoom_region_init(struct swift_params *params, struct space *s, int verbose) double shiftx = 0.; double shifty = 0.; double shiftz = 0.; - if ((ini_dim[0] > s->dim[0] / 2) || - (ini_dim[1] > s->dim[1] / 2) || + if ((ini_dim[0] > s->dim[0] / 2) || (ini_dim[1] > s->dim[1] / 2) || (ini_dim[2] > s->dim[2] / 2)) { if (ini_dim[0] > s->dim[0] / 2) shiftx = s->dim[0] / 2; if (ini_dim[1] > s->dim[1] / 2) shifty = s->dim[1] / 2; if (ini_dim[2] > s->dim[2] / 2) shiftz = s->dim[2] / 2; - error("Zoom region extends beyond the boundaries of the box. " - "Shift the ICs by [%f, %f, %f]", shiftx, shifty, shiftz); + error( + "Zoom region extends beyond the boundaries of the box. " + "Shift the ICs by [%f, %f, %f]", + shiftx, shifty, shiftz); } - /* Calculate the shift needed to place the mid point of the high res particles at the centre of the box. - * This shift is applied to the particles in space_init in space.c */ + /* Calculate the shift needed to place the mid point of the high res + * particles at the centre of the box. This shift is applied to the + * particles in space_init in space.c */ const double box_mid[3] = {s->dim[0] / 2, s->dim[1] / 2, s->dim[2] / 2}; for (int ijk = 0; ijk < 3; ijk++) { s->zoom_props->zoom_shift[ijk] = box_mid[ijk] - midpoint[ijk]; } if (verbose) { - message("box_mid = [%f %f %f] midpoint = [%f %f %f]", box_mid[0], box_mid[1], box_mid[2], - midpoint[0], midpoint[1], midpoint[2]); - message("Need to shift the box by [%e, %e, %e] to centre the zoom region", s->zoom_props->zoom_shift[0], - s->zoom_props->zoom_shift[1], s->zoom_props->zoom_shift[2]); + message("box_mid = [%f %f %f] midpoint = [%f %f %f]", box_mid[0], + box_mid[1], box_mid[2], midpoint[0], midpoint[1], midpoint[2]); + message("Need to shift the box by [%e, %e, %e] to centre the zoom region", + s->zoom_props->zoom_shift[0], s->zoom_props->zoom_shift[1], + s->zoom_props->zoom_shift[2]); } /* Let's shift the COM. * NOTE: boundaries are recalculated relative to box centre later. */ - for (int ijk = 0; ijk < 3; ijk++) s->zoom_props->com[ijk] += s->zoom_props->zoom_shift[ijk]; + for (int ijk = 0; ijk < 3; ijk++) + s->zoom_props->com[ijk] += s->zoom_props->zoom_shift[ijk]; - /* Compute maximum side length of the zoom region, we need zoom dim to be equal. */ - double max_dim = max3(ini_dim[0], ini_dim[1], ini_dim[2]) * s->zoom_props->zoom_boost_factor; + /* Compute maximum side length of the zoom region, we need zoom dim to be + * equal. */ + double max_dim = max3(ini_dim[0], ini_dim[1], ini_dim[2]) * + s->zoom_props->zoom_boost_factor; - /* This width has to divide the full parent box by an odd integer to ensure the two grids line up. - * NOTE: assumes box dimensions are equal! */ + /* This width has to divide the full parent box by an odd integer to ensure + * the two grids line up. NOTE: assumes box dimensions are equal! */ int nr_zoom_regions = (int)(s->dim[0] / max_dim); if (nr_zoom_regions % 2 == 0) nr_zoom_regions -= 1; max_dim = s->dim[0] / nr_zoom_regions; /* Do we want to refine the background cells? */ - if (s->zoom_props->refine_bkg && nr_zoom_regions >= s->zoom_props->cdim[0]) { + if (s->zoom_props->refine_bkg && + nr_zoom_regions >= s->zoom_props->cdim[0]) { /* Start with max_top_level_cells as a guess. */ nr_zoom_regions = s->zoom_props->cdim[0]; @@ -227,7 +240,9 @@ void zoom_region_init(struct swift_params *params, struct space *s, int verbose) if (nr_zoom_regions % 2 == 0) nr_zoom_regions -= 1; /* Compute the new boost factor and store it for reporting. */ - const float new_zoom_boost_factor = (s->dim[0] / nr_zoom_regions) / (max_dim / s->zoom_props->zoom_boost_factor); + const float new_zoom_boost_factor = + (s->dim[0] / nr_zoom_regions) / + (max_dim / s->zoom_props->zoom_boost_factor); if (verbose) message("Have increased zoom_boost_factor from %f to %f", @@ -242,17 +257,19 @@ void zoom_region_init(struct swift_params *params, struct space *s, int verbose) * The zoom region is already centred on the middle of the box */ for (int ijk = 0; ijk < 3; ijk++) { /* Set the new boundaries. */ - s->zoom_props->region_bounds[(ijk * 2)] = (s->dim[ijk] / 2) - (max_dim / 2); - s->zoom_props->region_bounds[(ijk * 2) + 1] = (s->dim[ijk] / 2) + (max_dim / 2); + s->zoom_props->region_bounds[(ijk * 2)] = + (s->dim[ijk] / 2) - (max_dim / 2); + s->zoom_props->region_bounds[(ijk * 2) + 1] = + (s->dim[ijk] / 2) + (max_dim / 2); /* Set the reigon dim. */ s->zoom_props->dim[ijk] = max_dim; } /* Set the minimum allowed zoom cell width. */ - const double zoom_dmax = max3(s->zoom_props->dim[0], s->zoom_props->dim[1], s->zoom_props->dim[2]); + const double zoom_dmax = max3(s->zoom_props->dim[0], s->zoom_props->dim[1], + s->zoom_props->dim[2]); s->zoom_props->cell_min = 0.99 * zoom_dmax / s->zoom_props->cdim[0]; - } #endif } @@ -260,9 +277,10 @@ void zoom_region_init(struct swift_params *params, struct space *s, int verbose) /** * @brief For a given particle location, what TL cell does it belong to? * - * Slightly more complicated in the zoom case, as there are now two embedded TL grids. - * - * First see if the particle is within the zoom bounds, then find its TL cell. + * Slightly more complicated in the zoom case, as there are now two embedded TL + * grids. + * + * First see if the particle is within the zoom bounds, then find its TL cell. * * @param s The space. * @param x, y, z Location of particle. @@ -278,12 +296,15 @@ int cell_getid_zoom(const struct space *s, const double x, const double y, /* Lets get some properties of the zoom region. */ const struct zoom_region_properties *zoom_props = s->zoom_props; - const int zoom_cdim[3] = {zoom_props->cdim[0], zoom_props->cdim[1], zoom_props->cdim[2]}; - const double zoom_iwidth[3] = {zoom_props->iwidth[0], zoom_props->iwidth[1], zoom_props->iwidth[2]}; + const int zoom_cdim[3] = {zoom_props->cdim[0], zoom_props->cdim[1], + zoom_props->cdim[2]}; + const double zoom_iwidth[3] = {zoom_props->iwidth[0], zoom_props->iwidth[1], + zoom_props->iwidth[2]}; const int bkg_cell_offset = zoom_props->tl_cell_offset; - const double zoom_region_bounds[6] = {zoom_props->region_bounds[0], zoom_props->region_bounds[1], - zoom_props->region_bounds[2], zoom_props->region_bounds[3], - zoom_props->region_bounds[4], zoom_props->region_bounds[5]}; + const double zoom_region_bounds[6] = { + zoom_props->region_bounds[0], zoom_props->region_bounds[1], + zoom_props->region_bounds[2], zoom_props->region_bounds[3], + zoom_props->region_bounds[4], zoom_props->region_bounds[5]}; /* Get the background cell ijk coordinates. */ const int bkg_i = x * iwidth[0]; @@ -306,7 +327,7 @@ int cell_getid_zoom(const struct space *s, const double x, const double y, error("cell_id out of range: %i (%f %f %f)", cell_id, x, y, z); #endif - /* If not then treat it like normal, and find the natural TL cell. */ + /* If not then treat it like normal, and find the natural TL cell. */ } else { cell_id = cell_getid(cdim, bkg_i, bkg_j, bkg_k) + bkg_cell_offset; @@ -314,7 +335,6 @@ int cell_getid_zoom(const struct space *s, const double x, const double y, if (cell_id < bkg_cell_offset || cell_id >= s->nr_cells) error("cell_id out of range: %i (%f %f %f)", cell_id, x, y, z); #endif - } return cell_id; @@ -326,54 +346,68 @@ int cell_getid_zoom(const struct space *s, const double x, const double y, #ifdef WITH_ZOOM_REGION #ifdef SWIFT_DEBUG_CHECKS /** -* @brief Run through all cells and ensure they have the correct cell type and width -* for their position in s->cells_top. -* -* @param s The space. -*/ + * @brief Run through all cells and ensure they have the correct cell type and + * width for their position in s->cells_top. + * + * @param s The space. + */ static void debug_cell_type(struct space *s) { - /* Get the cells array and cell properties */ - struct cell *cells = s->cells_top; - const int bkg_cell_offset = s->zoom_props->tl_cell_offset; - const double *zoom_width = s->zoom_props->width; - const double *width = s->width; - - /* Loop over all cells */ - for (int cid = 0; cid < s->nr_cells; cid++) { - - /* Check cell type */ - if (cid < bkg_cell_offset && cells[cid].tl_cell_type != 3) - error("Cell has the wrong cell type for it's array position (cid=%d, c->tl_cell_type=%d, " - "s->zoom_props->tl_cell_offset=%d)", cid, cells[cid].tl_cell_type, bkg_cell_offset); - if (cid >= bkg_cell_offset && cells[cid].tl_cell_type == 3) - error("Cell has the wrong cell type for it's array position (cid=%d, c->tl_cell_type=%d, " - "s->zoom_props->tl_cell_offset=%d)", cid, cells[cid].tl_cell_type, bkg_cell_offset); - - /* Check cell widths */ - for (int ijk = 0; ijk < 3; ijk++) { - if (cid < bkg_cell_offset && cells[cid].width[ijk] != zoom_width[ijk]) - error("Cell has the wrong cell width for it's array position (cid=%d, c->tl_cell_type=%d, " - "s->zoom_props->tl_cell_offset=%d, c->width=[%f %f %f], " - "s->zoom_props->width=[%f %f %f])", cid, cells[cid].tl_cell_type, bkg_cell_offset, - cells[cid].width[0], cells[cid].width[1], cells[cid].width[2], - s->zoom_props->width[0], s->zoom_props->width[1], s->zoom_props->width[2]); - if (cid >= bkg_cell_offset && cells[cid].width[ijk] != width[ijk]) - error("Cell has the wrong cell width for it's array position (cid=%d, c->tl_cell_type=%d, " - "s->zoom_props->tl_cell_offset=%d, c->width=[%f %f %f], " - "s->zoom_props->width=[%f %f %f])", cid, cells[cid].tl_cell_type, bkg_cell_offset, - cells[cid].width[0], cells[cid].width[1], cells[cid].width[2], - s->zoom_props->width[0], s->zoom_props->width[1], s->zoom_props->width[2]); - } - } + /* Get the cells array and cell properties */ + struct cell *cells = s->cells_top; + const int bkg_cell_offset = s->zoom_props->tl_cell_offset; + const double *zoom_width = s->zoom_props->width; + const double *width = s->width; + + /* Loop over all cells */ + for (int cid = 0; cid < s->nr_cells; cid++) { + + /* Check cell type */ + if (cid < bkg_cell_offset && cells[cid].tl_cell_type != 3) + error( + "Cell has the wrong cell type for it's array position (cid=%d, " + "c->tl_cell_type=%d, " + "s->zoom_props->tl_cell_offset=%d)", + cid, cells[cid].tl_cell_type, bkg_cell_offset); + if (cid >= bkg_cell_offset && cells[cid].tl_cell_type == 3) + error( + "Cell has the wrong cell type for it's array position (cid=%d, " + "c->tl_cell_type=%d, " + "s->zoom_props->tl_cell_offset=%d)", + cid, cells[cid].tl_cell_type, bkg_cell_offset); + + /* Check cell widths */ + for (int ijk = 0; ijk < 3; ijk++) { + if (cid < bkg_cell_offset && cells[cid].width[ijk] != zoom_width[ijk]) + error( + "Cell has the wrong cell width for it's array position (cid=%d, " + "c->tl_cell_type=%d, " + "s->zoom_props->tl_cell_offset=%d, c->width=[%f %f %f], " + "s->zoom_props->width=[%f %f %f])", + cid, cells[cid].tl_cell_type, bkg_cell_offset, cells[cid].width[0], + cells[cid].width[1], cells[cid].width[2], s->zoom_props->width[0], + s->zoom_props->width[1], s->zoom_props->width[2]); + if (cid >= bkg_cell_offset && cells[cid].width[ijk] != width[ijk]) + error( + "Cell has the wrong cell width for it's array position (cid=%d, " + "c->tl_cell_type=%d, " + "s->zoom_props->tl_cell_offset=%d, c->width=[%f %f %f], " + "s->zoom_props->width=[%f %f %f])", + cid, cells[cid].tl_cell_type, bkg_cell_offset, cells[cid].width[0], + cells[cid].width[1], cells[cid].width[2], s->zoom_props->width[0], + s->zoom_props->width[1], s->zoom_props->width[2]); + } + } } #endif #endif /** - * @brief Compute the extent/bounds of the zoom region using the high-res DM particles. + * @brief Compute the extent/bounds of the zoom region using the high-res DM + * particles. * - * The min/max [x,y,z] for each particle is found, and the CoM of these particles is computed. + * The min/max [x,y,z] for each particle is found, and the CoM of these + * particles is computed. * * @param s The space. * @param verbose Are we talking? @@ -382,22 +416,23 @@ void construct_zoom_region(struct space *s, int verbose) { #ifdef WITH_ZOOM_REGION /* Get the width of the zoom region, zoom dims are equal. */ const double zoom_dim = s->zoom_props->dim[0]; - + /* Let's set what we know about the zoom region. */ for (int ijk = 0; ijk < 3; ijk++) { - s->zoom_props->width[ijk] = zoom_dim / s->zoom_props->cdim[ijk]; - s->zoom_props->iwidth[ijk] = 1 / s->zoom_props->width[ijk]; - s->zoom_props->dim[ijk] = zoom_dim; + s->zoom_props->width[ijk] = zoom_dim / s->zoom_props->cdim[ijk]; + s->zoom_props->iwidth[ijk] = 1 / s->zoom_props->width[ijk]; + s->zoom_props->dim[ijk] = zoom_dim; } /* Overwrite the minimum allowed zoom cell width. */ s->zoom_props->cell_min = 0.99 * zoom_dim / s->zoom_props->cdim[0]; - + /* Now we can define the background grid and zoom region's background ijk. */ for (int ijk = 0; ijk < 3; ijk++) { s->width[ijk] = s->zoom_props->dim[ijk]; s->iwidth[ijk] = 1.0 / s->width[ijk]; - s->cdim[ijk] = (int)floor((s->dim[ijk] + 0.1 * s->width[ijk]) * s->iwidth[ijk]); + s->cdim[ijk] = + (int)floor((s->dim[ijk] + 0.1 * s->width[ijk]) * s->iwidth[ijk]); s->zoom_props->zoom_cell_ijk[ijk] = (int)floor(s->cdim[ijk] / 2); } @@ -407,48 +442,60 @@ void construct_zoom_region(struct space *s, int verbose) { s->cell_min = 0.99 * dmax / nr_zoom_regions; /* Check we have enough cells for periodicity. */ - if (s->periodic && (s->cdim[0] < 3 || s->cdim[1] < 3 || s->cdim[2] < 3)) - error( - "Must have at least 3 cells in each spatial dimension when periodicity " - "is switched on.\nThis error is often caused by any of the " - "followings:\n" - " - too few particles to generate a sensible grid,\n" - " - the initial value of 'Scheduler:max_top_level_cells' is too " - "small,\n" - " - the (minimal) time-step is too large leading to particles with " - "predicted smoothing lengths too large for the box size,\n" - " - particles with velocities so large that they move by more than two " - "box sizes per time-step.\n"); - - /* Store cell number information. */ - s->zoom_props->tl_cell_offset = s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * s->zoom_props->cdim[2]; - s->zoom_props->nr_zoom_cells = s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * s->zoom_props->cdim[2]; + if (s->periodic && (s->cdim[0] < 3 || s->cdim[1] < 3 || s->cdim[2] < 3)) + error( + "Must have at least 3 cells in each spatial dimension when periodicity " + "is switched on.\nThis error is often caused by any of the " + "followings:\n" + " - too few particles to generate a sensible grid,\n" + " - the initial value of 'Scheduler:max_top_level_cells' is too " + "small,\n" + " - the (minimal) time-step is too large leading to particles with " + "predicted smoothing lengths too large for the box size,\n" + " - particles with velocities so large that they move by more than two " + "box sizes per time-step.\n"); + + /* Store cell number information. */ + s->zoom_props->tl_cell_offset = + s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * s->zoom_props->cdim[2]; + s->zoom_props->nr_zoom_cells = + s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * s->zoom_props->cdim[2]; s->zoom_props->nr_bkg_cells = s->cdim[0] * s->cdim[1] * s->cdim[2]; /* Lets report what we have constructed. */ if (verbose) { - message("set cell dimensions to zoom_cdim=[%d %d %d] background_cdim=[%d %d %d]", s->zoom_props->cdim[0], - s->zoom_props->cdim[1], s->zoom_props->cdim[2], s->cdim[0], - s->cdim[1], s->cdim[2]); - message("zoom region located in cell [%d %d %d]", s->zoom_props->zoom_cell_ijk[0], - s->zoom_props->zoom_cell_ijk[1], s->zoom_props->zoom_cell_ijk[2]); - message("nr_zoom_cells: %d nr_bkg_cells: %d tl_cell_offset: %d", s->zoom_props->nr_zoom_cells, - s->zoom_props->nr_bkg_cells, s->zoom_props->tl_cell_offset); - message("zoom_boundary: [%f-%f %f-%f %f-%f]", + message( + "set cell dimensions to zoom_cdim=[%d %d %d] background_cdim=[%d %d " + "%d]", + s->zoom_props->cdim[0], s->zoom_props->cdim[1], s->zoom_props->cdim[2], + s->cdim[0], s->cdim[1], s->cdim[2]); + message("zoom region located in cell [%d %d %d]", + s->zoom_props->zoom_cell_ijk[0], s->zoom_props->zoom_cell_ijk[1], + s->zoom_props->zoom_cell_ijk[2]); + message("nr_zoom_cells: %d nr_bkg_cells: %d tl_cell_offset: %d", + s->zoom_props->nr_zoom_cells, s->zoom_props->nr_bkg_cells, + s->zoom_props->tl_cell_offset); + message("zoom_boundary: [%f-%f %f-%f %f-%f]", s->zoom_props->region_bounds[0], s->zoom_props->region_bounds[1], s->zoom_props->region_bounds[2], s->zoom_props->region_bounds[3], s->zoom_props->region_bounds[4], s->zoom_props->region_bounds[5]); - message("dim: [%f %f %f] tl_cell_width: [%f %f %f] zoom_cell_width: [%f %f %f]", - s->zoom_props->dim[0], s->zoom_props->dim[1], s->zoom_props->dim[2], - s->width[0], s->width[1], s->width[2], - s->zoom_props->width[0], s->zoom_props->width[1], s->zoom_props->width[2]); - message("nr_tl_cells_in_zoom_region: [%d %d %d] nr_zoom_cells_in_tl_cell: [%d %d %d]", - (int)floor((s->zoom_props->dim[0] + 0.5 * s->width[0]) * s->iwidth[0]), - (int)floor((s->zoom_props->dim[1] + 0.5 * s->width[1]) * s->iwidth[1]), - (int)floor((s->zoom_props->dim[2] + 0.5 * s->width[2]) * s->iwidth[2]), - (int)floor((s->width[0] + 0.5 * s->zoom_props->width[0]) * s->zoom_props->iwidth[0]), - (int)floor((s->width[1] + 0.5 * s->zoom_props->width[1]) * s->zoom_props->iwidth[1]), - (int)floor((s->width[2] + 0.5 * s->zoom_props->width[2]) * s->zoom_props->iwidth[2])); + message( + "dim: [%f %f %f] tl_cell_width: [%f %f %f] zoom_cell_width: [%f %f %f]", + s->zoom_props->dim[0], s->zoom_props->dim[1], s->zoom_props->dim[2], + s->width[0], s->width[1], s->width[2], s->zoom_props->width[0], + s->zoom_props->width[1], s->zoom_props->width[2]); + message( + "nr_tl_cells_in_zoom_region: [%d %d %d] nr_zoom_cells_in_tl_cell: [%d " + "%d %d]", + (int)floor((s->zoom_props->dim[0] + 0.5 * s->width[0]) * s->iwidth[0]), + (int)floor((s->zoom_props->dim[1] + 0.5 * s->width[1]) * s->iwidth[1]), + (int)floor((s->zoom_props->dim[2] + 0.5 * s->width[2]) * s->iwidth[2]), + (int)floor((s->width[0] + 0.5 * s->zoom_props->width[0]) * + s->zoom_props->iwidth[0]), + (int)floor((s->width[1] + 0.5 * s->zoom_props->width[1]) * + s->zoom_props->iwidth[1]), + (int)floor((s->width[2] + 0.5 * s->zoom_props->width[2]) * + s->zoom_props->iwidth[2])); } #endif @@ -458,27 +505,32 @@ void construct_zoom_region(struct space *s, int verbose) { * @brief Build the TL cells, with a zoom region. * * This replaces the loop in space_regrid when running with a zoom region. - * - * Construct an additional set of TL "zoom" cells embedded within the TL cell structure - * with the dimensions of each cell structure being the same (with differing widths). * - * Therefore the new TL cell structure is 2*cdim**3, with the "natural" TL cells occupying the - * first half of the TL cell list, and the "zoom" TL cells ocupying the second half. + * Construct an additional set of TL "zoom" cells embedded within the TL cell + * structure with the dimensions of each cell structure being the same (with + * differing widths). + * + * Therefore the new TL cell structure is 2*cdim**3, with the "natural" TL cells + * occupying the first half of the TL cell list, and the "zoom" TL cells + * ocupying the second half. * * @param s The space. * @param verbose Are we talking? */ -void construct_tl_cells_with_zoom_region(struct space *s, const int *cdim, const float dmin, - const integertime_t ti_current, struct gravity_props *gravity_properties, int verbose) { +void construct_tl_cells_with_zoom_region( + struct space *s, const int *cdim, const float dmin, + const integertime_t ti_current, struct gravity_props *gravity_properties, + int verbose) { #ifdef WITH_ZOOM_REGION /* Get some zoom region properties */ const float dmin_zoom = min3(s->zoom_props->width[0], s->zoom_props->width[1], s->zoom_props->width[2]); const int bkg_cell_offset = s->zoom_props->tl_cell_offset; - const double zoom_region_bounds[6] = {s->zoom_props->region_bounds[0], s->zoom_props->region_bounds[1], - s->zoom_props->region_bounds[2], s->zoom_props->region_bounds[3], - s->zoom_props->region_bounds[4], s->zoom_props->region_bounds[5]}; + const double zoom_region_bounds[6] = { + s->zoom_props->region_bounds[0], s->zoom_props->region_bounds[1], + s->zoom_props->region_bounds[2], s->zoom_props->region_bounds[3], + s->zoom_props->region_bounds[4], s->zoom_props->region_bounds[5]}; struct cell *restrict c; @@ -488,46 +540,45 @@ void construct_tl_cells_with_zoom_region(struct space *s, const int *cdim, const for (int k = 0; k < s->zoom_props->cdim[2]; k++) { const size_t cid = cell_getid(s->zoom_props->cdim, i, j, k); - /* Create the zoom cell and it's multipoles */ - c = &s->cells_top[cid]; - c->loc[0] = i * s->zoom_props->width[0] + zoom_region_bounds[0]; - c->loc[1] = j * s->zoom_props->width[1] + zoom_region_bounds[2]; - c->loc[2] = k * s->zoom_props->width[2] + zoom_region_bounds[4]; - c->parent_bkg_cid = cell_getid(s->cdim, - s->zoom_props->zoom_cell_ijk[0], - s->zoom_props->zoom_cell_ijk[1], - s->zoom_props->zoom_cell_ijk[2]) + bkg_cell_offset; - c->width[0] = s->zoom_props->width[0]; - c->width[1] = s->zoom_props->width[1]; - c->width[2] = s->zoom_props->width[2]; - if (s->with_self_gravity) - c->grav.multipole = &s->multipoles_top[cid]; - c->tl_cell_type = zoom_tl_cell; - c->dmin = dmin_zoom; - c->nr_zoom_per_bkg_cells = s->zoom_props->nr_zoom_per_bkg_cells; - c->depth = 0; - c->split = 0; - c->hydro.count = 0; - c->grav.count = 0; - c->stars.count = 0; - c->sinks.count = 0; - c->top = c; - c->super = c; - c->hydro.super = c; - c->grav.super = c; - c->hydro.ti_old_part = ti_current; - c->grav.ti_old_part = ti_current; - c->stars.ti_old_part = ti_current; - c->sinks.ti_old_part = ti_current; - c->black_holes.ti_old_part = ti_current; - c->grav.ti_old_multipole = ti_current; + /* Create the zoom cell and it's multipoles */ + c = &s->cells_top[cid]; + c->loc[0] = i * s->zoom_props->width[0] + zoom_region_bounds[0]; + c->loc[1] = j * s->zoom_props->width[1] + zoom_region_bounds[2]; + c->loc[2] = k * s->zoom_props->width[2] + zoom_region_bounds[4]; + c->parent_bkg_cid = cell_getid(s->cdim, s->zoom_props->zoom_cell_ijk[0], + s->zoom_props->zoom_cell_ijk[1], + s->zoom_props->zoom_cell_ijk[2]) + + bkg_cell_offset; + c->width[0] = s->zoom_props->width[0]; + c->width[1] = s->zoom_props->width[1]; + c->width[2] = s->zoom_props->width[2]; + if (s->with_self_gravity) c->grav.multipole = &s->multipoles_top[cid]; + c->tl_cell_type = zoom_tl_cell; + c->dmin = dmin_zoom; + c->nr_zoom_per_bkg_cells = s->zoom_props->nr_zoom_per_bkg_cells; + c->depth = 0; + c->split = 0; + c->hydro.count = 0; + c->grav.count = 0; + c->stars.count = 0; + c->sinks.count = 0; + c->top = c; + c->super = c; + c->hydro.super = c; + c->grav.super = c; + c->hydro.ti_old_part = ti_current; + c->grav.ti_old_part = ti_current; + c->stars.ti_old_part = ti_current; + c->sinks.ti_old_part = ti_current; + c->black_holes.ti_old_part = ti_current; + c->grav.ti_old_multipole = ti_current; #ifdef WITH_MPI - c->mpi.tag = -1; - c->mpi.recv = NULL; - c->mpi.send = NULL; + c->mpi.tag = -1; + c->mpi.recv = NULL; + c->mpi.send = NULL; #endif #if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_CELL_GRAPH) - cell_assign_top_level_cell_index(c, s); + cell_assign_top_level_cell_index(c, s); #endif } } @@ -538,52 +589,50 @@ void construct_tl_cells_with_zoom_region(struct space *s, const int *cdim, const for (int k = 0; k < cdim[2]; k++) { const size_t cid = cell_getid(cdim, i, j, k); - /* Natural top level cells. */ - c = &s->cells_top[cid + bkg_cell_offset]; - c->loc[0] = i * s->width[0]; - c->loc[1] = j * s->width[1]; - c->loc[2] = k * s->width[2]; - c->width[0] = s->width[0]; - c->width[1] = s->width[1]; - c->width[2] = s->width[2]; - c->dmin = dmin; - c->parent_bkg_cid = cid + bkg_cell_offset; - c->nr_zoom_per_bkg_cells = s->zoom_props->nr_zoom_per_bkg_cells; - if (s->with_self_gravity) - c->grav.multipole = &s->multipoles_top[cid + bkg_cell_offset]; - c->tl_cell_type = tl_cell; - c->depth = 0; - c->split = 0; - c->hydro.count = 0; - c->grav.count = 0; - c->stars.count = 0; - c->sinks.count = 0; - c->top = c; - c->super = c; - c->hydro.super = c; - c->grav.super = c; - c->hydro.ti_old_part = ti_current; - c->grav.ti_old_part = ti_current; - c->stars.ti_old_part = ti_current; - c->sinks.ti_old_part = ti_current; - c->black_holes.ti_old_part = ti_current; - c->grav.ti_old_multipole = ti_current; + /* Natural top level cells. */ + c = &s->cells_top[cid + bkg_cell_offset]; + c->loc[0] = i * s->width[0]; + c->loc[1] = j * s->width[1]; + c->loc[2] = k * s->width[2]; + c->width[0] = s->width[0]; + c->width[1] = s->width[1]; + c->width[2] = s->width[2]; + c->dmin = dmin; + c->parent_bkg_cid = cid + bkg_cell_offset; + c->nr_zoom_per_bkg_cells = s->zoom_props->nr_zoom_per_bkg_cells; + if (s->with_self_gravity) + c->grav.multipole = &s->multipoles_top[cid + bkg_cell_offset]; + c->tl_cell_type = tl_cell; + c->depth = 0; + c->split = 0; + c->hydro.count = 0; + c->grav.count = 0; + c->stars.count = 0; + c->sinks.count = 0; + c->top = c; + c->super = c; + c->hydro.super = c; + c->grav.super = c; + c->hydro.ti_old_part = ti_current; + c->grav.ti_old_part = ti_current; + c->stars.ti_old_part = ti_current; + c->sinks.ti_old_part = ti_current; + c->black_holes.ti_old_part = ti_current; + c->grav.ti_old_multipole = ti_current; #ifdef WITH_MPI - c->mpi.tag = -1; - c->mpi.recv = NULL; - c->mpi.send = NULL; + c->mpi.tag = -1; + c->mpi.recv = NULL; + c->mpi.send = NULL; #endif #if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_CELL_GRAPH) - cell_assign_top_level_cell_index(c, s); + cell_assign_top_level_cell_index(c, s); #endif - } } } /* We need to label the zoom region's background cell as void. */ - const size_t void_cid = cell_getid(cdim, - s->zoom_props->zoom_cell_ijk[0], + const size_t void_cid = cell_getid(cdim, s->zoom_props->zoom_cell_ijk[0], s->zoom_props->zoom_cell_ijk[1], s->zoom_props->zoom_cell_ijk[2]); c = &s->cells_top[void_cid + bkg_cell_offset]; @@ -595,7 +644,8 @@ void construct_tl_cells_with_zoom_region(struct space *s, const int *cdim, const #endif /* Now find what cells neighbour the zoom region. */ - if (s->with_zoom_region) find_neighbouring_cells(s, gravity_properties, verbose); + if (s->with_zoom_region) + find_neighbouring_cells(s, gravity_properties, verbose); #endif } @@ -603,33 +653,36 @@ void construct_tl_cells_with_zoom_region(struct space *s, const int *cdim, const /** * @brief Find what TL cells surround the zoom region. * - * When interacting "natural" TL cells and "zoom" TL cells, it helps to know what natural TL - * cells surround the zoom region. These cells then get tagged as "tl_cell_neighbour". + * When interacting "natural" TL cells and "zoom" TL cells, it helps to know + * what natural TL cells surround the zoom region. These cells then get tagged + * as "tl_cell_neighbour". * * @param s The space. * @param verbose Are we talking? */ -void find_neighbouring_cells(struct space *s, struct gravity_props *gravity_properties, const int verbose) { +void find_neighbouring_cells(struct space *s, + struct gravity_props *gravity_properties, + const int verbose) { #ifdef WITH_ZOOM_REGION const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; const int periodic = s->periodic; struct cell *cells = s->cells_top; /* Some info about the zoom domain */ - const int bkg_cell_offset = s->zoom_props->tl_cell_offset; + const int bkg_cell_offset = s->zoom_props->tl_cell_offset; /* Get some info about the physics */ - const double theta_crit_inv = 1. / gravity_properties->theta_crit; + const double theta_crit_inv = 1. / gravity_properties->theta_crit; - /* Maximal distance from shifted CoM to any corner */ - const double distance = 2. * cells[bkg_cell_offset].width[0] * theta_crit_inv; + /* Maximal distance from shifted CoM to any corner */ + const double distance = 2. * cells[bkg_cell_offset].width[0] * theta_crit_inv; - /* Compute how many cells away we need to walk */ - const int delta_cells = (int)(distance / cells[bkg_cell_offset].dmin) + 1; + /* Compute how many cells away we need to walk */ + const int delta_cells = (int)(distance / cells[bkg_cell_offset].dmin) + 1; - /* Turn this into upper and lower bounds for loops */ - const int delta_m = delta_cells; - const int delta_p = delta_cells; + /* Turn this into upper and lower bounds for loops */ + const int delta_m = delta_cells; + const int delta_p = delta_cells; int neighbour_count = 0; int void_count = 0; @@ -637,7 +690,8 @@ void find_neighbouring_cells(struct space *s, struct gravity_props *gravity_prop /* Let's be verbose about this choice */ if (verbose) message( - "Looking for neighbouring natural cells up to %d natural top-level cells away from the zoom region (delta_m=%d " + "Looking for neighbouring natural cells up to %d natural top-level " + "cells away from the zoom region (delta_m=%d " "delta_p=%d)", delta_cells, delta_m, delta_p); @@ -652,8 +706,8 @@ void find_neighbouring_cells(struct space *s, struct gravity_props *gravity_prop #ifdef SWIFT_DEBUG_CHECKS /* Ensure all background cells are actually background cells */ - if (cells[cid].width[0] == s->zoom_props->width[0]) - error("A zoom cell has been given a natural cell label!"); + if (cells[cid].width[0] == s->zoom_props->width[0]) + error("A zoom cell has been given a natural cell label!"); #endif /* Only interested in cells hosting zoom top level cells. */ @@ -692,8 +746,8 @@ void find_neighbouring_cells(struct space *s, struct gravity_props *gravity_prop } if (verbose) { - message("%i cells neighbouring the zoom region", neighbour_count); - message("%i void cells in the zoom region", void_count); + message("%i cells neighbouring the zoom region", neighbour_count); + message("%i void cells in the zoom region", void_count); } #endif } @@ -716,16 +770,20 @@ double cell_min_dist2_diff_size(const struct cell *restrict ci, if (ci->width[2] == cj->width[2]) error("z cells of same size!"); #endif - const double cix = ci->loc[0] + ci->width[0]/2.; - const double ciy = ci->loc[1] + ci->width[1]/2.; - const double ciz = ci->loc[2] + ci->width[2]/2.; + const double cix = ci->loc[0] + ci->width[0] / 2.; + const double ciy = ci->loc[1] + ci->width[1] / 2.; + const double ciz = ci->loc[2] + ci->width[2] / 2.; - const double cjx = cj->loc[0] + cj->width[0]/2.; - const double cjy = cj->loc[1] + cj->width[1]/2.; - const double cjz = cj->loc[2] + cj->width[2]/2.; + const double cjx = cj->loc[0] + cj->width[0] / 2.; + const double cjy = cj->loc[1] + cj->width[1] / 2.; + const double cjz = cj->loc[2] + cj->width[2] / 2.; - const double diag_ci2 = ci->width[0] * ci->width[0] + ci->width[1] * ci->width[1] + ci->width[2] * ci->width[2]; - const double diag_cj2 = cj->width[0] * cj->width[0] + cj->width[1] * cj->width[1] + cj->width[2] * cj->width[2]; + const double diag_ci2 = ci->width[0] * ci->width[0] + + ci->width[1] * ci->width[1] + + ci->width[2] * ci->width[2]; + const double diag_cj2 = cj->width[0] * cj->width[0] + + cj->width[1] * cj->width[1] + + cj->width[2] * cj->width[2]; /* Get the distance between the cells */ double dx = cix - cjx; @@ -741,9 +799,9 @@ double cell_min_dist2_diff_size(const struct cell *restrict ci, const double r2 = dx * dx + dy * dy + dz * dz; /* Minimal distance between any 2 particles in the two cells */ - const double dist2 = r2 - (diag_ci2/2. + diag_cj2/2.); + const double dist2 = r2 - (diag_ci2 / 2. + diag_cj2 / 2.); - return dist2; + return dist2; #else return 0; #endif @@ -752,7 +810,8 @@ double cell_min_dist2_diff_size(const struct cell *restrict ci, /** * @brief Minimum distance between two TL cells. * - * Generic wrapper, don't know if the TL cells are the same size or not at time of calling. + * Generic wrapper, don't know if the TL cells are the same size or not at time + * of calling. * * @param ci, cj The two TL cells. * @param periodic Account for periodicity? @@ -767,11 +826,11 @@ double cell_min_dist2(const struct cell *restrict ci, /* Two natural TL cells. */ if (ci->tl_cell_type <= 2 && cj->tl_cell_type <= 2) { dist2 = cell_min_dist2_same_size(ci, cj, periodic, dim); - /* Two zoom TL cells. */ + /* Two zoom TL cells. */ } else if (ci->tl_cell_type == zoom_tl_cell && cj->tl_cell_type == zoom_tl_cell) { dist2 = cell_min_dist2_same_size(ci, cj, periodic, dim); - /* A mix of natural and zoom TL cells. */ + /* A mix of natural and zoom TL cells. */ } else { dist2 = cell_min_dist2_diff_size(ci, cj, periodic, dim); } @@ -797,7 +856,7 @@ void engine_makeproxies_natural_cells(struct engine *e) { /* Handle on the cells and proxies */ struct cell *cells = s->cells_top; struct proxy *proxies = e->proxies; - + /* Some info about the zoom domain */ const int bkg_cell_offset = s->zoom_props->tl_cell_offset; @@ -805,7 +864,8 @@ void engine_makeproxies_natural_cells(struct engine *e) { const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]}; const int periodic = s->periodic; - const double cell_width[3] = {cells[bkg_cell_offset].width[0], cells[bkg_cell_offset].width[1], + const double cell_width[3] = {cells[bkg_cell_offset].width[0], + cells[bkg_cell_offset].width[1], cells[bkg_cell_offset].width[2]}; /* Get some info about the physics */ @@ -1056,7 +1116,8 @@ void engine_makeproxies_zoom_cells(struct engine *e) { struct proxy *proxies = e->proxies; /* Some info about the domain */ - const int cdim[3] = {s->zoom_props->cdim[0], s->zoom_props->cdim[1], s->zoom_props->cdim[2]}; + const int cdim[3] = {s->zoom_props->cdim[0], s->zoom_props->cdim[1], + s->zoom_props->cdim[2]}; const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]}; const int periodic = 0; const double cell_width[3] = {cells[0].width[0], cells[0].width[1], @@ -1152,7 +1213,8 @@ void engine_makeproxies_zoom_cells(struct engine *e) { if (with_hydro) { /* Check for direct neighbours without periodic BC */ - if ((abs(i - iii) <= 1) && (abs(j - jjj) <= 1) && (abs(k - kkk) <= 1)) + if ((abs(i - iii) <= 1) && (abs(j - jjj) <= 1) && + (abs(k - kkk) <= 1)) proxy_type |= (int)proxy_cell_type_hydro; } @@ -1163,7 +1225,8 @@ void engine_makeproxies_zoom_cells(struct engine *e) { some further out if the opening angle demands it */ /* Check for direct neighbours without periodic BC */ - if ((abs(i - iii) <= 1) && (abs(j - jjj) <= 1) && (abs(k - kkk) <= 1)) { + if ((abs(i - iii) <= 1) && (abs(j - jjj) <= 1) && + (abs(k - kkk) <= 1)) { proxy_type |= (int)proxy_cell_type_gravity; @@ -1282,7 +1345,7 @@ void engine_makeproxies_zoom_cells(struct engine *e) { /** * @brief Create and fill the proxies for relations between cell grids. - * + * * This is done "lazily" by just making proxies for all neighbour cells. * * @param e The #engine. @@ -1296,7 +1359,7 @@ void engine_makeproxies_between_grids(struct engine *e) { /* Handle on the cells and proxies */ struct cell *cells = s->cells_top; struct proxy *proxies = e->proxies; - + /* Some info about the zoom domain */ const int bkg_cell_offset = s->zoom_props->tl_cell_offset; @@ -1304,7 +1367,8 @@ void engine_makeproxies_between_grids(struct engine *e) { const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]}; const int periodic = s->periodic; - const double cell_width[3] = {cells[bkg_cell_offset].width[0], cells[bkg_cell_offset].width[1], + const double cell_width[3] = {cells[bkg_cell_offset].width[0], + cells[bkg_cell_offset].width[1], cells[bkg_cell_offset].width[2]}; /* Get some info about the physics */ @@ -1325,31 +1389,29 @@ void engine_makeproxies_between_grids(struct engine *e) { /* Loop over each zoom cell in the space. */ for (int cid = 0; cid < bkg_cell_offset; cid++) { - + /* Integer indices of this cell in the natural parent */ const int natural_tl_cid = cells[cid].parent_bkg_cid - bkg_cell_offset; const int i = natural_tl_cid / (cdim[1] * cdim[2]); const int j = (natural_tl_cid / cdim[2]) % cdim[1]; const int k = natural_tl_cid % cdim[2]; - /* Loop over every cell in the natural grid */ - for (int cjd = bkg_cell_offset; cjd < s->nr_cells; cjd++) { + /* Loop over every cell in the natural grid */ + for (int cjd = bkg_cell_offset; cjd < s->nr_cells; cjd++) { - /* We only want to consider background cells if they are neighbours */ - if (cells[cjd].tl_cell_type != tl_cell_neighbour) continue; + /* We only want to consider background cells if they are neighbours */ + if (cells[cjd].tl_cell_type != tl_cell_neighbour) continue; /* Early abort (both same node) */ - if (cells[cid].nodeID == nodeID && cells[cjd].nodeID == nodeID) - continue; + if (cells[cid].nodeID == nodeID && cells[cjd].nodeID == nodeID) continue; /* Early abort (both foreign node) */ - if (cells[cid].nodeID != nodeID && cells[cjd].nodeID != nodeID) - continue; - + if (cells[cid].nodeID != nodeID && cells[cjd].nodeID != nodeID) continue; + /* Integer indices of the cell in the top-level grid */ - const int ii = (cjd - bkg_cell_offset) / (cdim[1] * cdim[2]); - const int jj = ((cjd - bkg_cell_offset) / cdim[2]) % cdim[1]; - const int kk = (cjd - bkg_cell_offset) % cdim[2]; + const int ii = (cjd - bkg_cell_offset) / (cdim[1] * cdim[2]); + const int jj = ((cjd - bkg_cell_offset) / cdim[2]) % cdim[1]; + const int kk = (cjd - bkg_cell_offset) % cdim[2]; int proxy_type = 0; @@ -1398,16 +1460,15 @@ void engine_makeproxies_between_grids(struct engine *e) { cannot rely on just an M2L calculation. */ /* Minimal distance between any two points in the cells */ - const double min_dist_CoM2 = cell_min_dist2_diff_size( - &cells[cid], &cells[cjd], periodic, dim); + const double min_dist_CoM2 = + cell_min_dist2_diff_size(&cells[cid], &cells[cjd], periodic, dim); /* Are we beyond the distance where the truncated forces are 0 * but not too far such that M2L can be used? */ if (periodic) { if ((min_dist_CoM2 < max_mesh_dist2) && - !(4. * r_max * r_max < - theta_crit * theta_crit * min_dist_CoM2)) + !(4. * r_max * r_max < theta_crit * theta_crit * min_dist_CoM2)) proxy_type |= (int)proxy_cell_type_gravity; } else { @@ -1433,8 +1494,7 @@ void engine_makeproxies_between_grids(struct engine *e) { error("Maximum number of proxies exceeded."); /* Ok, start a new proxy for this pair of nodes */ - proxy_init(&proxies[e->nr_proxies], e->nodeID, - cells[cjd].nodeID); + proxy_init(&proxies[e->nr_proxies], e->nodeID, cells[cjd].nodeID); /* Store the information */ e->proxy_ind[cells[cjd].nodeID] = e->nr_proxies; @@ -1467,8 +1527,7 @@ void engine_makeproxies_between_grids(struct engine *e) { error("Maximum number of proxies exceeded."); /* Ok, start a new proxy for this pair of nodes */ - proxy_init(&proxies[e->nr_proxies], e->nodeID, - cells[cid].nodeID); + proxy_init(&proxies[e->nr_proxies], e->nodeID, cells[cid].nodeID); /* Store the information */ e->proxy_ind[cells[cid].nodeID] = e->nr_proxies; @@ -1505,17 +1564,17 @@ void engine_makeproxies_between_grids(struct engine *e) { void engine_makeproxies_with_zoom_region(struct engine *e) { #ifdef WITH_MPI - /* Let's time this */ + /* Let's time this */ const ticks tic = getticks(); - /* Prepare the proxies and the proxy index. */ + /* Prepare the proxies and the proxy index. */ if (e->proxy_ind == NULL) if ((e->proxy_ind = (int *)malloc(sizeof(int) * e->nr_nodes)) == NULL) error("Failed to allocate proxy index."); for (int k = 0; k < e->nr_nodes; k++) e->proxy_ind[k] = -1; e->nr_proxies = 0; - engine_makeproxies_zoom_cells(e); + engine_makeproxies_zoom_cells(e); engine_makeproxies_natural_cells(e); engine_makeproxies_between_grids(e); @@ -1525,7 +1584,7 @@ void engine_makeproxies_with_zoom_region(struct engine *e) { clocks_getunit()); #else - error("SWIFT was not compiled with MPI support."); + error("SWIFT was not compiled with MPI support."); #endif } @@ -1537,122 +1596,131 @@ void engine_makeproxies_with_zoom_region(struct engine *e) { * - All pairs within range according to the multipole acceptance * criterion get a pair task. */ -void engine_make_self_gravity_tasks_mapper_natural_cells(void *map_data, int num_elements, - void *extra_data) { - - struct engine *e = (struct engine *)extra_data; - struct space *s = e->s; - struct scheduler *sched = &e->sched; - const int nodeID = e->nodeID; - const int periodic = s->periodic; - const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]}; - const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; - struct cell *cells = s->cells_top; - const double theta_crit = e->gravity_properties->theta_crit; - const double max_distance = e->mesh->r_cut_max; - const double max_distance2 = max_distance * max_distance; - - /* Some info about the zoom domain */ - const int bkg_cell_offset = s->zoom_props->tl_cell_offset; - - /* Compute how many cells away we need to walk */ - const double distance = 2.5 * cells[bkg_cell_offset].width[0] / theta_crit; - int delta = (int)(distance / cells[bkg_cell_offset].width[0]) + 1; - int delta_m = delta; - int delta_p = delta; - - /* Special case where every cell is in range of every other one */ - if (delta >= cdim[0] / 2) { - if (cdim[0] % 2 == 0) { - delta_m = cdim[0] / 2; - delta_p = cdim[0] / 2 - 1; - } else { - delta_m = cdim[0] / 2; - delta_p = cdim[0] / 2; - } - } - - /* Loop through the elements, which are just byte offsets from NULL. */ - for (int ind = 0; ind < num_elements; ind++) { - - /* Get the cell index, including background cell offset. */ - const int cid = (size_t)(map_data) + ind + bkg_cell_offset; - - /* Integer indices of the cell in the top-level grid */ - const int i = (cid - bkg_cell_offset) / (cdim[1] * cdim[2]); - const int j = ((cid - bkg_cell_offset) / cdim[2]) % cdim[1]; - const int k = (cid - bkg_cell_offset) % cdim[2]; - - /* Get the cell */ - struct cell *ci = &cells[cid]; - - /* Skip cells without gravity particles */ - if (ci->grav.count == 0) continue; - - /* Ensure we haven't found a void cell with particles */ - if (ci->tl_cell_type == void_tl_cell) - error("This void cell (cid=%d) has got particles!", cid); - - /* If the cell is local build a self-interaction */ - if (ci->nodeID == nodeID) { - scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, - NULL); - } - - /* Loop over every other cell within (Manhattan) range delta */ - for (int ii = -delta_m; ii <= delta_p; ii++) { - int iii = i + ii; - if (!periodic && (iii < 0 || iii >= cdim[0])) continue; - iii = (iii + cdim[0]) % cdim[0]; - for (int jj = -delta_m; jj <= delta_p; jj++) { - int jjj = j + jj; - if (!periodic && (jjj < 0 || jjj >= cdim[1])) continue; - jjj = (jjj + cdim[1]) % cdim[1]; - for (int kk = -delta_m; kk <= delta_p; kk++) { - int kkk = k + kk; - if (!periodic && (kkk < 0 || kkk >= cdim[2])) continue; - kkk = (kkk + cdim[2]) % cdim[2]; - - /* Get the cell */ - const int cjd = cell_getid(cdim, iii, jjj, kkk) + bkg_cell_offset; - struct cell *cj = &cells[cjd]; - - /* Avoid duplicates, empty cells and completely foreign pairs */ - if (cid >= cjd || cj->grav.count == 0 || - (ci->nodeID != nodeID && cj->nodeID != nodeID)) - continue; - - /* Recover the multipole information */ - const struct gravity_tensors *multi_i = ci->grav.multipole; - const struct gravity_tensors *multi_j = cj->grav.multipole; - - if (multi_i == NULL && ci->nodeID != nodeID) - error("Multipole of ci was not exchanged properly via the proxies (nodeID=%d, ci->nodeID=%d)", nodeID, ci->nodeID); - if (multi_j == NULL && cj->nodeID != nodeID) - error("Multipole of cj was not exchanged properly via the proxies (nodeID=%d, cj->nodeID=%d)", nodeID, cj->nodeID); - - /* Minimal distance between any pair of particles */ - const double min_radius2 = - cell_min_dist2_same_size(ci, cj, periodic, dim); - - /* Are we beyond the distance where the truncated forces are 0 ?*/ - if (periodic && min_radius2 > max_distance2) continue; - - /* Are the cells too close for a MM interaction ? */ - if (!cell_can_use_pair_mm(ci, cj, e, s, /*use_rebuild_data=*/1, - /*is_tree_walk=*/0)) { - - /* Ok, we need to add a direct pair calculation */ - scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, - ci, cj); +void engine_make_self_gravity_tasks_mapper_natural_cells(void *map_data, + int num_elements, + void *extra_data) { + + struct engine *e = (struct engine *)extra_data; + struct space *s = e->s; + struct scheduler *sched = &e->sched; + const int nodeID = e->nodeID; + const int periodic = s->periodic; + const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]}; + const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; + struct cell *cells = s->cells_top; + const double theta_crit = e->gravity_properties->theta_crit; + const double max_distance = e->mesh->r_cut_max; + const double max_distance2 = max_distance * max_distance; + + /* Some info about the zoom domain */ + const int bkg_cell_offset = s->zoom_props->tl_cell_offset; + + /* Compute how many cells away we need to walk */ + const double distance = 2.5 * cells[bkg_cell_offset].width[0] / theta_crit; + int delta = (int)(distance / cells[bkg_cell_offset].width[0]) + 1; + int delta_m = delta; + int delta_p = delta; + + /* Special case where every cell is in range of every other one */ + if (delta >= cdim[0] / 2) { + if (cdim[0] % 2 == 0) { + delta_m = cdim[0] / 2; + delta_p = cdim[0] / 2 - 1; + } else { + delta_m = cdim[0] / 2; + delta_p = cdim[0] / 2; + } + } + + /* Loop through the elements, which are just byte offsets from NULL. */ + for (int ind = 0; ind < num_elements; ind++) { + + /* Get the cell index, including background cell offset. */ + const int cid = (size_t)(map_data) + ind + bkg_cell_offset; + + /* Integer indices of the cell in the top-level grid */ + const int i = (cid - bkg_cell_offset) / (cdim[1] * cdim[2]); + const int j = ((cid - bkg_cell_offset) / cdim[2]) % cdim[1]; + const int k = (cid - bkg_cell_offset) % cdim[2]; + + /* Get the cell */ + struct cell *ci = &cells[cid]; + + /* Skip cells without gravity particles */ + if (ci->grav.count == 0) continue; + + /* Ensure we haven't found a void cell with particles */ + if (ci->tl_cell_type == void_tl_cell) + error("This void cell (cid=%d) has got particles!", cid); + + /* If the cell is local build a self-interaction */ + if (ci->nodeID == nodeID) { + scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, + NULL); + } + + /* Loop over every other cell within (Manhattan) range delta */ + for (int ii = -delta_m; ii <= delta_p; ii++) { + int iii = i + ii; + if (!periodic && (iii < 0 || iii >= cdim[0])) continue; + iii = (iii + cdim[0]) % cdim[0]; + for (int jj = -delta_m; jj <= delta_p; jj++) { + int jjj = j + jj; + if (!periodic && (jjj < 0 || jjj >= cdim[1])) continue; + jjj = (jjj + cdim[1]) % cdim[1]; + for (int kk = -delta_m; kk <= delta_p; kk++) { + int kkk = k + kk; + if (!periodic && (kkk < 0 || kkk >= cdim[2])) continue; + kkk = (kkk + cdim[2]) % cdim[2]; + + /* Get the cell */ + const int cjd = cell_getid(cdim, iii, jjj, kkk) + bkg_cell_offset; + struct cell *cj = &cells[cjd]; + + /* Avoid duplicates, empty cells and completely foreign pairs */ + if (cid >= cjd || cj->grav.count == 0 || + (ci->nodeID != nodeID && cj->nodeID != nodeID)) + continue; + + /* Recover the multipole information */ + const struct gravity_tensors *multi_i = ci->grav.multipole; + const struct gravity_tensors *multi_j = cj->grav.multipole; + + if (multi_i == NULL && ci->nodeID != nodeID) + error( + "Multipole of ci was not exchanged properly via the proxies " + "(nodeID=%d, ci->nodeID=%d)", + nodeID, ci->nodeID); + if (multi_j == NULL && cj->nodeID != nodeID) + error( + "Multipole of cj was not exchanged properly via the proxies " + "(nodeID=%d, cj->nodeID=%d)", + nodeID, cj->nodeID); + + /* Minimal distance between any pair of particles */ + const double min_radius2 = + cell_min_dist2_same_size(ci, cj, periodic, dim); + + /* Are we beyond the distance where the truncated forces are 0 ?*/ + if (periodic && min_radius2 > max_distance2) continue; + + /* Are the cells too close for a MM interaction ? */ + if (!cell_can_use_pair_mm(ci, cj, e, s, /*use_rebuild_data=*/1, + /*is_tree_walk=*/0)) { + + /* Ok, we need to add a direct pair calculation */ + scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, + ci, cj); #ifdef SWIFT_DEBUG_CHECKS - /* Ensure both cells are background cells */ - if (ci->tl_cell_type == 3 || cj->tl_cell_type == 3) { - error("Cell %d and cell %d are not background cells! (ci->tl_cell_type=%d, cj->tl_cell_type=%d)", - cid, cjd, ci->tl_cell_type, cj->tl_cell_type); - } - #ifdef WITH_MPI + /* Ensure both cells are background cells */ + if (ci->tl_cell_type == 3 || cj->tl_cell_type == 3) { + error( + "Cell %d and cell %d are not background cells! " + "(ci->tl_cell_type=%d, cj->tl_cell_type=%d)", + cid, cjd, ci->tl_cell_type, cj->tl_cell_type); + } +#ifdef WITH_MPI /* Let's cross-check that we had a proxy for that cell */ if (ci->nodeID == nodeID && cj->nodeID != engine_rank) { @@ -1698,11 +1766,11 @@ void engine_make_self_gravity_tasks_mapper_natural_cells(void *map_data, int num } #endif /* WITH_MPI */ #endif /* SWIFT_DEBUG_CHECKS */ - } - } - } - } - } + } + } + } + } + } } /** @@ -1713,105 +1781,115 @@ void engine_make_self_gravity_tasks_mapper_natural_cells(void *map_data, int num * - All pairs within range according to the multipole acceptance * criterion get a pair task. */ -void engine_make_self_gravity_tasks_mapper_zoom_cells(void *map_data, int num_elements, - void *extra_data) { - - struct engine *e = (struct engine *)extra_data; - struct space *s = e->s; - struct scheduler *sched = &e->sched; - const int nodeID = e->nodeID; - const int cdim[3] = {s->zoom_props->cdim[0], s->zoom_props->cdim[1], s->zoom_props->cdim[2]}; - struct cell *cells = s->cells_top; - const double theta_crit = e->gravity_properties->theta_crit; - - /* Compute how many cells away we need to walk */ - const double distance = 2.5 * cells[0].width[0] / theta_crit; - int delta = (int)(distance / cells[0].width[0]) + 1; - int delta_m = delta; - int delta_p = delta; - - /* Special case where every cell is in range of every other one */ - if (delta >= cdim[0] / 2) { - if (cdim[0] % 2 == 0) { - delta_m = cdim[0] / 2; - delta_p = cdim[0] / 2 - 1; - } else { - delta_m = cdim[0] / 2; - delta_p = cdim[0] / 2; - } - } - - /* Loop through the elements, which are just byte offsets from NULL. */ - for (int ind = 0; ind < num_elements; ind++) { - - /* Get the cell index. */ - const int cid = (size_t)(map_data) + ind; - - /* Integer indices of the cell in the top-level grid */ - const int i = cid / (cdim[1] * cdim[2]); - const int j = (cid / cdim[2]) % cdim[1]; - const int k = cid % cdim[2]; - - /* Get the cell */ - struct cell *ci = &cells[cid]; - - /* Skip cells without gravity particles */ - if (ci->grav.count == 0) continue; - - /* If the cell is local build a self-interaction */ - if (ci->nodeID == nodeID) { - scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, - NULL); - } - - /* Loop over every other cell within (Manhattan) range delta - * NOTE: Zoom cells are never periodic */ - for (int ii = -delta_m; ii <= delta_p; ii++) { - int iii = i + ii; - if (iii < 0 || iii >= cdim[0]) continue; - iii = (iii + cdim[0]) % cdim[0]; - for (int jj = -delta_m; jj <= delta_p; jj++) { - int jjj = j + jj; - if (jjj < 0 || jjj >= cdim[1]) continue; - jjj = (jjj + cdim[1]) % cdim[1]; - for (int kk = -delta_m; kk <= delta_p; kk++) { - int kkk = k + kk; - if (kkk < 0 || kkk >= cdim[2]) continue; - kkk = (kkk + cdim[2]) % cdim[2]; - - /* Get the cell */ - const int cjd = cell_getid(cdim, iii, jjj, kkk); - struct cell *cj = &cells[cjd]; - - /* Avoid duplicates, empty cells and completely foreign pairs */ - if (cid >= cjd || cj->grav.count == 0 || - (ci->nodeID != nodeID && cj->nodeID != nodeID)) - continue; - - /* Recover the multipole information */ - const struct gravity_tensors *multi_i = ci->grav.multipole; - const struct gravity_tensors *multi_j = cj->grav.multipole; - - if (multi_i == NULL && ci->nodeID != nodeID) - error("Multipole of ci was not exchanged properly via the proxies (nodeID=%d, ci->nodeID=%d)", nodeID, ci->nodeID); - if (multi_j == NULL && cj->nodeID != nodeID) - error("Multipole of cj was not exchanged properly via the proxies (nodeID=%d, cj->nodeID=%d)", nodeID, cj->nodeID); - - /* Are the cells too close for a MM interaction ? */ - if (!cell_can_use_pair_mm(ci, cj, e, s, /*use_rebuild_data=*/1, - /*is_tree_walk=*/0)) { - - /* Ok, we need to add a direct pair calculation */ - scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, - ci, cj); +void engine_make_self_gravity_tasks_mapper_zoom_cells(void *map_data, + int num_elements, + void *extra_data) { + + struct engine *e = (struct engine *)extra_data; + struct space *s = e->s; + struct scheduler *sched = &e->sched; + const int nodeID = e->nodeID; + const int cdim[3] = {s->zoom_props->cdim[0], s->zoom_props->cdim[1], + s->zoom_props->cdim[2]}; + struct cell *cells = s->cells_top; + const double theta_crit = e->gravity_properties->theta_crit; + + /* Compute how many cells away we need to walk */ + const double distance = 2.5 * cells[0].width[0] / theta_crit; + int delta = (int)(distance / cells[0].width[0]) + 1; + int delta_m = delta; + int delta_p = delta; + + /* Special case where every cell is in range of every other one */ + if (delta >= cdim[0] / 2) { + if (cdim[0] % 2 == 0) { + delta_m = cdim[0] / 2; + delta_p = cdim[0] / 2 - 1; + } else { + delta_m = cdim[0] / 2; + delta_p = cdim[0] / 2; + } + } + + /* Loop through the elements, which are just byte offsets from NULL. */ + for (int ind = 0; ind < num_elements; ind++) { + + /* Get the cell index. */ + const int cid = (size_t)(map_data) + ind; + + /* Integer indices of the cell in the top-level grid */ + const int i = cid / (cdim[1] * cdim[2]); + const int j = (cid / cdim[2]) % cdim[1]; + const int k = cid % cdim[2]; + + /* Get the cell */ + struct cell *ci = &cells[cid]; + + /* Skip cells without gravity particles */ + if (ci->grav.count == 0) continue; + + /* If the cell is local build a self-interaction */ + if (ci->nodeID == nodeID) { + scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, + NULL); + } + + /* Loop over every other cell within (Manhattan) range delta + * NOTE: Zoom cells are never periodic */ + for (int ii = -delta_m; ii <= delta_p; ii++) { + int iii = i + ii; + if (iii < 0 || iii >= cdim[0]) continue; + iii = (iii + cdim[0]) % cdim[0]; + for (int jj = -delta_m; jj <= delta_p; jj++) { + int jjj = j + jj; + if (jjj < 0 || jjj >= cdim[1]) continue; + jjj = (jjj + cdim[1]) % cdim[1]; + for (int kk = -delta_m; kk <= delta_p; kk++) { + int kkk = k + kk; + if (kkk < 0 || kkk >= cdim[2]) continue; + kkk = (kkk + cdim[2]) % cdim[2]; + + /* Get the cell */ + const int cjd = cell_getid(cdim, iii, jjj, kkk); + struct cell *cj = &cells[cjd]; + + /* Avoid duplicates, empty cells and completely foreign pairs */ + if (cid >= cjd || cj->grav.count == 0 || + (ci->nodeID != nodeID && cj->nodeID != nodeID)) + continue; + + /* Recover the multipole information */ + const struct gravity_tensors *multi_i = ci->grav.multipole; + const struct gravity_tensors *multi_j = cj->grav.multipole; + + if (multi_i == NULL && ci->nodeID != nodeID) + error( + "Multipole of ci was not exchanged properly via the proxies " + "(nodeID=%d, ci->nodeID=%d)", + nodeID, ci->nodeID); + if (multi_j == NULL && cj->nodeID != nodeID) + error( + "Multipole of cj was not exchanged properly via the proxies " + "(nodeID=%d, cj->nodeID=%d)", + nodeID, cj->nodeID); + + /* Are the cells too close for a MM interaction ? */ + if (!cell_can_use_pair_mm(ci, cj, e, s, /*use_rebuild_data=*/1, + /*is_tree_walk=*/0)) { + + /* Ok, we need to add a direct pair calculation */ + scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, + ci, cj); #ifdef SWIFT_DEBUG_CHECKS - /* Ensure both cells are zoom cells */ - if (ci->tl_cell_type <= 2 || cj->tl_cell_type <= 2) { - error("Cell %d and cell %d are not zoom cells! (ci->tl_cell_type=%d, cj->tl_cell_type=%d)", - cid, cjd, ci->tl_cell_type, cj->tl_cell_type); - } - #ifdef WITH_MPI + /* Ensure both cells are zoom cells */ + if (ci->tl_cell_type <= 2 || cj->tl_cell_type <= 2) { + error( + "Cell %d and cell %d are not zoom cells! " + "(ci->tl_cell_type=%d, cj->tl_cell_type=%d)", + cid, cjd, ci->tl_cell_type, cj->tl_cell_type); + } +#ifdef WITH_MPI /* Let's cross-check that we had a proxy for that cell */ if (ci->nodeID == nodeID && cj->nodeID != engine_rank) { @@ -1857,11 +1935,11 @@ void engine_make_self_gravity_tasks_mapper_zoom_cells(void *map_data, int num_el } #endif /* WITH_MPI */ #endif /* SWIFT_DEBUG_CHECKS */ - } - } - } - } - } + } + } + } + } + } } /** @@ -1869,7 +1947,8 @@ void engine_make_self_gravity_tasks_mapper_zoom_cells(void *map_data, int num_el * and long-range gravity interactions between natural level cells * and zoom level cells. * - * This replaces the function in engine_maketasks when running with a zoom region. + * This replaces the function in engine_maketasks when running with a zoom + region. * * - All top-cells get a self task. * - All pairs of differing sized cells within range according to @@ -1880,167 +1959,176 @@ void engine_make_self_gravity_tasks_mapper_zoom_cells(void *map_data, int num_el * @param extra_data The #engine. */ -void engine_make_self_gravity_tasks_mapper_with_zoom_diffsize(void *map_data, - int num_elements, - void *extra_data) { +void engine_make_self_gravity_tasks_mapper_with_zoom_diffsize( + void *map_data, int num_elements, void *extra_data) { - /* Useful local information */ - struct engine *e = (struct engine *)extra_data; - struct space *s = e->s; - struct scheduler *sched = &e->sched; - const int nodeID = e->nodeID; + /* Useful local information */ + struct engine *e = (struct engine *)extra_data; + struct space *s = e->s; + struct scheduler *sched = &e->sched; + const int nodeID = e->nodeID; - /* Handle on the cells and proxies */ - struct cell *cells = s->cells_top; + /* Handle on the cells and proxies */ + struct cell *cells = s->cells_top; - /* Some info about the zoom domain */ - const int bkg_cell_offset = s->zoom_props->tl_cell_offset; + /* Some info about the zoom domain */ + const int bkg_cell_offset = s->zoom_props->tl_cell_offset; - /* Some info about the domain */ - const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]}; - int periodic = s->periodic; + /* Some info about the domain */ + const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]}; + int periodic = s->periodic; - /* Get some info about the physics */ - const double max_mesh_dist = e->mesh->r_cut_max; - const double max_mesh_dist2 = max_mesh_dist * max_mesh_dist; + /* Get some info about the physics */ + const double max_mesh_dist = e->mesh->r_cut_max; + const double max_mesh_dist2 = max_mesh_dist * max_mesh_dist; - /* Define neighbour loop variables */ - int cjd_low = 0; - int cjd_high = bkg_cell_offset; + /* Define neighbour loop variables */ + int cjd_low = 0; + int cjd_high = bkg_cell_offset; - /* Loop through the elements, which are just byte offsets from NULL. */ - for (int ind = 0; ind < num_elements; ind++) { + /* Loop through the elements, which are just byte offsets from NULL. */ + for (int ind = 0; ind < num_elements; ind++) { - /* Get the cell index. */ - const int cid = (size_t)(map_data) + ind; + /* Get the cell index. */ + const int cid = (size_t)(map_data) + ind; - /* Get the cell */ - struct cell *ci = &cells[cid]; + /* Get the cell */ + struct cell *ci = &cells[cid]; /* If this cell is on this node and is a background cell * then we have to avoid duplicating tasks */ if (ci->nodeID == nodeID && ci->tl_cell_type <= 2) { - continue; - } - - /* Skip cells without gravity particles */ - if (ci->grav.count == 0) continue; - - /* If the cell is a natural cell and not a neighbour cell - * we don't need to do anything */ - if ((ci->tl_cell_type <= 2) && (ci->tl_cell_type != tl_cell_neighbour)) { - continue; - } - - /* Get the loop range for the neighbouring cells */ - if (ci->tl_cell_type <= 2) { - cjd_low = 0; - cjd_high = bkg_cell_offset; - } else { - cjd_low = bkg_cell_offset; - cjd_high = s->nr_cells; - } - - /* Loop over every other cell */ - for (int cjd = cjd_low; cjd < cjd_high; cjd++) { - - /* Get the cell */ - struct cell *cj = &cells[cjd]; - - /* If the cell is a natural cell and not a neighbour cell - * we don't need to do anything */ - if ((cj->tl_cell_type <= 2) && (cj->tl_cell_type != tl_cell_neighbour)) { - continue; - } - - /* Avoid empty cells and completely foreign pairs */ - if (cj->grav.count == 0 || (ci->nodeID != nodeID && cj->nodeID != nodeID)) - continue; - - /* Explictly avoid duplicates */ - if ((ci->nodeID == cj->nodeID && cid >= cjd) || (cj->nodeID == nodeID && cj->tl_cell_type == zoom_tl_cell)) { - continue; - } - - /* Recover the multipole information */ - const struct gravity_tensors *multi_i = ci->grav.multipole; - const struct gravity_tensors *multi_j = cj->grav.multipole; - - if (multi_i == NULL && ci->nodeID != nodeID) - error("Multipole of ci was not exchanged properly via the proxies (nodeID=%d, ci->nodeID=%d)", nodeID, ci->nodeID); - if (multi_j == NULL && cj->nodeID != nodeID) - error("Multipole of cj was not exchanged properly via the proxies (nodeID=%d, cj->nodeID=%d)", nodeID, cj->nodeID); - - /* Minimal distance between any pair of particles */ - const double min_radius2 = cell_min_dist2_diff_size(ci, cj, periodic, dim); - - /* Are we beyond the distance where the truncated forces are 0 ?*/ - if (periodic && min_radius2 > max_mesh_dist2) continue; - - /* Are the cells too close for a MM interaction ? */ - if (!cell_can_use_pair_mm(ci, cj, e, s, /*use_rebuild_data=*/1, - /*is_tree_walk=*/0)) { - - /* Ok, we need to add a direct pair calculation */ - scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, - ci, cj); + continue; + } + + /* Skip cells without gravity particles */ + if (ci->grav.count == 0) continue; + + /* If the cell is a natural cell and not a neighbour cell + * we don't need to do anything */ + if ((ci->tl_cell_type <= 2) && (ci->tl_cell_type != tl_cell_neighbour)) { + continue; + } + + /* Get the loop range for the neighbouring cells */ + if (ci->tl_cell_type <= 2) { + cjd_low = 0; + cjd_high = bkg_cell_offset; + } else { + cjd_low = bkg_cell_offset; + cjd_high = s->nr_cells; + } + + /* Loop over every other cell */ + for (int cjd = cjd_low; cjd < cjd_high; cjd++) { + + /* Get the cell */ + struct cell *cj = &cells[cjd]; + + /* If the cell is a natural cell and not a neighbour cell + * we don't need to do anything */ + if ((cj->tl_cell_type <= 2) && (cj->tl_cell_type != tl_cell_neighbour)) { + continue; + } + + /* Avoid empty cells and completely foreign pairs */ + if (cj->grav.count == 0 || (ci->nodeID != nodeID && cj->nodeID != nodeID)) + continue; + + /* Explictly avoid duplicates */ + if ((ci->nodeID == cj->nodeID && cid >= cjd) || + (cj->nodeID == nodeID && cj->tl_cell_type == zoom_tl_cell)) { + continue; + } + + /* Recover the multipole information */ + const struct gravity_tensors *multi_i = ci->grav.multipole; + const struct gravity_tensors *multi_j = cj->grav.multipole; + + if (multi_i == NULL && ci->nodeID != nodeID) + error( + "Multipole of ci was not exchanged properly via the proxies " + "(nodeID=%d, ci->nodeID=%d)", + nodeID, ci->nodeID); + if (multi_j == NULL && cj->nodeID != nodeID) + error( + "Multipole of cj was not exchanged properly via the proxies " + "(nodeID=%d, cj->nodeID=%d)", + nodeID, cj->nodeID); + + /* Minimal distance between any pair of particles */ + const double min_radius2 = + cell_min_dist2_diff_size(ci, cj, periodic, dim); + + /* Are we beyond the distance where the truncated forces are 0 ?*/ + if (periodic && min_radius2 > max_mesh_dist2) continue; + + /* Are the cells too close for a MM interaction ? */ + if (!cell_can_use_pair_mm(ci, cj, e, s, /*use_rebuild_data=*/1, + /*is_tree_walk=*/0)) { + + /* Ok, we need to add a direct pair calculation */ + scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci, + cj); #ifdef SWIFT_DEBUG_CHECKS - /* Ensure both cells are not in the same level */ - if (((ci->tl_cell_type <= 2 && cj->tl_cell_type <= 2) || - (ci->tl_cell_type == cj->tl_cell_type))) { - error("Cell %d and cell %d are the same cell type! (ci->tl_cell_type=%d, cj->tl_cell_type=%d)", - cid, cjd, ci->tl_cell_type, cj->tl_cell_type); - } + /* Ensure both cells are not in the same level */ + if (((ci->tl_cell_type <= 2 && cj->tl_cell_type <= 2) || + (ci->tl_cell_type == cj->tl_cell_type))) { + error( + "Cell %d and cell %d are the same cell type! " + "(ci->tl_cell_type=%d, cj->tl_cell_type=%d)", + cid, cjd, ci->tl_cell_type, cj->tl_cell_type); + } #ifdef WITH_MPI - /* Let's cross-check that we had a proxy for that cell */ - if (ci->nodeID == nodeID && cj->nodeID != engine_rank) { - - /* Find the proxy for this node */ - const int proxy_id = e->proxy_ind[cj->nodeID]; - if (proxy_id < 0) - error("No proxy exists for that foreign node %d!", cj->nodeID); - - const struct proxy *p = &e->proxies[proxy_id]; - - /* Check whether the cell exists in the proxy */ - int n = 0; - for (; n < p->nr_cells_in; n++) - if (p->cells_in[n] == cj) { - break; - } - if (n == p->nr_cells_in) - error( - "Cell %d not found in the proxy but trying to construct " - "grav task!", - cjd); - } else if (cj->nodeID == nodeID && ci->nodeID != engine_rank) { - - /* Find the proxy for this node */ - const int proxy_id = e->proxy_ind[ci->nodeID]; - if (proxy_id < 0) - error("No proxy exists for that foreign node %d!", ci->nodeID); - - const struct proxy *p = &e->proxies[proxy_id]; - - /* Check whether the cell exists in the proxy */ - int n = 0; - for (; n < p->nr_cells_in; n++) - if (p->cells_in[n] == ci) { - break; - } - if (n == p->nr_cells_in) - error( - "Cell %d not found in the proxy but trying to construct " - "grav task!", - cid); - } + /* Let's cross-check that we had a proxy for that cell */ + if (ci->nodeID == nodeID && cj->nodeID != engine_rank) { + + /* Find the proxy for this node */ + const int proxy_id = e->proxy_ind[cj->nodeID]; + if (proxy_id < 0) + error("No proxy exists for that foreign node %d!", cj->nodeID); + + const struct proxy *p = &e->proxies[proxy_id]; + + /* Check whether the cell exists in the proxy */ + int n = 0; + for (; n < p->nr_cells_in; n++) + if (p->cells_in[n] == cj) { + break; + } + if (n == p->nr_cells_in) + error( + "Cell %d not found in the proxy but trying to construct " + "grav task!", + cjd); + } else if (cj->nodeID == nodeID && ci->nodeID != engine_rank) { + + /* Find the proxy for this node */ + const int proxy_id = e->proxy_ind[ci->nodeID]; + if (proxy_id < 0) + error("No proxy exists for that foreign node %d!", ci->nodeID); + + const struct proxy *p = &e->proxies[proxy_id]; + + /* Check whether the cell exists in the proxy */ + int n = 0; + for (; n < p->nr_cells_in; n++) + if (p->cells_in[n] == ci) { + break; + } + if (n == p->nr_cells_in) + error( + "Cell %d not found in the proxy but trying to construct " + "grav task!", + cid); + } #endif /* WITH_MPI */ #endif /* SWIFT_DEBUG_CHECKS */ - } - } - } + } + } + } } /** @@ -2059,8 +2147,9 @@ void engine_make_self_gravity_tasks_mapper_with_zoom_diffsize(void *map_data, * @param num_elements Number of cells to traverse. * @param extra_data The #engine. */ -void engine_make_hydroloop_tasks_mapper_with_zoom(void *map_data, int num_elements, - void *extra_data) { +void engine_make_hydroloop_tasks_mapper_with_zoom(void *map_data, + int num_elements, + void *extra_data) { /* Extract the engine pointer. */ struct engine *e = (struct engine *)extra_data; @@ -2198,8 +2287,9 @@ void engine_make_hydroloop_tasks_mapper_with_zoom(void *map_data, int num_elemen * @param num_elements Number of cells to traverse. * @param extra_data The #engine. */ -void engine_make_fofloop_tasks_mapper_with_zoom(void *map_data, int num_elements, - void *extra_data) { +void engine_make_fofloop_tasks_mapper_with_zoom(void *map_data, + int num_elements, + void *extra_data) { /* Extract the engine pointer. */ struct engine *e = (struct engine *)extra_data; @@ -2263,5 +2353,4 @@ void engine_make_fofloop_tasks_mapper_with_zoom(void *map_data, int num_elements } } - #endif /* WITH_ZOOM_REGION */ diff --git a/src/zoom_region.h b/src/zoom_region.h index 0a7cac0870..6c6c68b73f 100644 --- a/src/zoom_region.h +++ b/src/zoom_region.h @@ -5,13 +5,18 @@ #ifndef SWIFT_ZOOM_H #define SWIFT_ZOOM_H -void zoom_region_init(struct swift_params *params, struct space *s, int verbose); +void zoom_region_init(struct swift_params *params, struct space *s, + int verbose); int cell_getid_zoom(const struct space *s, const double x, const double y, const double z); void construct_zoom_region(struct space *s, int verbose); -void construct_tl_cells_with_zoom_region(struct space *s, const int *cdim, const float dmin, - const integertime_t ti_current, struct gravity_props *gravity_properties, int verbose); -void find_neighbouring_cells(struct space *s, struct gravity_props *gravity_properties, const int verbose); +void construct_tl_cells_with_zoom_region( + struct space *s, const int *cdim, const float dmin, + const integertime_t ti_current, struct gravity_props *gravity_properties, + int verbose); +void find_neighbouring_cells(struct space *s, + struct gravity_props *gravity_properties, + const int verbose); double cell_min_dist2_diff_size(const struct cell *restrict ci, const struct cell *restrict cj, const int periodic, const double dim[3]); @@ -23,25 +28,29 @@ void engine_makeproxies_with_zoom_cells(struct engine *e); void engine_makeproxies_with_between_grids(struct engine *e); void engine_makeproxies_with_zoom_region(struct engine *e); void engine_make_self_gravity_tasks_mapper_natural_cells(void *map_data, - int num_elements, - void *extra_data); + int num_elements, + void *extra_data); void engine_make_self_gravity_tasks_mapper_zoom_cells(void *map_data, int num_elements, void *extra_data); void engine_make_self_gravity_tasks_mapper_with_zoom_diffsize(void *map_data, int num_elements, void *extra_data); -void engine_make_hydroloop_tasks_mapper_with_zoom(void *map_data, int num_elements, +void engine_make_hydroloop_tasks_mapper_with_zoom(void *map_data, + int num_elements, void *extra_data); void engine_make_drift_tasks_for_wanderers(struct engine *e, struct cell *c); -void engine_make_drift_tasks_for_wanderers_mapper(void *map_data, int num_elements, +void engine_make_drift_tasks_for_wanderers_mapper(void *map_data, + int num_elements, void *extra_data); -void engine_make_fofloop_tasks_mapper_with_zoom(void *map_data, int num_elements, +void engine_make_fofloop_tasks_mapper_with_zoom(void *map_data, + int num_elements, void *extra_data); /* Parition prototypes */ -int partition_space_to_space_zoom(double *oldh, double *oldcdim, double *oldzoomh, - double *oldzoomcdim, int *oldnodeIDs, struct space *s); +int partition_space_to_space_zoom(double *oldh, double *oldcdim, + double *oldzoomh, double *oldzoomcdim, + int *oldnodeIDs, struct space *s); void pick_vector_zoom(struct space *s, int nregions, int *samplecells); void split_vector_zoom(struct space *s, int nregions, int *samplecells); diff --git a/src/zoom_regrid.c b/src/zoom_regrid.c index 29257d7f45..6da01565bb 100644 --- a/src/zoom_regrid.c +++ b/src/zoom_regrid.c @@ -29,8 +29,8 @@ /* Local headers. */ #include "cell.h" -#include "gravity_properties.h" #include "engine.h" +#include "gravity_properties.h" #include "scheduler.h" #include "zoom_region.h" @@ -41,86 +41,89 @@ * @param s The #space. * @param verbose Print messages to stdout or not. */ -void space_regrid_zoom(struct space *s, struct gravity_props *gravity_properties, int verbose) { - - const size_t nr_parts = s->nr_parts; - const size_t nr_sparts = s->nr_sparts; - const size_t nr_bparts = s->nr_bparts; - const size_t nr_sinks = s->nr_sinks; - const ticks tic = getticks(); - const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0; - - /* Run through the cells and get the current h_max, when using a zoom region - * h_max needs to be set by the zoom cells. */ - // tic = getticks(); - const double zoom_cell_min = s->zoom_props->cell_min; - float h_max = zoom_cell_min / kernel_gamma / space_stretch; - if (nr_parts > 0) { - - /* Can we use the list of local non-empty top-level cells? */ - if (s->local_cells_with_particles_top != NULL) { - for (int k = 0; k < s->nr_local_cells_with_particles; ++k) { - const struct cell *c = - &s->cells_top[s->local_cells_with_particles_top[k]]; - - /* We only want to consider zoom cells to avoid over inflating the smoothing length */ - if (c->tl_cell_type != zoom_tl_cell) continue; - - if (c->hydro.h_max > h_max) { - h_max = c->hydro.h_max; - } - if (c->stars.h_max > h_max) { - h_max = c->stars.h_max; - } - if (c->black_holes.h_max > h_max) { - h_max = c->black_holes.h_max; - } - if (c->sinks.r_cut_max > h_max) { - h_max = c->sinks.r_cut_max / kernel_gamma; - } - } - - /* Can we instead use all the top-level cells? */ - } else if (s->cells_top != NULL) { - /* We only want to consider zoom cells to avoid over inflating the smoothing length */ - for (int k = s->zoom_props->tl_cell_offset; k < s->nr_cells; k++) { - const struct cell *c = &s->cells_top[k]; - - if (c->nodeID == engine_rank && c->hydro.h_max > h_max) { - h_max = c->hydro.h_max; - } - if (c->nodeID == engine_rank && c->stars.h_max > h_max) { - h_max = c->stars.h_max; - } - if (c->nodeID == engine_rank && c->black_holes.h_max > h_max) { - h_max = c->black_holes.h_max; - } - if (c->nodeID == engine_rank && c->sinks.r_cut_max > h_max) { - h_max = c->sinks.r_cut_max / kernel_gamma; - } - } - - /* Last option: run through the particles */ - } else { - for (size_t k = 0; k < nr_parts; k++) { - if (s->parts[k].h > h_max) h_max = s->parts[k].h; - } - for (size_t k = 0; k < nr_sparts; k++) { - if (s->sparts[k].h > h_max) h_max = s->sparts[k].h; - } - for (size_t k = 0; k < nr_bparts; k++) { - if (s->bparts[k].h > h_max) h_max = s->bparts[k].h; - } - for (size_t k = 0; k < nr_sinks; k++) { - if (s->sinks[k].r_cut > h_max) h_max = s->sinks[k].r_cut / kernel_gamma; - } - } - } +void space_regrid_zoom(struct space *s, + struct gravity_props *gravity_properties, int verbose) { + + const size_t nr_parts = s->nr_parts; + const size_t nr_sparts = s->nr_sparts; + const size_t nr_bparts = s->nr_bparts; + const size_t nr_sinks = s->nr_sinks; + const ticks tic = getticks(); + const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0; + + /* Run through the cells and get the current h_max, when using a zoom region + * h_max needs to be set by the zoom cells. */ + // tic = getticks(); + const double zoom_cell_min = s->zoom_props->cell_min; + float h_max = zoom_cell_min / kernel_gamma / space_stretch; + if (nr_parts > 0) { + + /* Can we use the list of local non-empty top-level cells? */ + if (s->local_cells_with_particles_top != NULL) { + for (int k = 0; k < s->nr_local_cells_with_particles; ++k) { + const struct cell *c = + &s->cells_top[s->local_cells_with_particles_top[k]]; + + /* We only want to consider zoom cells to avoid over inflating the + * smoothing length */ + if (c->tl_cell_type != zoom_tl_cell) continue; + + if (c->hydro.h_max > h_max) { + h_max = c->hydro.h_max; + } + if (c->stars.h_max > h_max) { + h_max = c->stars.h_max; + } + if (c->black_holes.h_max > h_max) { + h_max = c->black_holes.h_max; + } + if (c->sinks.r_cut_max > h_max) { + h_max = c->sinks.r_cut_max / kernel_gamma; + } + } + + /* Can we instead use all the top-level cells? */ + } else if (s->cells_top != NULL) { + /* We only want to consider zoom cells to avoid over inflating the + * smoothing length */ + for (int k = s->zoom_props->tl_cell_offset; k < s->nr_cells; k++) { + const struct cell *c = &s->cells_top[k]; + + if (c->nodeID == engine_rank && c->hydro.h_max > h_max) { + h_max = c->hydro.h_max; + } + if (c->nodeID == engine_rank && c->stars.h_max > h_max) { + h_max = c->stars.h_max; + } + if (c->nodeID == engine_rank && c->black_holes.h_max > h_max) { + h_max = c->black_holes.h_max; + } + if (c->nodeID == engine_rank && c->sinks.r_cut_max > h_max) { + h_max = c->sinks.r_cut_max / kernel_gamma; + } + } + + /* Last option: run through the particles */ + } else { + for (size_t k = 0; k < nr_parts; k++) { + if (s->parts[k].h > h_max) h_max = s->parts[k].h; + } + for (size_t k = 0; k < nr_sparts; k++) { + if (s->sparts[k].h > h_max) h_max = s->sparts[k].h; + } + for (size_t k = 0; k < nr_bparts; k++) { + if (s->bparts[k].h > h_max) h_max = s->bparts[k].h; + } + for (size_t k = 0; k < nr_sinks; k++) { + if (s->sinks[k].r_cut > h_max) h_max = s->sinks[k].r_cut / kernel_gamma; + } + } + } /* If we are running in parallel, make sure everybody agrees on how large the largest cell should be. */ #ifdef WITH_MPI - { + { float buff; if (MPI_Allreduce(&h_max, &buff, 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD) != MPI_SUCCESS) @@ -128,18 +131,19 @@ void space_regrid_zoom(struct space *s, struct gravity_props *gravity_properties h_max = buff; } #endif - if (verbose) message("h_max is %.3e (zoom_cell_min=%.3e).", h_max, zoom_cell_min); - - /* Get the new putative zoom cell dimensions. */ - const double dmax = max3(s->zoom_props->dim[0], - s->zoom_props->dim[1], - s->zoom_props->dim[2]); - const int zoom_cdim[3] = {(int)floor((dmax + 0.01 * zoom_cell_min) / - fmax(h_max * kernel_gamma * space_stretch, zoom_cell_min)), - (int)floor((dmax + 0.01 * zoom_cell_min) / - fmax(h_max * kernel_gamma * space_stretch, zoom_cell_min)), - (int)floor((dmax + 0.01 * zoom_cell_min) / - fmax(h_max * kernel_gamma * space_stretch, zoom_cell_min))}; + if (verbose) + message("h_max is %.3e (zoom_cell_min=%.3e).", h_max, zoom_cell_min); + + /* Get the new putative zoom cell dimensions. */ + const double dmax = + max3(s->zoom_props->dim[0], s->zoom_props->dim[1], s->zoom_props->dim[2]); + const int zoom_cdim[3] = { + (int)floor((dmax + 0.01 * zoom_cell_min) / + fmax(h_max * kernel_gamma * space_stretch, zoom_cell_min)), + (int)floor((dmax + 0.01 * zoom_cell_min) / + fmax(h_max * kernel_gamma * space_stretch, zoom_cell_min)), + (int)floor((dmax + 0.01 * zoom_cell_min) / + fmax(h_max * kernel_gamma * space_stretch, zoom_cell_min))}; /* In MPI-Land, changing the top-level cell size requires that the * global partition is recomputed and the particles redistributed. @@ -200,149 +204,172 @@ void space_regrid_zoom(struct space *s, struct gravity_props *gravity_properties const int no_regrid = (s->cells_top == NULL && oldnodeIDs == NULL); #endif /* WITH_MPI */ - /* Do we need to re-build the upper-level cells? */ - // tic = getticks(); - if (s->cells_top == NULL || - zoom_cdim[0] < s->zoom_props->cdim[0] || + /* Do we need to re-build the upper-level cells? */ + // tic = getticks(); + if (s->cells_top == NULL || zoom_cdim[0] < s->zoom_props->cdim[0] || zoom_cdim[1] < s->zoom_props->cdim[1] || zoom_cdim[2] < s->zoom_props->cdim[2]) { - /* Free the old cells, if they were allocated. */ - if (s->cells_top != NULL) { - space_free_cells(s); - swift_free("local_cells_with_tasks_top", s->local_cells_with_tasks_top); - swift_free("local_cells_top", s->local_cells_top); - swift_free("local_zoom_cells_top", s->zoom_props->local_zoom_cells_top); - swift_free("local_bkg_cells_top", s->zoom_props->local_bkg_cells_top); - swift_free("local_zoom_cells_with_particles_top", s->zoom_props->local_zoom_cells_with_particles_top); - swift_free("local_bkg_cell_with_particless_top", s->zoom_props->local_bkg_cells_with_particles_top); - swift_free("cells_with_particles_top", s->cells_with_particles_top); - swift_free("local_cells_with_particles_top", - s->local_cells_with_particles_top); - swift_free("cells_top", s->cells_top); - swift_free("multipoles_top", s->multipoles_top); - } - - /* Also free the task arrays, these will be regenerated and we can use the - * memory while copying the particle arrays. */ - if (s->e != NULL) scheduler_free_tasks(&s->e->sched); - - /* Setting the new zoom cdim. */ + /* Free the old cells, if they were allocated. */ + if (s->cells_top != NULL) { + space_free_cells(s); + swift_free("local_cells_with_tasks_top", s->local_cells_with_tasks_top); + swift_free("local_cells_top", s->local_cells_top); + swift_free("local_zoom_cells_top", s->zoom_props->local_zoom_cells_top); + swift_free("local_bkg_cells_top", s->zoom_props->local_bkg_cells_top); + swift_free("local_zoom_cells_with_particles_top", + s->zoom_props->local_zoom_cells_with_particles_top); + swift_free("local_bkg_cell_with_particless_top", + s->zoom_props->local_bkg_cells_with_particles_top); + swift_free("cells_with_particles_top", s->cells_with_particles_top); + swift_free("local_cells_with_particles_top", + s->local_cells_with_particles_top); + swift_free("cells_top", s->cells_top); + swift_free("multipoles_top", s->multipoles_top); + } + + /* Also free the task arrays, these will be regenerated and we can use the + * memory while copying the particle arrays. */ + if (s->e != NULL) scheduler_free_tasks(&s->e->sched); + + /* Setting the new zoom cdim. */ for (int ijk = 0; ijk < 3; ijk++) { s->zoom_props->cdim[ijk] = zoom_cdim[ijk]; } - message("Constructing zoom region."); + message("Constructing zoom region."); /* Compute the bounds of the zoom region from the DM particles * and define the background cells based on this. */ construct_zoom_region(s, verbose); /* Be verbose about this. */ #ifdef SWIFT_DEBUG_CHECKS - message("(re)griding space cdim=(%d %d %d) zoom_cdim=(%d %d %d)", s->cdim[0], s->cdim[1], s->cdim[2], - s->zoom_props->cdim[0], s->zoom_props->cdim[1], s->zoom_props->cdim[2]); + message("(re)griding space cdim=(%d %d %d) zoom_cdim=(%d %d %d)", + s->cdim[0], s->cdim[1], s->cdim[2], s->zoom_props->cdim[0], + s->zoom_props->cdim[1], s->zoom_props->cdim[2]); fflush(stdout); #endif - /* Allocate the highest level of cells. */ - s->tot_cells = s->nr_cells = (s->cdim[0] * s->cdim[1] * s->cdim[2]) + - (s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * s->zoom_props->cdim[2]); - - if (swift_memalign("cells_top", (void **)&s->cells_top, cell_align, - s->nr_cells * sizeof(struct cell)) != 0) - error("Failed to allocate top-level cells."); - bzero(s->cells_top, s->nr_cells * sizeof(struct cell)); - - /* Allocate the multipoles for the top-level cells. */ - if (s->with_self_gravity) { - if (swift_memalign("multipoles_top", (void **)&s->multipoles_top, - multipole_align, - s->nr_cells * sizeof(struct gravity_tensors)) != 0) - error("Failed to allocate top-level multipoles."); - bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors)); - } - - /* Allocate the indices of local cells */ - if (swift_memalign("local_cells_top", (void **)&s->local_cells_top, - SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) - error("Failed to allocate indices of local top-level cells."); - bzero(s->local_cells_top, s->nr_cells * sizeof(int)); - - /* Allocate the indices of local zoom cells */ - if (swift_memalign("local_zoom_cells_top", (void **)&s->zoom_props->local_zoom_cells_top, - SWIFT_STRUCT_ALIGNMENT, s->zoom_props->nr_zoom_cells * sizeof(int)) != 0) - error("Failed to allocate indices of local top-level zoom cells."); - bzero(s->zoom_props->local_zoom_cells_top, s->zoom_props->nr_zoom_cells * sizeof(int)); - - /* Allocate the indices of local bkg cells */ - if (swift_memalign("local_bkg_cells_top", (void **)&s->zoom_props->local_bkg_cells_top, - SWIFT_STRUCT_ALIGNMENT, s->zoom_props->nr_bkg_cells * sizeof(int)) != 0) - error("Failed to allocate indices of local top-level background cells."); - bzero(s->zoom_props->local_bkg_cells_top, s->zoom_props->nr_bkg_cells * sizeof(int)); - - /* Allocate the indices of local zoom cells with particles */ - if (swift_memalign("local_zoom_cells_with_particles_top", (void **)&s->zoom_props->local_zoom_cells_with_particles_top, - SWIFT_STRUCT_ALIGNMENT, s->zoom_props->nr_zoom_cells * sizeof(int)) != 0) - error("Failed to allocate indices of local top-level zoom cells with particles."); - bzero(s->zoom_props->local_zoom_cells_with_particles_top, s->zoom_props->nr_zoom_cells * sizeof(int)); - - /* Allocate the indices of local bkg cells with particles */ - if (swift_memalign("local_bkg_cells_with_particles_top", (void **)&s->zoom_props->local_bkg_cells_with_particles_top, - SWIFT_STRUCT_ALIGNMENT, s->zoom_props->nr_bkg_cells * sizeof(int)) != 0) - error("Failed to allocate indices of local top-level background cells with particles."); - bzero(s->zoom_props->local_bkg_cells_with_particles_top, s->zoom_props->nr_bkg_cells * sizeof(int)); - - /* Allocate the indices of local cells with tasks */ - if (swift_memalign("local_cells_with_tasks_top", - (void **)&s->local_cells_with_tasks_top, - SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) - error("Failed to allocate indices of local top-level cells with tasks."); - bzero(s->local_cells_with_tasks_top, s->nr_cells * sizeof(int)); - - /* Allocate the indices of cells with particles */ - if (swift_memalign("cells_with_particles_top", - (void **)&s->cells_with_particles_top, - SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) - error("Failed to allocate indices of top-level cells with particles."); - bzero(s->cells_with_particles_top, s->nr_cells * sizeof(int)); - - /* Allocate the indices of local cells with particles */ - if (swift_memalign("local_cells_with_particles_top", - (void **)&s->local_cells_with_particles_top, - SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) - error( - "Failed to allocate indices of local top-level cells with " - "particles."); - bzero(s->local_cells_with_particles_top, s->nr_cells * sizeof(int)); - - /* Set the cells' locks */ - for (int k = 0; k < s->nr_cells; k++) { - if (lock_init(&s->cells_top[k].hydro.lock) != 0) - error("Failed to init spinlock for hydro."); - if (lock_init(&s->cells_top[k].grav.plock) != 0) - error("Failed to init spinlock for gravity."); - if (lock_init(&s->cells_top[k].grav.mlock) != 0) - error("Failed to init spinlock for multipoles."); - if (lock_init(&s->cells_top[k].grav.star_formation_lock) != 0) - error("Failed to init spinlock for star formation (gpart)."); - if (lock_init(&s->cells_top[k].stars.lock) != 0) - error("Failed to init spinlock for stars."); - if (lock_init(&s->cells_top[k].sinks.lock) != 0) - error("Failed to init spinlock for sinks."); - if (lock_init(&s->cells_top[k].sinks.sink_formation_lock) != 0) - error("Failed to init spinlock for sink formation."); - if (lock_init(&s->cells_top[k].black_holes.lock) != 0) - error("Failed to init spinlock for black holes."); - if (lock_init(&s->cells_top[k].stars.star_formation_lock) != 0) - error("Failed to init spinlock for star formation (spart)."); - } - - /* Construct both grids of cells */ - const double dmin = min3(s->width[0], s->width[1], s->width[2]); - construct_tl_cells_with_zoom_region(s, s->cdim, dmin, ti_current, gravity_properties, verbose); + /* Allocate the highest level of cells. */ + s->tot_cells = s->nr_cells = + (s->cdim[0] * s->cdim[1] * s->cdim[2]) + + (s->zoom_props->cdim[0] * s->zoom_props->cdim[1] * + s->zoom_props->cdim[2]); + + if (swift_memalign("cells_top", (void **)&s->cells_top, cell_align, + s->nr_cells * sizeof(struct cell)) != 0) + error("Failed to allocate top-level cells."); + bzero(s->cells_top, s->nr_cells * sizeof(struct cell)); + + /* Allocate the multipoles for the top-level cells. */ + if (s->with_self_gravity) { + if (swift_memalign("multipoles_top", (void **)&s->multipoles_top, + multipole_align, + s->nr_cells * sizeof(struct gravity_tensors)) != 0) + error("Failed to allocate top-level multipoles."); + bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors)); + } + + /* Allocate the indices of local cells */ + if (swift_memalign("local_cells_top", (void **)&s->local_cells_top, + SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) + error("Failed to allocate indices of local top-level cells."); + bzero(s->local_cells_top, s->nr_cells * sizeof(int)); + + /* Allocate the indices of local zoom cells */ + if (swift_memalign("local_zoom_cells_top", + (void **)&s->zoom_props->local_zoom_cells_top, + SWIFT_STRUCT_ALIGNMENT, + s->zoom_props->nr_zoom_cells * sizeof(int)) != 0) + error("Failed to allocate indices of local top-level zoom cells."); + bzero(s->zoom_props->local_zoom_cells_top, + s->zoom_props->nr_zoom_cells * sizeof(int)); + + /* Allocate the indices of local bkg cells */ + if (swift_memalign("local_bkg_cells_top", + (void **)&s->zoom_props->local_bkg_cells_top, + SWIFT_STRUCT_ALIGNMENT, + s->zoom_props->nr_bkg_cells * sizeof(int)) != 0) + error("Failed to allocate indices of local top-level background cells."); + bzero(s->zoom_props->local_bkg_cells_top, + s->zoom_props->nr_bkg_cells * sizeof(int)); + + /* Allocate the indices of local zoom cells with particles */ + if (swift_memalign( + "local_zoom_cells_with_particles_top", + (void **)&s->zoom_props->local_zoom_cells_with_particles_top, + SWIFT_STRUCT_ALIGNMENT, + s->zoom_props->nr_zoom_cells * sizeof(int)) != 0) + error( + "Failed to allocate indices of local top-level zoom cells with " + "particles."); + bzero(s->zoom_props->local_zoom_cells_with_particles_top, + s->zoom_props->nr_zoom_cells * sizeof(int)); + + /* Allocate the indices of local bkg cells with particles */ + if (swift_memalign( + "local_bkg_cells_with_particles_top", + (void **)&s->zoom_props->local_bkg_cells_with_particles_top, + SWIFT_STRUCT_ALIGNMENT, + s->zoom_props->nr_bkg_cells * sizeof(int)) != 0) + error( + "Failed to allocate indices of local top-level background cells with " + "particles."); + bzero(s->zoom_props->local_bkg_cells_with_particles_top, + s->zoom_props->nr_bkg_cells * sizeof(int)); + + /* Allocate the indices of local cells with tasks */ + if (swift_memalign("local_cells_with_tasks_top", + (void **)&s->local_cells_with_tasks_top, + SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) + error("Failed to allocate indices of local top-level cells with tasks."); + bzero(s->local_cells_with_tasks_top, s->nr_cells * sizeof(int)); + + /* Allocate the indices of cells with particles */ + if (swift_memalign("cells_with_particles_top", + (void **)&s->cells_with_particles_top, + SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) + error("Failed to allocate indices of top-level cells with particles."); + bzero(s->cells_with_particles_top, s->nr_cells * sizeof(int)); + + /* Allocate the indices of local cells with particles */ + if (swift_memalign("local_cells_with_particles_top", + (void **)&s->local_cells_with_particles_top, + SWIFT_STRUCT_ALIGNMENT, s->nr_cells * sizeof(int)) != 0) + error( + "Failed to allocate indices of local top-level cells with " + "particles."); + bzero(s->local_cells_with_particles_top, s->nr_cells * sizeof(int)); + + /* Set the cells' locks */ + for (int k = 0; k < s->nr_cells; k++) { + if (lock_init(&s->cells_top[k].hydro.lock) != 0) + error("Failed to init spinlock for hydro."); + if (lock_init(&s->cells_top[k].grav.plock) != 0) + error("Failed to init spinlock for gravity."); + if (lock_init(&s->cells_top[k].grav.mlock) != 0) + error("Failed to init spinlock for multipoles."); + if (lock_init(&s->cells_top[k].grav.star_formation_lock) != 0) + error("Failed to init spinlock for star formation (gpart)."); + if (lock_init(&s->cells_top[k].stars.lock) != 0) + error("Failed to init spinlock for stars."); + if (lock_init(&s->cells_top[k].sinks.lock) != 0) + error("Failed to init spinlock for sinks."); + if (lock_init(&s->cells_top[k].sinks.sink_formation_lock) != 0) + error("Failed to init spinlock for sink formation."); + if (lock_init(&s->cells_top[k].black_holes.lock) != 0) + error("Failed to init spinlock for black holes."); + if (lock_init(&s->cells_top[k].stars.star_formation_lock) != 0) + error("Failed to init spinlock for star formation (spart)."); + } + + /* Construct both grids of cells */ + const double dmin = min3(s->width[0], s->width[1], s->width[2]); + construct_tl_cells_with_zoom_region(s, s->cdim, dmin, ti_current, + gravity_properties, verbose); #ifdef WITH_MPI - if (oldnodeIDs != NULL) { + if (oldnodeIDs != NULL) { /* We have changed the top-level cell dimension, so need to redistribute * cells around the nodes. We repartition using the old space node * positions as a grid to resample. */ @@ -352,7 +379,7 @@ void space_regrid_zoom(struct space *s, struct gravity_props *gravity_properties "global partition."); if (!partition_space_to_space_zoom(oldwidth, oldcdim, oldzoomwidth, - oldzoomcdim, oldnodeIDs, s)) { + oldzoomcdim, oldnodeIDs, s)) { /* Failed, try another technique that requires no settings. */ message("Failed to get a new partition, trying less optimal method"); @@ -392,18 +419,18 @@ void space_regrid_zoom(struct space *s, struct gravity_props *gravity_properties } #endif /* WITH_MPI */ - // message( "rebuilding upper-level cells took %.3f %s." , - // clocks_from_ticks(double)(getticks() - tic), clocks_getunit()); + // message( "rebuilding upper-level cells took %.3f %s." , + // clocks_from_ticks(double)(getticks() - tic), clocks_getunit()); - } /* re-build upper-level cells? */ - else { /* Otherwise, just clean up the cells. */ + } /* re-build upper-level cells? */ + else { /* Otherwise, just clean up the cells. */ - /* Free the old cells, if they were allocated. */ - space_free_cells(s); - } + /* Free the old cells, if they were allocated. */ + space_free_cells(s); + } - if (verbose) - message("took %.3f %s.", clocks_from_ticks(getticks() - tic), - clocks_getunit()); + if (verbose) + message("took %.3f %s.", clocks_from_ticks(getticks() - tic), + clocks_getunit()); } #endif /* WITH_ZOOM_REGION */ diff --git a/tests/testFFT.c b/tests/testFFT.c index 5668c92612..ab6d6c4260 100644 --- a/tests/testFFT.c +++ b/tests/testFFT.c @@ -83,8 +83,8 @@ int main(int argc, char *argv[]) { /* Build the infrastructure */ struct space space; double dim[3] = {1., 1., 1.}; - space_init(&space, params, &cosmo, dim, NULL, gparts, NULL, NULL, 0, nr_gparts, 0, - 1, 1, 0, 0, 1, 1, 0); + space_init(&space, params, &cosmo, dim, NULL, gparts, NULL, NULL, 0, + nr_gparts, 0, 1, 1, 0, 0, 1, 1, 0); struct engine engine; engine.s = &space; -- GitLab