diff --git a/doc/RTD/source/Snapshots/index.rst b/doc/RTD/source/Snapshots/index.rst
index f3e81c6410da14f859eb8843e91f7f45889ddfd1..ae82ab48c0da639ef96ed803f6828b94373c5e37 100644
--- a/doc/RTD/source/Snapshots/index.rst
+++ b/doc/RTD/source/Snapshots/index.rst
@@ -24,10 +24,18 @@ that is always set to the string ``SWIFT``, which can be used to identify
SWIFT-generated snapshots and hence make use of all the extensions to the file
format described below.
-The most important quantity of the header is the array ``NumPart_ThisFile``
-which contains the number of particles of each type in this snapshot. This is an
-array of 6 numbers; one for each of the 5 supported types and a dummy "type 3"
-field only used for compatibility reasons but always containing a zero.
+The most important quantity of the header is the array ``NumPart_Total`` which
+contains the number of particles of each type in this snapshot. This is an array
+of 6 numbers; one for each of the supported types. The field
+``NumPart_ThisFile`` contains the number of particles in this sub-snapshot file
+when the user asked for distributed snapshots (see :ref:`Parameters_snapshots`);
+otherwise it contains the same information as ``NumPart_Total``. The field
+``NumFilesPerSnapshot`` specifies the number of sub-snapshot files (always 1
+unless a distributed snapshot was asked).
+
+The field ``InitialMassTable`` contains the *mean* initial mass of each of the
+particle types present in the initial conditions. This can be used as estimator
+of the mass resolution of the run. The masses are expressed in internal units.
The ``RunName`` field contains the name of the simulation that was specified as
the ``run_name`` in the :ref:`Parameters_meta_data` section of the YAML
diff --git a/src/distributed_io.c b/src/distributed_io.c
index f5cf8b3d3ce64a2c7a085fca544b31cb129cdccf..3d5a550a62b9afef0fdc27b3ebc1d3c3dbab803d 100644
--- a/src/distributed_io.c
+++ b/src/distributed_io.c
@@ -434,6 +434,8 @@ void write_output_distributed(struct engine* e,
numParticlesHighWord, swift_type_count);
double MassTable[swift_type_count] = {0};
io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
+ io_write_attribute(h_grp, "InitialMassTable", DOUBLE,
+ e->s->initial_mean_mass_particles, swift_type_count);
unsigned int flagEntropy[swift_type_count] = {0};
flagEntropy[0] = writeEntropyFlag();
io_write_attribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
diff --git a/src/engine.c b/src/engine.c
index d0de2005b58101a1cbd5200b2f4f277b330d70a1..40e69458ce5eecd08f59f78c1388eee8dcd06203 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1771,6 +1771,9 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
}
}
+ /* Collect initial mean mass of each particle type */
+ space_collect_mean_masses(e->s, e->verbose);
+
#ifdef SWIFT_DEBUG_CHECKS
/* Check that we have the correct total mass in the top-level multipoles */
long long num_gpart_mpole = 0;
diff --git a/src/line_of_sight.c b/src/line_of_sight.c
index 3c686f99280f951a5d466d9b172c49e7c7c21a89..dec2e7c947b335c2a741a55a37be01d6bf1695be 100644
--- a/src/line_of_sight.c
+++ b/src/line_of_sight.c
@@ -531,6 +531,14 @@ void write_hdf5_header(hid_t h_file, const struct engine *e,
swift_type_count);
io_write_attribute(h_grp, "NumPart_Total_HighWord", UINT,
numParticlesHighWord, swift_type_count);
+ double MassTable[swift_type_count] = {0};
+ io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
+ io_write_attribute(h_grp, "InitialMassTable", DOUBLE,
+ e->s->initial_mean_mass_particles, swift_type_count);
+ unsigned int flagEntropy[swift_type_count] = {0};
+ flagEntropy[0] = writeEntropyFlag();
+ io_write_attribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
+ swift_type_count);
io_write_attribute_i(h_grp, "NumFilesPerSnapshot", 1);
io_write_attribute_i(h_grp, "ThisFile", 0);
io_write_attribute_s(h_grp, "OutputType", "LineOfSight");
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 55b345912da6b57060bf587179089afcf3c88f9a..de756d3e2fa35f4fb7e65260a3f2afaf049c723d 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1208,6 +1208,8 @@ void prepare_file(struct engine* e, const char* fileName,
numParticlesHighWord, swift_type_count);
double MassTable[6] = {0., 0., 0., 0., 0., 0.};
io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
+ io_write_attribute(h_grp, "InitialMassTable", DOUBLE,
+ e->s->initial_mean_mass_particles, swift_type_count);
unsigned int flagEntropy[swift_type_count] = {0};
flagEntropy[0] = writeEntropyFlag();
io_write_attribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
diff --git a/src/serial_io.c b/src/serial_io.c
index eb801c8f10a21f1336e6dc13897f57925a0905bb..c76b0c23ea777beba7500ec523838a87ba522426 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -1094,6 +1094,8 @@ void write_output_serial(struct engine* e,
numParticlesHighWord, swift_type_count);
double MassTable[6] = {0., 0., 0., 0., 0., 0.};
io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
+ io_write_attribute(h_grp, "InitialMassTable", DOUBLE,
+ e->s->initial_mean_mass_particles, swift_type_count);
unsigned int flagEntropy[swift_type_count] = {0};
flagEntropy[0] = writeEntropyFlag();
io_write_attribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
diff --git a/src/single_io.c b/src/single_io.c
index e0912ba45e53945a3ee98d9d945343df77292b53..b46e331e6cc506d5a3faa895f2518f2d3a94f9fc 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -933,6 +933,8 @@ void write_output_single(struct engine* e,
numParticlesHighWord, swift_type_count);
double MassTable[swift_type_count] = {0};
io_write_attribute(h_grp, "MassTable", DOUBLE, MassTable, swift_type_count);
+ io_write_attribute(h_grp, "InitialMassTable", DOUBLE,
+ e->s->initial_mean_mass_particles, swift_type_count);
unsigned int flagEntropy[swift_type_count] = {0};
flagEntropy[0] = writeEntropyFlag();
io_write_attribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
diff --git a/src/space.c b/src/space.c
index 9c8afa141aed1d638c94da3fea4ac5d390074cfd..6747524f4fc0c33d5eb9b0392ba88cb8a85a78ae 100644
--- a/src/space.c
+++ b/src/space.c
@@ -842,6 +842,134 @@ void space_convert_quantities(struct space *s, int verbose) {
clocks_getunit());
}
+void space_collect_sum_part_mass(void *restrict map_data, int count,
+ void *restrict extra_data) {
+
+ struct space *s = (struct space *)extra_data;
+ const struct part *parts = (const struct part *)map_data;
+
+ /* Local collection */
+ double sum = 0.;
+ for (int i = 0; i < count; ++i) sum += hydro_get_mass(&parts[i]);
+
+ /* Store back */
+ atomic_add_d(&s->initial_mean_mass_particles[0], sum);
+ atomic_add(&s->initial_count_particles[0], count);
+}
+
+void space_collect_sum_gpart_mass(void *restrict map_data, int count,
+ void *restrict extra_data) {
+
+ struct space *s = (struct space *)extra_data;
+ const struct gpart *gparts = (const struct gpart *)map_data;
+
+ /* Local collection */
+ double sum_DM = 0., sum_DM_background = 0.;
+ long long count_DM = 0, count_DM_background = 0;
+ for (int i = 0; i < count; ++i) {
+ if (gparts[i].type == swift_type_dark_matter) {
+ sum_DM += gparts[i].mass;
+ count_DM++;
+ }
+ if (gparts[i].type == swift_type_dark_matter_background) {
+ sum_DM_background += gparts[i].mass;
+ count_DM_background++;
+ }
+ }
+
+ /* Store back */
+ atomic_add_d(&s->initial_mean_mass_particles[1], sum_DM);
+ atomic_add(&s->initial_count_particles[1], count_DM);
+ atomic_add_d(&s->initial_mean_mass_particles[2], sum_DM_background);
+ atomic_add(&s->initial_count_particles[2], count_DM_background);
+}
+
+void space_collect_sum_sink_mass(void *restrict map_data, int count,
+ void *restrict extra_data) {
+
+ struct space *s = (struct space *)extra_data;
+ const struct sink *sinks = (const struct sink *)map_data;
+
+ /* Local collection */
+ double sum = 0.;
+ for (int i = 0; i < count; ++i) sum += sinks[i].mass;
+
+ /* Store back */
+ atomic_add_d(&s->initial_mean_mass_particles[3], sum);
+ atomic_add(&s->initial_count_particles[3], count);
+}
+
+void space_collect_sum_spart_mass(void *restrict map_data, int count,
+ void *restrict extra_data) {
+
+ struct space *s = (struct space *)extra_data;
+ const struct spart *sparts = (const struct spart *)map_data;
+
+ /* Local collection */
+ double sum = 0.;
+ for (int i = 0; i < count; ++i) sum += sparts[i].mass;
+
+ /* Store back */
+ atomic_add_d(&s->initial_mean_mass_particles[4], sum);
+ atomic_add(&s->initial_count_particles[4], count);
+}
+
+void space_collect_sum_bpart_mass(void *restrict map_data, int count,
+ void *restrict extra_data) {
+
+ struct space *s = (struct space *)extra_data;
+ const struct bpart *bparts = (const struct bpart *)map_data;
+
+ /* Local collection */
+ double sum = 0.;
+ for (int i = 0; i < count; ++i) sum += bparts[i].mass;
+
+ /* Store back */
+ atomic_add_d(&s->initial_mean_mass_particles[5], sum);
+ atomic_add(&s->initial_count_particles[5], count);
+}
+
+/**
+ * @breif Collect the mean mass of each particle type in the #space.
+ */
+void space_collect_mean_masses(struct space *s, int verbose) {
+
+ /* Init counters */
+ for (int i = 0; i < swift_type_count; ++i)
+ s->initial_mean_mass_particles[i] = 0.;
+ for (int i = 0; i < swift_type_count; ++i) s->initial_count_particles[i] = 0;
+
+ /* Collect each particle type */
+ threadpool_map(&s->e->threadpool, space_collect_sum_part_mass, s->parts,
+ s->nr_parts, sizeof(struct part), threadpool_auto_chunk_size,
+ s);
+ threadpool_map(&s->e->threadpool, space_collect_sum_gpart_mass, s->gparts,
+ s->nr_gparts, sizeof(struct gpart), threadpool_auto_chunk_size,
+ s);
+ threadpool_map(&s->e->threadpool, space_collect_sum_spart_mass, s->sparts,
+ s->nr_sparts, sizeof(struct spart), threadpool_auto_chunk_size,
+ s);
+ threadpool_map(&s->e->threadpool, space_collect_sum_sink_mass, s->sinks,
+ s->nr_sinks, sizeof(struct sink), threadpool_auto_chunk_size,
+ s);
+ threadpool_map(&s->e->threadpool, space_collect_sum_bpart_mass, s->bparts,
+ s->nr_bparts, sizeof(struct bpart), threadpool_auto_chunk_size,
+ s);
+
+#ifdef WITH_MPI
+ MPI_Allreduce(MPI_IN_PLACE, s->initial_mean_mass_particles, swift_type_count,
+ MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, s->initial_count_particles, swift_type_count,
+ MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
+#endif
+
+ /* Get means */
+ for (int i = 0; i < swift_type_count; ++i)
+ if (s->initial_count_particles[i] > 0)
+ s->initial_mean_mass_particles[i] /=
+ (double)s->initial_count_particles[i];
+}
+
/**
* @brief Split the space into cells given the array of particles.
*
diff --git a/src/space.h b/src/space.h
index 17949392f46f63ff8d980a5fcc8da1236dad1c45..c726a68d9873f7473c191f649fb3d7fa53b10400 100644
--- a/src/space.h
+++ b/src/space.h
@@ -276,6 +276,12 @@ struct space {
/*! Sum of the norm of the velocity of all the #bpart */
float sum_bpart_vel_norm;
+ /* Initial mean mass of each particle type in the system. */
+ double initial_mean_mass_particles[swift_type_count];
+
+ /* Initial count of each particle type in the system. */
+ long long initial_count_particles[swift_type_count];
+
/*! Initial value of the smoothing length read from the parameter file */
float initial_spart_h;
@@ -385,6 +391,7 @@ void space_first_init_gparts(struct space *s, int verbose);
void space_first_init_sparts(struct space *s, int verbose);
void space_first_init_bparts(struct space *s, int verbose);
void space_first_init_sinks(struct space *s, int verbose);
+void space_collect_mean_masses(struct space *s, int verbose);
void space_init_parts(struct space *s, int verbose);
void space_init_gparts(struct space *s, int verbose);
void space_init_sparts(struct space *s, int verbose);