Commit 43682ecb authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Also allow black hole particles when running over MPI.

parent 75e14a5d
......@@ -499,7 +499,6 @@ int main(int argc, char *argv[]) {
if (with_star_formation && with_feedback)
error("Can't run with star formation and feedback over MPI (yet)");
if (with_limiter) error("Can't run with time-step limiter over MPI (yet)");
if (with_black_holes) error("Can't run with black holes over MPI (yet)");
#endif
/* Temporary early aborts for modes not supported with hand-vec. */
......@@ -799,16 +798,18 @@ int main(int argc, char *argv[]) {
#if defined(HAVE_HDF5)
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
&Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas,
&Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
with_black_holes,
cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
dry_run);
#else
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
&Ngpart, &Nspart, &flag_entropy_ICs, with_hydro,
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas,
&Ngpart, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
with_black_holes,
cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
dry_run);
......@@ -834,6 +835,10 @@ int main(int argc, char *argv[]) {
for (size_t k = 0; k < Ngpart; ++k)
if (gparts[k].type == swift_type_stars) error("Linking problem");
}
if (!with_black_holes && !dry_run) {
for (size_t k = 0; k < Ngpart; ++k)
if (gparts[k].type == swift_type_black_hole) error("Linking problem");
}
if (!with_hydro && !dry_run) {
for (size_t k = 0; k < Ngpart; ++k)
if (gparts[k].type == swift_type_gas) error("Linking problem");
......@@ -1058,9 +1063,9 @@ int main(int argc, char *argv[]) {
/* Legend */
if (myrank == 0) {
printf("# %6s %14s %12s %12s %14s %9s %12s %12s %12s %16s [%s] %6s\n",
printf("# %6s %14s %12s %12s %14s %9s %12s %12s %12s %12s %16s [%s] %6s\n",
"Step", "Time", "Scale-factor", "Redshift", "Time-step", "Time-bins",
"Updates", "g-Updates", "s-Updates", "Wall-clock time",
"Updates", "g-Updates", "s-Updates", "b-Updates", "Wall-clock time",
clocks_getunit(), "Props");
fflush(stdout);
}
......
......@@ -188,6 +188,39 @@ int cell_link_sparts(struct cell *c, struct spart *sparts) {
return c->stars.count;
}
/**
* @brief Link the cells recursively to the given #bpart array.
*
* @param c The #cell.
* @param bparts The #bpart array.
*
* @return The number of particles linked.
*/
int cell_link_bparts(struct cell *c, struct bpart *bparts) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->nodeID == engine_rank)
error("Linking foreign particles in a local cell!");
if (c->black_holes.parts != NULL)
error("Linking bparts into a cell that was already linked");
#endif
c->black_holes.parts = bparts;
/* Fill the progeny recursively, depth-first. */
if (c->split) {
int offset = 0;
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL)
offset += cell_link_bparts(c->progeny[k], &bparts[offset]);
}
}
/* Return the total number of linked particles. */
return c->black_holes.count;
}
/**
* @brief Recurse down foreign cells until reaching one with hydro
* tasks; then trigger the linking of the #part array from that
......
......@@ -824,6 +824,7 @@ int cell_getsize(struct cell *c);
int cell_link_parts(struct cell *c, struct part *parts);
int cell_link_gparts(struct cell *c, struct gpart *gparts);
int cell_link_sparts(struct cell *c, struct spart *sparts);
int cell_link_bparts(struct cell *c, struct bpart *bparts);
int cell_link_foreign_parts(struct cell *c, struct part *parts);
int cell_link_foreign_gparts(struct cell *c, struct gpart *gparts);
int cell_count_parts_for_tasks(const struct cell *c);
......
......@@ -36,17 +36,20 @@
/* Local collections for MPI reduces. */
struct mpicollectgroup1 {
long long updated, g_updated, s_updated;
long long inhibited, g_inhibited, s_inhibited;
long long updated, g_updated, s_updated, b_updated;
long long inhibited, g_inhibited, s_inhibited, b_inhibited;
integertime_t ti_hydro_end_min;
integertime_t ti_gravity_end_min;
integertime_t ti_stars_end_min;
integertime_t ti_black_holes_end_min;
integertime_t ti_hydro_end_max;
integertime_t ti_gravity_end_max;
integertime_t ti_stars_end_max;
integertime_t ti_black_holes_end_max;
integertime_t ti_hydro_beg_max;
integertime_t ti_gravity_beg_max;
integertime_t ti_stars_beg_max;
integertime_t ti_black_holes_beg_max;
int forcerebuild;
long long total_nr_cells;
long long total_nr_tasks;
......@@ -96,18 +99,23 @@ void collectgroup1_apply(struct collectgroup1 *grp1, struct engine *e) {
e->ti_stars_end_min = grp1->ti_stars_end_min;
e->ti_stars_end_max = grp1->ti_stars_end_max;
e->ti_stars_beg_max = grp1->ti_stars_beg_max;
e->ti_black_holes_end_min = grp1->ti_black_holes_end_min;
e->ti_black_holes_end_max = grp1->ti_black_holes_end_max;
e->ti_black_holes_beg_max = grp1->ti_black_holes_beg_max;
e->ti_end_min =
min3(e->ti_hydro_end_min, e->ti_gravity_end_min, e->ti_stars_end_min);
min4(e->ti_hydro_end_min, e->ti_gravity_end_min, e->ti_stars_end_min, e->ti_black_holes_end_min);
e->ti_end_max =
max3(e->ti_hydro_end_max, e->ti_gravity_end_max, e->ti_stars_end_max);
max4(e->ti_hydro_end_max, e->ti_gravity_end_max, e->ti_stars_end_max , e->ti_black_holes_end_max);
e->ti_beg_max =
max3(e->ti_hydro_beg_max, e->ti_gravity_beg_max, e->ti_stars_beg_max);
max4(e->ti_hydro_beg_max, e->ti_gravity_beg_max, e->ti_stars_beg_max, e->ti_black_holes_beg_max);
e->updates = grp1->updated;
e->g_updates = grp1->g_updated;
e->s_updates = grp1->s_updated;
e->b_updates = grp1->b_updated;
e->nr_inhibited_parts = grp1->inhibited;
e->nr_inhibited_gparts = grp1->g_inhibited;
e->nr_inhibited_sparts = grp1->s_inhibited;
e->nr_inhibited_bparts = grp1->b_inhibited;
e->forcerebuild = grp1->forcerebuild;
e->total_nr_cells = grp1->total_nr_cells;
e->total_nr_tasks = grp1->total_nr_tasks;
......@@ -153,20 +161,24 @@ void collectgroup1_apply(struct collectgroup1 *grp1, struct engine *e) {
*/
void collectgroup1_init(
struct collectgroup1 *grp1, size_t updated, size_t g_updated,
size_t s_updated, size_t inhibited, size_t g_inhibited, size_t s_inhibited,
size_t s_updated, size_t b_updated, size_t inhibited, size_t g_inhibited, size_t s_inhibited,
size_t b_inhibited,
integertime_t ti_hydro_end_min, integertime_t ti_hydro_end_max,
integertime_t ti_hydro_beg_max, integertime_t ti_gravity_end_min,
integertime_t ti_gravity_end_max, integertime_t ti_gravity_beg_max,
integertime_t ti_stars_end_min, integertime_t ti_stars_end_max,
integertime_t ti_stars_beg_max, int forcerebuild, long long total_nr_cells,
integertime_t ti_stars_beg_max, integertime_t ti_black_holes_end_min, integertime_t ti_black_holes_end_max,
integertime_t ti_black_holes_beg_max, int forcerebuild, long long total_nr_cells,
long long total_nr_tasks, float tasks_per_cell) {
grp1->updated = updated;
grp1->g_updated = g_updated;
grp1->s_updated = s_updated;
grp1->b_updated = b_updated;
grp1->inhibited = inhibited;
grp1->g_inhibited = g_inhibited;
grp1->s_inhibited = s_inhibited;
grp1->b_inhibited = b_inhibited;
grp1->ti_hydro_end_min = ti_hydro_end_min;
grp1->ti_hydro_end_max = ti_hydro_end_max;
grp1->ti_hydro_beg_max = ti_hydro_beg_max;
......@@ -176,6 +188,9 @@ void collectgroup1_init(
grp1->ti_stars_end_min = ti_stars_end_min;
grp1->ti_stars_end_max = ti_stars_end_max;
grp1->ti_stars_beg_max = ti_stars_beg_max;
grp1->ti_black_holes_end_min = ti_black_holes_end_min;
grp1->ti_black_holes_end_max = ti_black_holes_end_max;
grp1->ti_black_holes_beg_max = ti_black_holes_beg_max;
grp1->forcerebuild = forcerebuild;
grp1->total_nr_cells = total_nr_cells;
grp1->total_nr_tasks = total_nr_tasks;
......@@ -199,18 +214,23 @@ void collectgroup1_reduce(struct collectgroup1 *grp1) {
mpigrp11.updated = grp1->updated;
mpigrp11.g_updated = grp1->g_updated;
mpigrp11.s_updated = grp1->s_updated;
mpigrp11.b_updated = grp1->b_updated;
mpigrp11.inhibited = grp1->inhibited;
mpigrp11.g_inhibited = grp1->g_inhibited;
mpigrp11.s_inhibited = grp1->s_inhibited;
mpigrp11.b_inhibited = grp1->b_inhibited;
mpigrp11.ti_hydro_end_min = grp1->ti_hydro_end_min;
mpigrp11.ti_gravity_end_min = grp1->ti_gravity_end_min;
mpigrp11.ti_stars_end_min = grp1->ti_stars_end_min;
mpigrp11.ti_black_holes_end_min = grp1->ti_black_holes_end_min;
mpigrp11.ti_hydro_end_max = grp1->ti_hydro_end_max;
mpigrp11.ti_gravity_end_max = grp1->ti_gravity_end_max;
mpigrp11.ti_stars_end_max = grp1->ti_stars_end_max;
mpigrp11.ti_black_holes_end_max = grp1->ti_black_holes_end_max;
mpigrp11.ti_hydro_beg_max = grp1->ti_hydro_beg_max;
mpigrp11.ti_gravity_beg_max = grp1->ti_gravity_beg_max;
mpigrp11.ti_stars_beg_max = grp1->ti_stars_beg_max;
mpigrp11.ti_black_holes_beg_max = grp1->ti_black_holes_beg_max;
mpigrp11.forcerebuild = grp1->forcerebuild;
mpigrp11.total_nr_cells = grp1->total_nr_cells;
mpigrp11.total_nr_tasks = grp1->total_nr_tasks;
......@@ -225,18 +245,23 @@ void collectgroup1_reduce(struct collectgroup1 *grp1) {
grp1->updated = mpigrp12.updated;
grp1->g_updated = mpigrp12.g_updated;
grp1->s_updated = mpigrp12.s_updated;
grp1->b_updated = mpigrp12.b_updated;
grp1->inhibited = mpigrp12.inhibited;
grp1->g_inhibited = mpigrp12.g_inhibited;
grp1->s_inhibited = mpigrp12.s_inhibited;
grp1->b_inhibited = mpigrp12.b_inhibited;
grp1->ti_hydro_end_min = mpigrp12.ti_hydro_end_min;
grp1->ti_gravity_end_min = mpigrp12.ti_gravity_end_min;
grp1->ti_stars_end_min = mpigrp12.ti_stars_end_min;
grp1->ti_black_holes_end_min = mpigrp12.ti_black_holes_end_min;
grp1->ti_hydro_end_max = mpigrp12.ti_hydro_end_max;
grp1->ti_gravity_end_max = mpigrp12.ti_gravity_end_max;
grp1->ti_stars_end_max = mpigrp12.ti_stars_end_max;
grp1->ti_black_holes_end_max = mpigrp12.ti_black_holes_end_max;
grp1->ti_hydro_beg_max = mpigrp12.ti_hydro_beg_max;
grp1->ti_gravity_beg_max = mpigrp12.ti_gravity_beg_max;
grp1->ti_stars_beg_max = mpigrp12.ti_stars_beg_max;
grp1->ti_black_holes_beg_max = mpigrp12.ti_black_holes_beg_max;
grp1->forcerebuild = mpigrp12.forcerebuild;
grp1->total_nr_cells = mpigrp12.total_nr_cells;
grp1->total_nr_tasks = mpigrp12.total_nr_tasks;
......@@ -260,11 +285,13 @@ static void doreduce1(struct mpicollectgroup1 *mpigrp11,
mpigrp11->updated += mpigrp12->updated;
mpigrp11->g_updated += mpigrp12->g_updated;
mpigrp11->s_updated += mpigrp12->s_updated;
mpigrp11->b_updated += mpigrp12->b_updated;
/* Sum of inhibited */
mpigrp11->inhibited += mpigrp12->inhibited;
mpigrp11->g_inhibited += mpigrp12->g_inhibited;
mpigrp11->s_inhibited += mpigrp12->s_inhibited;
mpigrp11->b_inhibited += mpigrp12->b_inhibited;
/* Minimum end time. */
mpigrp11->ti_hydro_end_min =
......@@ -273,6 +300,8 @@ static void doreduce1(struct mpicollectgroup1 *mpigrp11,
min(mpigrp11->ti_gravity_end_min, mpigrp12->ti_gravity_end_min);
mpigrp11->ti_stars_end_min =
min(mpigrp11->ti_stars_end_min, mpigrp12->ti_stars_end_min);
mpigrp11->ti_black_holes_end_min =
min(mpigrp11->ti_black_holes_end_min, mpigrp12->ti_black_holes_end_min);
/* Maximum end time. */
mpigrp11->ti_hydro_end_max =
......@@ -281,6 +310,8 @@ static void doreduce1(struct mpicollectgroup1 *mpigrp11,
max(mpigrp11->ti_gravity_end_max, mpigrp12->ti_gravity_end_max);
mpigrp11->ti_stars_end_max =
max(mpigrp11->ti_stars_end_max, mpigrp12->ti_stars_end_max);
mpigrp11->ti_black_holes_end_max =
max(mpigrp11->ti_black_holes_end_max, mpigrp12->ti_black_holes_end_max);
/* Maximum beg time. */
mpigrp11->ti_hydro_beg_max =
......@@ -289,6 +320,8 @@ static void doreduce1(struct mpicollectgroup1 *mpigrp11,
max(mpigrp11->ti_gravity_beg_max, mpigrp12->ti_gravity_beg_max);
mpigrp11->ti_stars_beg_max =
max(mpigrp11->ti_stars_beg_max, mpigrp12->ti_stars_beg_max);
mpigrp11->ti_black_holes_beg_max =
max(mpigrp11->ti_black_holes_beg_max, mpigrp12->ti_black_holes_beg_max);
/* Everyone must agree to not rebuild. */
if (mpigrp11->forcerebuild || mpigrp12->forcerebuild)
......
......@@ -35,15 +35,16 @@ struct engine;
struct collectgroup1 {
/* Number of particles updated */
long long updated, g_updated, s_updated;
long long updated, g_updated, s_updated, b_updated;
/* Number of particles inhibited */
long long inhibited, g_inhibited, s_inhibited;
long long inhibited, g_inhibited, s_inhibited, b_inhibited;
/* Times for the time-step */
integertime_t ti_hydro_end_min, ti_hydro_end_max, ti_hydro_beg_max;
integertime_t ti_gravity_end_min, ti_gravity_end_max, ti_gravity_beg_max;
integertime_t ti_stars_end_min, ti_stars_end_max, ti_stars_beg_max;
integertime_t ti_black_holes_end_min, ti_black_holes_end_max, ti_black_holes_beg_max;
/* Force the engine to rebuild? */
int forcerebuild;
......@@ -59,14 +60,16 @@ struct collectgroup1 {
void collectgroup_init(void);
void collectgroup1_apply(struct collectgroup1 *grp1, struct engine *e);
void collectgroup1_init(
struct collectgroup1 *grp1, size_t updated, size_t g_updated,
size_t s_updated, size_t inhibited, size_t g_inhibited, size_t s_inhibited,
struct collectgroup1 *grp1, size_t updated, size_t g_updated,
size_t s_updated, size_t b_updated, size_t inhibited, size_t g_inhibited, size_t s_inhibited,
size_t b_inhibited,
integertime_t ti_hydro_end_min, integertime_t ti_hydro_end_max,
integertime_t ti_hydro_beg_max, integertime_t ti_gravity_end_min,
integertime_t ti_gravity_end_max, integertime_t ti_gravity_beg_max,
integertime_t ti_stars_end_min, integertime_t ti_stars_end_max,
integertime_t ti_stars_beg_max, int forcerebuild, long long total_nr_cells,
long long total_nr_tasks, float tasks_per_cell);
integertime_t ti_stars_beg_max, integertime_t ti_black_holes_end_min, integertime_t ti_black_holes_end_max,
integertime_t ti_black_holes_beg_max, int forcerebuild, long long total_nr_cells,
long long total_nr_tasks, float tasks_per_cell);
void collectgroup1_reduce(struct collectgroup1 *grp1);
#endif /* SWIFT_COLLECTGROUP_H */
This diff is collapsed.
......@@ -204,6 +204,15 @@ struct engine {
/* Maximal stars ti_beg for the next time-step */
integertime_t ti_stars_beg_max;
/* Minimal black holes ti_end for the next time-step */
integertime_t ti_black_holes_end_min;
/* Maximal black holes ti_end for the next time-step */
integertime_t ti_black_holes_end_max;
/* Maximal black holes ti_beg for the next time-step */
integertime_t ti_black_holes_beg_max;
/* Minimal overall ti_end for the next time-step */
integertime_t ti_end_min;
......@@ -473,7 +482,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
const int *ind_part, size_t *Npart,
const size_t offset_gparts, const int *ind_gpart,
size_t *Ngpart, const size_t offset_sparts,
const int *ind_spart, size_t *Nspart);
const int *ind_spart, size_t *Nspart, const size_t offset_bparts,
const int *ind_bpart, size_t *Nbpart);
void engine_rebuild(struct engine *e, int redistributed, int clean_h_values);
void engine_repartition(struct engine *e);
void engine_repartition_trigger(struct engine *e);
......
......@@ -37,6 +37,7 @@
#include "parallel_io.h"
/* Local includes. */
#include "black_holes_io.h"
#include "chemistry_io.h"
#include "common_io.h"
#include "cooling_io.h"
......@@ -651,14 +652,17 @@ void writeArray(struct engine* e, hid_t grp, char* fileName,
* @param parts (output) The array of #part read from the file.
* @param gparts (output) The array of #gpart read from the file.
* @param sparts (output) The array of #spart read from the file.
* @param bparts (output) The array of #bpart read from the file.
* @param Ngas (output) The number of particles read from the file.
* @param Ngparts (output) The number of particles read from the file.
* @param Nstars (output) The number of particles read from the file.
* @param Nblackholes (output) The number of particles read from the file.
* @param flag_entropy (output) 1 if the ICs contained Entropy in the
* InternalEnergy field
* @param with_hydro Are we running with hydro ?
* @param with_gravity Are we running with gravity ?
* @param with_stars Are we running with stars ?
* @param with_black_holes Are we running with black holes ?
* @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
* @param cleanup_sqrt_a Are we cleaning-up the sqrt(a) factors in the Gadget
* IC velocities?
......@@ -674,9 +678,9 @@ void writeArray(struct engine* e, hid_t grp, char* fileName,
*/
void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, int* flag_entropy, int with_hydro,
int with_gravity, int with_stars, int cleanup_h,
struct spart** sparts, struct bpart** bparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, size_t* Nblackholes, int* flag_entropy, int with_hydro,
int with_gravity, int with_stars, int with_black_holes, int cleanup_h,
int cleanup_sqrt_a, double h, double a, int mpi_rank,
int mpi_size, MPI_Comm comm, MPI_Info info, int n_threads,
int dry_run) {
......@@ -692,6 +696,9 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
int dimension = 3; /* Assume 3D if nothing is specified */
size_t Ndm = 0;
/* Initialise counters */
*Ngas = 0, *Ngparts = 0, *Nstars = 0, *Nblackholes = 0;
/* Open file */
/* message("Opening file '%s' as IC.", fileName); */
hid_t h_plist_id = H5Pcreate(H5P_FILE_ACCESS);
......@@ -831,12 +838,22 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
bzero(*sparts, *Nstars * sizeof(struct spart));
}
/* Allocate memory to store black hole particles */
if (with_black_holes) {
*Nblackholes = N[swift_type_black_hole];
if (swift_memalign("sparts", (void**)bparts, bpart_align,
*Nblackholes * sizeof(struct bpart)) != 0)
error("Error while allocating memory for black_holes particles");
bzero(*bparts, *Nblackholes * sizeof(struct bpart));
}
/* Allocate memory to store gravity particles */
if (with_gravity) {
Ndm = N[1];
*Ngparts = (with_hydro ? N[swift_type_gas] : 0) +
N[swift_type_dark_matter] +
(with_stars ? N[swift_type_stars] : 0);
(with_stars ? N[swift_type_stars] : 0) +
(with_black_holes ? N[swift_type_black_hole] : 0);
if (swift_memalign("gparts", (void**)gparts, gpart_align,
*Ngparts * sizeof(struct gpart)) != 0)
error("Error while allocating memory for gravity particles");
......@@ -892,6 +909,13 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
}
break;
case swift_type_black_hole:
if (with_black_holes) {
Nparticles = *Nblackholes;
black_holes_read_particles(*bparts, list, &num_fields);
}
break;
default:
if (mpi_rank == 0)
message("Particle Type %d not yet supported. Particles ignored",
......@@ -925,6 +949,10 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
if (with_stars)
io_duplicate_stars_gparts(&tp, *sparts, *gparts, *Nstars, Ndm + *Ngas);
/* Duplicate the stars particles into gparts */
if (with_black_holes)
io_duplicate_black_holes_gparts(&tp, *bparts, *gparts, *Nblackholes, Ndm + *Ngas + *Nstars);
threadpool_clean(&tp);
}
......@@ -957,6 +985,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
const struct xpart* xparts = e->s->xparts;
const struct gpart* gparts = e->s->gparts;
const struct spart* sparts = e->s->sparts;
const struct bpart* bparts = e->s->bparts;
struct swift_params* params = e->parameter_file;
const int with_cosmology = e->policy & engine_policy_cosmology;
const int with_cooling = e->policy & engine_policy_cooling;
......@@ -1182,6 +1211,13 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
}
break;
case swift_type_black_hole:
black_holes_write_particles(bparts, list, &num_fields);
if (with_stf) {
num_fields += velociraptor_write_bparts(bparts, list + num_fields);
}
break;
default:
error("Particle Type %d not yet supported. Aborting", ptype);
}
......@@ -1245,6 +1281,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
const struct xpart* xparts = e->s->xparts;
const struct gpart* gparts = e->s->gparts;
const struct spart* sparts = e->s->sparts;
const struct bpart* bparts = e->s->bparts;
struct swift_params* params = e->parameter_file;
const int with_cosmology = e->policy & engine_policy_cosmology;
const int with_cooling = e->policy & engine_policy_cooling;
......@@ -1260,6 +1297,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
const size_t Ntot = e->s->nr_gparts;
const size_t Ngas = e->s->nr_parts;
const size_t Nstars = e->s->nr_sparts;
const size_t Nblackholes = e->s->nr_bparts;
// const size_t Nbaryons = Ngas + Nstars;
// const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
......@@ -1270,6 +1308,8 @@ void write_output_parallel(struct engine* e, const char* baseName,
e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
const size_t Nstars_written =
e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
const size_t Nblackholes_written =
e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
const size_t Nbaryons_written = Ngas_written + Nstars_written;
const size_t Ndm_written =
Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
......@@ -1458,6 +1498,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
struct gpart* gparts_written = NULL;
struct velociraptor_gpart_data* gpart_group_data_written = NULL;
struct spart* sparts_written = NULL;
struct bpart* bparts_written = NULL;
/* Write particle fields from the particle structure */
switch (ptype) {
......@@ -1608,6 +1649,39 @@ void write_output_parallel(struct engine* e, const char* baseName,
}
} break;
case swift_type_black_hole: {
if (Nblackholes == Nblackholes_written) {
/* No inhibted particles: easy case */
Nparticles = Nblackholes;
black_holes_write_particles(bparts, list, &num_fields);
if (with_stf) {
num_fields += velociraptor_write_bparts(bparts, list + num_fields);
}
} else {
/* Ok, we need to fish out the particles we want */
Nparticles = Nblackholes_written;
/* Allocate temporary arrays */
if (swift_memalign("bparts_written", (void**)&bparts_written,
bpart_align,
Nblackholes_written * sizeof(struct bpart)) != 0)
error("Error while allocating temporart memory for bparts");
/* Collect the particles we want to write */
io_collect_bparts_to_write(bparts, bparts_written, Nblackholes,
Nblackholes_written);
/* Select the fields to write */
black_holes_write_particles(bparts_written, list, &num_fields);
if (with_stf) {
num_fields +=
velociraptor_write_bparts(bparts_written, list + num_fields);
}
}
} break;
default:
error("Particle Type %d not yet supported. Aborting", ptype);
}
......@@ -1634,6 +1708,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
if (gpart_group_data_written)
swift_free("gpart_group_written", gpart_group_data_written);
if (sparts_written) swift_free("sparts_written", sparts_written);
if (bparts_written) swift_free("bparts_written", bparts_written);
#ifdef IO_SPEED_MEASUREMENT
MPI_Barrier(MPI_COMM_WORLD);
......
......@@ -35,9 +35,10 @@
void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nsparts, int* flag_entropy, int with_hydro,
int with_gravity, int with_stars, int cleanup_h,
struct spart** sparts, struct bpart** bparts, size_t* Ngas,
size_t* Ngparts,
size_t* Nsparts, size_t *Nbparts, int* flag_entropy, int with_hydro,
int with_gravity, int with_stars, int with_black_holes, int cleanup_h,
int cleanup_sqrt_a, double h, double a, int mpi_rank,
int mpi_size, MPI_Comm comm, MPI_Info info,
int nr_threads, int dry_run);
......@@ -47,13 +48,6 @@ void write_output_parallel(struct engine* e, const char* baseName,
const struct unit_system* snapshot_units,
int mpi_rank, int mpi_size, MPI_Comm comm,
MPI_Info info);
void writeArray(struct engine* e, hid_t grp, char* fileName,
char* partTypeGroupName, struct io_props props, size_t N,
long long N_total, int mpi_rank, long long offset,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units);
#endif
#endif /* SWIFT_PARALLEL_IO_H */
......@@ -64,6 +64,22 @@ void part_relink_gparts_to_sparts(struct spart *sparts, size_t N,
}
}
/**
* @brief Re-link the #gpart%s associated with the list of #bpart%s.
*
* @param bparts The list of #bpart.
* @param N The number of s-particles to re-link;
* @param offset The offset of #bpart%s relative to the global bparts list.
*/
void part_relink_gparts_to_bparts(struct bpart *bparts, size_t N,
ptrdiff_t offset) {
for (size_t k = 0; k < N; k++) {
if (bparts[k].gpart) {
bparts[k].gpart->id_or_neg_offset = -(k + offset);
}
}
}
/**
* @brief Re-link the #part%s associated with the list of #gpart%s.
*
......@@ -243,7 +259,7 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
gparts[k].x[2] != bpart->x[2])
error(
"Linked particles are not at the same position !\n"
"gp->x=[%e %e %e] sp->x=[%e %e %e] diff=[%e %e %e]",
"gp->x=[%e %e %e] bp->x=[%e %e %e] diff=[%e %e %e]",
gparts[k].x[0], gparts[k].x[1], gparts[k].x[2], bpart->x[0],
bpart->x[1], bpart->x[2], gparts[k].x[0] - bpart->x[0],
gparts[k].x[1] - bpart->x[1], gparts[k].x[2] - bpart->x[2]);
......
......@@ -115,6 +115,8 @@ void part_relink_gparts_to_parts(struct part *parts, size_t N,
ptrdiff_t offset);
void part_relink_gparts_to_sparts(struct spart *sparts, size_t N,
ptrdiff_t offset);
void part_relink_gparts_to_bparts(struct bpart *bparts, size_t N,
ptrdiff_t offset);
void part_relink_parts_to_gparts(struct gpart *gparts, size_t N,
struct part *parts);
void part_relink_sparts_to_gparts(struct gpart *gparts, size_t N,
......
......@@ -584,7 +584,8 @@ void proxy_parts_exchange_first(struct proxy *p) {
p->buff_out[0] = p->nr_parts_out;