diff --git a/doc/RTD/source/Snapshots/index.rst b/doc/RTD/source/Snapshots/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..86284d520838bcfaf081add9fca782f208cf0bf3
--- /dev/null
+++ b/doc/RTD/source/Snapshots/index.rst
@@ -0,0 +1,120 @@
+.. Snapshots
+ Matthieu Schaller, 5th January 2019
+
+.. _snapshots:
+
+Snapshots
+=========
+
+The snapshots are stored using the HDF5 format and are fully compatible with
+Gadget-2. They do, however, contain a large set of extensions including units,
+meta-data about the code and runs as well as facilities to quickly access the
+particles in a specific region of the simulation volume.
+
+Meta-data about the code and run
+--------------------------------
+
+Unit system
+-----------
+
+Used and unused run-time parameters
+-----------------------------------
+
+Structure of the particle arrays
+--------------------------------
+
+There are several groups that contain 'auxiliary' information, such as
+``Header``. Particle data is placed in separate groups depending of the type of
+the particles. The type use the naming convention of Gadget-2 (with
+the OWLS and EAGLE extensions).
+
++---------------------+------------------------+----------------------------+
+| HDF5 Group Name | Physical Particle Type | In code ``enum part_type`` |
++=====================+========================+============================+
+| ``/PartType0/`` | Gas | ``swift_type_gas`` |
++---------------------+------------------------+----------------------------+
+| ``/PartType1/`` | Dark Matter | ``swift_type_dark_matter`` |
++---------------------+------------------------+----------------------------+
+| ``/PartType4/`` | Stars | ``swift_type_star`` |
++---------------------+------------------------+----------------------------+
+| ``/PartType5/`` | Black Holes | ``swift_type_black_hole`` |
++---------------------+------------------------+----------------------------+
+
+The last column in the table gives the ``enum`` value from ``part_type.h``
+corresponding to a given entry in the files.
+
+Quick access to particles via hash-tables
+-----------------------------------------
+
+The particles are not sorted in a specific order when they are written to the
+snapshots. However, the particles are sorted into the top-level cell structure
+used internally by the code every time a tree rebuild is triggered. The
+top-level cells are a coarse-grained mesh but knowing which particle belongs to
+which cell can nevertheless be useful to rapidly access particles in a given
+region only.
+
+One important caveat is that particles are free to drift out of their cells
+between rebuilds of the tree (but not by more than one cell-length). If one
+wants to have all the particles in a given cell, one has to read all the
+neighbouring cells as well. We note that for image making purposes, for instance
+to generate a slice, this is typically not necessary and reading just the cells
+of interest is sufficient.
+
+At the root of the HDF5 file, the ``Cells`` group contains all the relevant
+information. The dimension of the top-level grid (a triplet of integers) is
+given by the attribute ``Cells/Meta-data/dimension`` and the size of each cell (a
+triplet of floating-point numbers) is given by the attribute
+``Cells/Meta-data/size``. All the cells have the same size but for non-cubic
+simulation volumes the cells themselves can have different sizes along each
+axis.
+
+The ``/Cells/Centres`` array gives the centre of each of the top-level cells in the
+simulation volume. Both the cell sizes and positions of the centres are
+expressed in the unit system used for the snapshots (see above) and are hence
+consistent with the particle positions themselves.
+
+Once the cell(s) containing the region of interest has been located, users can
+use the ``/Cells/Offsets/PartTypeN/Counts`` and
+``/Cells/Offsets/PartTypeN/Offsets`` to retrieve the location of the particles
+of type ``N`` in the ``/PartTypeN`` arrays. For instance, if one is interested
+in retriving all the densities of the gas particles in the cell around the
+position `[1, 1, 1]` one could use a piece of code similar to:
+
+.. code-block:: python
+ :linenos:
+
+ import numpy as np
+ import h5py
+
+ snapshot_file = h5py.File("snapshot.hdf5", "r")
+
+ my_pos = [1, 1, 1]
+
+ # Read in the cell centres and size
+ nr_cells = f["/Cells/Meta-data"].attrs["nr_cells"]
+ centres = f["/Cells/Centres"][:,:]
+ size = f["/Cells/Meta-data"].attrs["size"]
+ half_size = size / 2.
+
+ # Look for the cell containing the position of interest
+ my_cell = -1
+ for i in range(nr_cells):
+ if my_pos[0] > centres[i, 0] - half_size[0] and my_pos[0] < centres[i, 0] + half_size[0] and
+ my_pos[1] > centres[i, 1] - half_size[1] and my_pos[1] < centres[i, 1] + half_size[1] and
+ my_pos[2] > centres[i, 2] - half_size[2] and my_pos[2] < centres[i, 2] + half_size[2]:
+ my_cell = i
+ break
+
+ # Print the position of the centre of the cell of interest
+ centre = snapshot_file["/Cells/Centres"][my_cell, :]
+ print("Centre of the cell:", centre)
+
+ # Retrieve the offset and counts
+ my_offset = snapshot_file["/Cells/Offsets/PartType0"][my_cell]
+ my_count = snapshot_file["/Cells/Counts/PartType0"][my_cell]
+
+ # Get the densities of the particles in this cell
+ rho = snapshot_file["/PartType0/Density"][my_offset:my_offset + my_count]
+
+For large simulations, this vastly reduces the amount of data that needs to be read
+from the disk.
diff --git a/doc/RTD/source/index.rst b/doc/RTD/source/index.rst
index b9370c3f24b2ffb3c5174f2fe99fb9ec610e18f6..e04efe8c889fb8a005c88f691f1e01a387f19ebb 100644
--- a/doc/RTD/source/index.rst
+++ b/doc/RTD/source/index.rst
@@ -18,6 +18,7 @@ difference is the parameter file that will need to be adapted for SWIFT.
CommandLineOptions/index
ParameterFiles/index
InitialConditions/index
+ Snapshots/index
HydroSchemes/index
SubgridModels/index
EquationOfState/index
diff --git a/src/common_io.c b/src/common_io.c
index 24e74014fd52936023b5c7a41378faf3268bfdb3..6cce7e21aa9d0eb0e8841050e4b5f6d612417a12 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -142,7 +142,7 @@ void io_read_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
* Calls #error() if an error occurs.
*/
void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
- void* data, int num) {
+ const void* data, int num) {
const hid_t h_space = H5Screate(H5S_SIMPLE);
if (h_space < 0)
@@ -387,6 +387,332 @@ void io_write_engine_policy(hid_t h_file, const struct engine* e) {
H5Gclose(h_grp);
}
+void io_write_cell_offsets(hid_t h_grp, const int cdim[3],
+ const struct cell* cells_top, const int nr_cells,
+ const double width[3], const int nodeID,
+ const long long global_counts[swift_type_count],
+ const long long global_offsets[swift_type_count],
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units) {
+
+ double cell_width[3] = {width[0], width[1], width[2]};
+
+ /* Temporary memory for the cell-by-cell information */
+ double* centres = NULL;
+ centres = malloc(3 * nr_cells * sizeof(double));
+
+ /* Count of particles in each cell */
+ long long *count_part = NULL, *count_gpart = NULL, *count_spart = NULL;
+ count_part = malloc(nr_cells * sizeof(long long));
+ count_gpart = malloc(nr_cells * sizeof(long long));
+ count_spart = malloc(nr_cells * sizeof(long long));
+
+ /* Global offsets of particles in each cell */
+ long long *offset_part = NULL, *offset_gpart = NULL, *offset_spart = NULL;
+ offset_part = malloc(nr_cells * sizeof(long long));
+ offset_gpart = malloc(nr_cells * sizeof(long long));
+ offset_spart = malloc(nr_cells * sizeof(long long));
+
+ /* Offsets of the 0^th element */
+ offset_part[0] = 0;
+ offset_gpart[0] = 0;
+ offset_spart[0] = 0;
+
+ /* Collect the cell information of *local* cells */
+ long long local_offset_part = 0;
+ long long local_offset_gpart = 0;
+ long long local_offset_spart = 0;
+ for (int i = 0; i < nr_cells; ++i) {
+
+ if (cells_top[i].nodeID == nodeID) {
+
+ /* Centre of each cell */
+ centres[i * 3 + 0] = cells_top[i].loc[0] + cell_width[0] * 0.5;
+ centres[i * 3 + 1] = cells_top[i].loc[1] + cell_width[1] * 0.5;
+ centres[i * 3 + 2] = cells_top[i].loc[2] + cell_width[2] * 0.5;
+
+ /* Count real particles that will be written */
+ count_part[i] = cells_top[i].hydro.count - cells_top[i].hydro.inhibited;
+ count_gpart[i] = cells_top[i].grav.count - cells_top[i].grav.inhibited;
+ count_spart[i] = cells_top[i].stars.count - cells_top[i].stars.inhibited;
+
+ /* Only count DM gpart (gpart without friends) */
+ count_gpart[i] -= count_part[i];
+ count_gpart[i] -= count_spart[i];
+
+ /* Offsets including the global offset of all particles on this MPI rank
+ */
+ offset_part[i] = local_offset_part + global_offsets[swift_type_gas];
+ offset_gpart[i] =
+ local_offset_gpart + global_offsets[swift_type_dark_matter];
+ offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
+
+ local_offset_part += count_part[i];
+ local_offset_gpart += count_gpart[i];
+ local_offset_spart += count_spart[i];
+
+ } else {
+
+ /* Just zero everything for the foregin cells */
+
+ centres[i * 3 + 0] = 0.;
+ centres[i * 3 + 1] = 0.;
+ centres[i * 3 + 2] = 0.;
+
+ count_part[i] = 0;
+ count_gpart[i] = 0;
+ count_spart[i] = 0;
+
+ offset_part[i] = 0;
+ offset_gpart[i] = 0;
+ offset_spart[i] = 0;
+ }
+ }
+
+#ifdef WITH_MPI
+ /* Now, reduce all the arrays. Note that we use a bit-wise OR here. This
+ is safe as we made sure only local cells have non-zero values. */
+ if (nodeID == 0) {
+ MPI_Reduce(MPI_IN_PLACE, count_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ 0, MPI_COMM_WORLD);
+ } else {
+ MPI_Reduce(count_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
+ MPI_COMM_WORLD);
+ }
+ if (nodeID == 0) {
+ MPI_Reduce(MPI_IN_PLACE, count_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ 0, MPI_COMM_WORLD);
+ } else {
+ MPI_Reduce(count_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
+ MPI_COMM_WORLD);
+ }
+ if (nodeID == 0) {
+ MPI_Reduce(MPI_IN_PLACE, count_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ 0, MPI_COMM_WORLD);
+ } else {
+ MPI_Reduce(count_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
+ MPI_COMM_WORLD);
+ }
+ if (nodeID == 0) {
+ MPI_Reduce(MPI_IN_PLACE, offset_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ 0, MPI_COMM_WORLD);
+ } else {
+ MPI_Reduce(offset_part, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
+ MPI_COMM_WORLD);
+ }
+ if (nodeID == 0) {
+ MPI_Reduce(MPI_IN_PLACE, offset_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ 0, MPI_COMM_WORLD);
+ } else {
+ MPI_Reduce(offset_gpart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
+ MPI_COMM_WORLD);
+ }
+ if (nodeID == 0) {
+ MPI_Reduce(MPI_IN_PLACE, offset_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ 0, MPI_COMM_WORLD);
+ } else {
+ MPI_Reduce(offset_spart, NULL, nr_cells, MPI_LONG_LONG_INT, MPI_BOR, 0,
+ MPI_COMM_WORLD);
+ }
+
+ /* For the centres we use a sum as MPI does not like bit-wise operations
+ on floating point numbers */
+ if (nodeID == 0) {
+ MPI_Reduce(MPI_IN_PLACE, centres, 3 * nr_cells, MPI_DOUBLE, MPI_SUM, 0,
+ MPI_COMM_WORLD);
+ } else {
+ MPI_Reduce(centres, NULL, 3 * nr_cells, MPI_DOUBLE, MPI_SUM, 0,
+ MPI_COMM_WORLD);
+ }
+#endif
+
+ /* Only rank 0 actually writes */
+ if (nodeID == 0) {
+
+ /* Unit conversion if necessary */
+ const double factor = units_conversion_factor(
+ internal_units, snapshot_units, UNIT_CONV_LENGTH);
+ if (factor != 1.) {
+
+ /* Convert the cell centres */
+ for (int i = 0; i < nr_cells; ++i) {
+ centres[i * 3 + 0] *= factor;
+ centres[i * 3 + 1] *= factor;
+ centres[i * 3 + 2] *= factor;
+ }
+
+ /* Convert the cell widths */
+ cell_width[0] *= factor;
+ cell_width[1] *= factor;
+ cell_width[2] *= factor;
+ }
+
+ /* Write some meta-information first */
+ hid_t h_subgrp =
+ H5Gcreate(h_grp, "Meta-data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_subgrp < 0) error("Error while creating meta-data sub-group");
+ io_write_attribute(h_subgrp, "nr_cells", INT, &nr_cells, 1);
+ io_write_attribute(h_subgrp, "size", DOUBLE, cell_width, 3);
+ io_write_attribute(h_subgrp, "dimension", INT, cdim, 3);
+ H5Gclose(h_subgrp);
+
+ /* Write the centres to the group */
+ hsize_t shape[2] = {nr_cells, 3};
+ hid_t h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0) error("Error while creating data space for cell centres");
+ hid_t h_err = H5Sset_extent_simple(h_space, 2, shape, shape);
+ if (h_err < 0)
+ error("Error while changing shape of gas offsets data space.");
+ hid_t h_data = H5Dcreate(h_grp, "Centres", io_hdf5_type(DOUBLE), h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace for gas offsets.");
+ h_err = H5Dwrite(h_data, io_hdf5_type(DOUBLE), h_space, H5S_ALL,
+ H5P_DEFAULT, centres);
+ if (h_err < 0) error("Error while writing centres.");
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+
+ /* Group containing the offsets for each particle type */
+ h_subgrp =
+ H5Gcreate(h_grp, "Offsets", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_subgrp < 0) error("Error while creating offsets sub-group");
+
+ if (global_counts[swift_type_gas] > 0) {
+
+ shape[0] = nr_cells;
+ shape[1] = 1;
+ h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0) error("Error while creating data space for gas offsets");
+ h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
+ if (h_err < 0)
+ error("Error while changing shape of gas offsets data space.");
+ h_data = H5Dcreate(h_subgrp, "PartType0", io_hdf5_type(LONGLONG), h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace for gas offsets.");
+ h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
+ H5P_DEFAULT, offset_part);
+ if (h_err < 0) error("Error while writing gas offsets.");
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+ }
+
+ if (global_counts[swift_type_dark_matter] > 0) {
+
+ shape[0] = nr_cells;
+ shape[1] = 1;
+ h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0) error("Error while creating data space for DM offsets");
+ h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
+ if (h_err < 0)
+ error("Error while changing shape of DM offsets data space.");
+ h_data = H5Dcreate(h_subgrp, "PartType1", io_hdf5_type(LONGLONG), h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace for DM offsets.");
+ h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
+ H5P_DEFAULT, offset_gpart);
+ if (h_err < 0) error("Error while writing DM offsets.");
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+ }
+
+ if (global_counts[swift_type_stars] > 0) {
+
+ shape[0] = nr_cells;
+ shape[1] = 1;
+ h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0)
+ error("Error while creating data space for stars offsets");
+ h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
+ if (h_err < 0)
+ error("Error while changing shape of stars offsets data space.");
+ h_data = H5Dcreate(h_subgrp, "PartType4", io_hdf5_type(LONGLONG), h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace for star offsets.");
+ h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
+ H5P_DEFAULT, offset_spart);
+ if (h_err < 0) error("Error while writing star offsets.");
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+ }
+
+ H5Gclose(h_subgrp);
+
+ /* Group containing the counts for each particle type */
+ h_subgrp =
+ H5Gcreate(h_grp, "Counts", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_subgrp < 0) error("Error while creating counts sub-group");
+
+ if (global_counts[swift_type_gas] > 0) {
+
+ shape[0] = nr_cells;
+ shape[1] = 1;
+ h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0) error("Error while creating data space for gas counts");
+ h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
+ if (h_err < 0)
+ error("Error while changing shape of gas counts data space.");
+ h_data = H5Dcreate(h_subgrp, "PartType0", io_hdf5_type(LONGLONG), h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace for gas counts.");
+ h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
+ H5P_DEFAULT, count_part);
+ if (h_err < 0) error("Error while writing gas counts.");
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+ }
+
+ if (global_counts[swift_type_dark_matter] > 0) {
+
+ shape[0] = nr_cells;
+ shape[1] = 1;
+ h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0) error("Error while creating data space for DM counts");
+ h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
+ if (h_err < 0)
+ error("Error while changing shape of DM counts data space.");
+ h_data = H5Dcreate(h_subgrp, "PartType1", io_hdf5_type(LONGLONG), h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace for DM counts.");
+ h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
+ H5P_DEFAULT, count_gpart);
+ if (h_err < 0) error("Error while writing DM counts.");
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+ }
+
+ if (global_counts[swift_type_stars] > 0) {
+
+ shape[0] = nr_cells;
+ shape[1] = 1;
+ h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0)
+ error("Error while creating data space for stars counts");
+ h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
+ if (h_err < 0)
+ error("Error while changing shape of stars counts data space.");
+ h_data = H5Dcreate(h_subgrp, "PartType4", io_hdf5_type(LONGLONG), h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace for star counts.");
+ h_err = H5Dwrite(h_data, io_hdf5_type(LONGLONG), h_space, H5S_ALL,
+ H5P_DEFAULT, count_spart);
+ if (h_err < 0) error("Error while writing star counts.");
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+ }
+
+ H5Gclose(h_subgrp);
+ }
+
+ /* Free everything we allocated */
+ free(centres);
+ free(count_part);
+ free(count_gpart);
+ free(count_spart);
+ free(offset_part);
+ free(offset_gpart);
+ free(offset_spart);
+}
+
#endif /* HAVE_HDF5 */
/**
diff --git a/src/common_io.h b/src/common_io.h
index 016c5138e18ae8636834c35d659e07d8fcd46e36..bdff3e37d11c59a71c08dcf3e19758019adc3a73 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -24,6 +24,7 @@
#include "../config.h"
/* Local includes. */
+#include "part_type.h"
#include "units.h"
#define FIELD_BUFFER_SIZE 200
@@ -32,6 +33,7 @@
#define IO_BUFFER_ALIGNMENT 1024
/* Avoid cyclic inclusion problems */
+struct cell;
struct part;
struct gpart;
struct spart;
@@ -65,7 +67,7 @@ void io_read_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
void* data);
void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
- void* data, int num);
+ const void* data, int num);
void io_write_attribute_d(hid_t grp, const char* name, double data);
void io_write_attribute_f(hid_t grp, const char* name, float data);
@@ -76,6 +78,14 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str);
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_cell_offsets(hid_t h_grp, const int cdim[3],
+ const struct cell* cells_top, const int nr_cells,
+ const double width[3], const int nodeID,
+ const long long global_counts[swift_type_count],
+ const long long global_offsets[swift_type_count],
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units);
+
void io_read_unit_system(hid_t h_file, struct unit_system* ic_units,
const struct unit_system* internal_units,
int mpi_rank);
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index fffbf22ec187f179f0e80b7121beaa3a96de0260..e548e3010f3b46065a2510723b5bde97121b4c02 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -170,20 +170,22 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
io_write_attribute_s(h_grpgrav, "Softening style",
kernel_gravity_softening_name);
io_write_attribute_f(
- h_grpgrav, "Comoving softening length",
+ h_grpgrav, "Comoving softening length [internal units]",
p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent);
- io_write_attribute_f(h_grpgrav,
- "Comoving Softening length (Plummer equivalent)",
- p->epsilon_comoving);
io_write_attribute_f(
- h_grpgrav, "Maximal physical softening length",
+ h_grpgrav,
+ "Comoving Softening length (Plummer equivalent) [internal units]",
+ p->epsilon_comoving);
+ io_write_attribute_f(
+ h_grpgrav, "Maximal physical softening length [internal units]",
p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent);
io_write_attribute_f(h_grpgrav,
- "Maximal physical softening length (Plummer equivalent)",
+ "Maximal physical softening length (Plummer equivalent) "
+ " [internal units]",
p->epsilon_max_physical);
io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
io_write_attribute_s(h_grpgrav, "Scheme", GRAVITY_IMPLEMENTATION);
- io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
+ io_write_attribute_i(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
io_write_attribute_f(h_grpgrav, "Mesh r_cut_max ratio", p->r_cut_max_ratio);
io_write_attribute_f(h_grpgrav, "Mesh r_cut_min ratio", p->r_cut_min_ratio);
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 85f88d418bd46354f7a1cd3dd89b0e77b556b7d9..d864eb00367a76c2287f53dbded0fb6feb60ef26 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -239,7 +239,8 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
io_write_attribute_f(h_grpsph, "Kernel delta N_ngb", p->delta_neighbours);
io_write_attribute_f(h_grpsph, "Kernel eta", p->eta_neighbours);
io_write_attribute_f(h_grpsph, "Smoothing length tolerance", p->h_tolerance);
- io_write_attribute_f(h_grpsph, "Maximal smoothing length", p->h_max);
+ io_write_attribute_f(h_grpsph, "Maximal smoothing length [internal units]",
+ p->h_max);
io_write_attribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
io_write_attribute_f(h_grpsph, "Volume log(max(delta h))",
p->log_max_h_change);
@@ -248,8 +249,12 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
io_write_attribute_i(h_grpsph, "Max ghost iterations",
p->max_smoothing_iterations);
io_write_attribute_f(h_grpsph, "Minimal temperature", p->minimal_temperature);
+ io_write_attribute_f(h_grpsph,
+ "Minimal energy per unit mass [internal units]",
+ p->minimal_internal_energy);
io_write_attribute_f(h_grpsph, "Initial temperature", p->initial_temperature);
- io_write_attribute_f(h_grpsph, "Initial energy per unit mass",
+ io_write_attribute_f(h_grpsph,
+ "Initial energy per unit mass [internal units]",
p->initial_internal_energy);
io_write_attribute_f(h_grpsph, "Hydrogen mass fraction",
p->hydrogen_mass_fraction);
@@ -260,7 +265,8 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
p->viscosity.alpha_max);
io_write_attribute_f(h_grpsph, "Alpha viscosity (min)",
p->viscosity.alpha_min);
- io_write_attribute_f(h_grpsph, "Viscosity decay length", p->viscosity.length);
+ io_write_attribute_f(h_grpsph, "Viscosity decay length [internal units]",
+ p->viscosity.length);
io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
}
#endif
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 611e4af99f138a557719b505f0891c345025e41e..b62422d22fb099e0b4da31a0bd1063b952268575 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -957,7 +957,6 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
const int with_cosmology = e->policy & engine_policy_cosmology;
FILE* xmfFile = 0;
- int periodic = e->s->periodic;
int numFiles = 1;
/* First time, we need to create the XMF file */
@@ -983,28 +982,25 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
* specific output */
xmf_write_outputheader(xmfFile, fileName, e->time);
- /* Open header to write simulation properties */
- /* message("Writing runtime parameters..."); */
- hid_t h_grp =
- H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_grp < 0) error("Error while creating runtime parameters group\n");
-
- /* Write the relevant information */
- io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
-
- /* Close runtime parameters */
- H5Gclose(h_grp);
-
/* Open header to write simulation properties */
/* message("Writing file header..."); */
h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (h_grp < 0) error("Error while creating file header\n");
+ /* Convert basic output information to snapshot units */
+ const double factor_time =
+ units_conversion_factor(internal_units, snapshot_units, UNIT_CONV_TIME);
+ const double factor_length =
+ units_conversion_factor(internal_units, snapshot_units, UNIT_CONV_LENGTH);
+ const double dblTime = e->time * factor_time;
+ const double dim[3] = {e->s->dim[0] * factor_length,
+ e->s->dim[1] * factor_length,
+ e->s->dim[2] * factor_length};
+
/* Print the relevant information and print status */
- io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
- double dblTime = e->time;
+ io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
- int dimension = (int)hydro_dimension;
+ const int dimension = (int)hydro_dimension;
io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
@@ -1282,6 +1278,32 @@ void write_output_parallel(struct engine* e, const char* baseName,
snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%04i.hdf5", baseName,
e->snapshot_output_count);
+ /* Now write the top-level cell structure */
+ hid_t h_file_cells = 0, h_grp_cells = 0;
+ if (mpi_rank == 0) {
+
+ /* Open the snapshot on rank 0 */
+ h_file_cells = H5Fopen(fileName, H5F_ACC_RDWR, H5P_DEFAULT);
+ if (h_file_cells < 0)
+ error("Error while opening file '%s' on rank %d.", fileName, mpi_rank);
+
+ /* Create the group we want in the file */
+ h_grp_cells = H5Gcreate(h_file_cells, "/Cells", H5P_DEFAULT, H5P_DEFAULT,
+ H5P_DEFAULT);
+ if (h_grp_cells < 0) error("Error while creating cells group");
+ }
+
+ /* Write the location of the particles in the arrays */
+ io_write_cell_offsets(h_grp_cells, e->s->cdim, e->s->cells_top,
+ e->s->nr_cells, e->s->width, mpi_rank, N_total, offset,
+ internal_units, snapshot_units);
+
+ /* Close everything */
+ if (mpi_rank == 0) {
+ H5Gclose(h_grp_cells);
+ H5Fclose(h_file_cells);
+ }
+
/* Prepare some file-access properties */
hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
diff --git a/src/serial_io.c b/src/serial_io.c
index a273d033323702ad22a0ac63532f5f81765d5fca..c05c43a6521d46f5145d5bc29ad544c82e39ff04 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -776,7 +776,6 @@ void write_output_serial(struct engine* e, const char* baseName,
int mpi_size, MPI_Comm comm, MPI_Info info) {
hid_t h_file = 0, h_grp = 0;
- int periodic = e->s->periodic;
int numFiles = 1;
const struct part* parts = e->s->parts;
const struct xpart* xparts = e->s->xparts;
@@ -845,28 +844,25 @@ void write_output_serial(struct engine* e, const char* baseName,
h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
if (h_file < 0) error("Error while opening file '%s'.", fileName);
- /* Open header to write simulation properties */
- /* message("Writing runtime parameters..."); */
- h_grp = H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT,
- H5P_DEFAULT);
- if (h_grp < 0) error("Error while creating runtime parameters group\n");
-
- /* Write the relevant information */
- io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
-
- /* Close runtime parameters */
- H5Gclose(h_grp);
-
/* Open header to write simulation properties */
/* message("Writing file header..."); */
h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (h_grp < 0) error("Error while creating file header\n");
+ /* Convert basic output information to snapshot units */
+ const double factor_time =
+ units_conversion_factor(internal_units, snapshot_units, UNIT_CONV_TIME);
+ const double factor_length = units_conversion_factor(
+ internal_units, snapshot_units, UNIT_CONV_LENGTH);
+ const double dblTime = e->time * factor_time;
+ const double dim[3] = {e->s->dim[0] * factor_length,
+ e->s->dim[1] * factor_length,
+ e->s->dim[2] * factor_length};
+
/* Print the relevant information and print status */
- io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
- double dblTime = e->time;
+ io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
- int dimension = (int)hydro_dimension;
+ const int dimension = (int)hydro_dimension;
io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
@@ -1028,6 +1024,32 @@ void write_output_serial(struct engine* e, const char* baseName,
H5Fclose(h_file);
}
+ /* Now write the top-level cell structure */
+ hid_t h_file_cells = 0, h_grp_cells = 0;
+ if (mpi_rank == 0) {
+
+ /* Open the snapshot on rank 0 */
+ h_file_cells = H5Fopen(fileName, H5F_ACC_RDWR, H5P_DEFAULT);
+ if (h_file_cells < 0)
+ error("Error while opening file '%s' on rank %d.", fileName, mpi_rank);
+
+ /* Create the group we want in the file */
+ h_grp_cells = H5Gcreate(h_file_cells, "/Cells", H5P_DEFAULT, H5P_DEFAULT,
+ H5P_DEFAULT);
+ if (h_grp_cells < 0) error("Error while creating cells group");
+ }
+
+ /* Write the location of the particles in the arrays */
+ io_write_cell_offsets(h_grp_cells, e->s->cdim, e->s->cells_top,
+ e->s->nr_cells, e->s->width, mpi_rank, N_total, offset,
+ internal_units, snapshot_units);
+
+ /* Close everything */
+ if (mpi_rank == 0) {
+ H5Gclose(h_grp_cells);
+ H5Fclose(h_file_cells);
+ }
+
/* Now loop over ranks and write the data */
for (int rank = 0; rank < mpi_size; ++rank) {
diff --git a/src/single_io.c b/src/single_io.c
index 93d5f7f0e5f200668a3e4ffe0e53b5d37b024cf6..178f59a528016ddf6e642ac4e81deda72f508f85 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -639,7 +639,6 @@ void write_output_single(struct engine* e, const char* baseName,
const struct unit_system* snapshot_units) {
hid_t h_file = 0, h_grp = 0;
- int periodic = e->s->periodic;
int numFiles = 1;
const struct part* parts = e->s->parts;
const struct xpart* xparts = e->s->xparts;
@@ -698,28 +697,25 @@ void write_output_single(struct engine* e, const char* baseName,
h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
if (h_file < 0) error("Error while opening file '%s'.", fileName);
- /* Open header to write simulation properties */
- /* message("Writing runtime parameters..."); */
- h_grp =
- H5Gcreate(h_file, "/RuntimePars", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_grp < 0) error("Error while creating runtime parameters group\n");
-
- /* Write the relevant information */
- io_write_attribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
-
- /* Close runtime parameters */
- H5Gclose(h_grp);
-
/* Open header to write simulation properties */
/* message("Writing file header..."); */
h_grp = H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (h_grp < 0) error("Error while creating file header\n");
+ /* Convert basic output information to snapshot units */
+ const double factor_time =
+ units_conversion_factor(internal_units, snapshot_units, UNIT_CONV_TIME);
+ const double factor_length =
+ units_conversion_factor(internal_units, snapshot_units, UNIT_CONV_LENGTH);
+ const double dblTime = e->time * factor_time;
+ const double dim[3] = {e->s->dim[0] * factor_length,
+ e->s->dim[1] * factor_length,
+ e->s->dim[2] * factor_length};
+
/* Print the relevant information and print status */
- io_write_attribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
- double dblTime = e->time;
+ io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
- int dimension = (int)hydro_dimension;
+ const int dimension = (int)hydro_dimension;
io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
@@ -826,6 +822,17 @@ void write_output_single(struct engine* e, const char* baseName,
/* Print the system of Units used internally */
io_write_unit_system(h_file, internal_units, "InternalCodeUnits");
+ /* Now write the top-level cell structure */
+ long long global_offsets[swift_type_count] = {0};
+ 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 */
+ io_write_cell_offsets(h_grp, e->s->cdim, e->s->cells_top, e->s->nr_cells,
+ e->s->width, e->nodeID, N_total, global_offsets,
+ internal_units, snapshot_units);
+ H5Gclose(h_grp);
+
/* Tell the user if a conversion will be needed */
if (e->verbose) {
if (units_are_equal(snapshot_units, internal_units)) {