diff --git a/src/serial_io.c b/src/serial_io.c
index 2e4b9076dbd3677f2bd4c7c353890226759b16ab..fbf967da81b7d63668ba438ece3e9f6da03cb8a7 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -321,6 +321,64 @@ void read_ic_serial ( char* fileName, double dim[3], struct part **parts, int*
* Routines writing an output file
*-----------------------------------------------------------------------------*/
+void prepareArray(hid_t grp, char* fileName, FILE* xmfFile, char* name, enum DATA_TYPE type, long long N_total, int dim, struct UnitSystem* us, enum UnitConversionFactor convFactor)
+{
+ hid_t h_data=0, h_err=0, h_space=0;
+ void* temp = 0;
+ int i=0, rank=0;
+ const size_t typeSize = sizeOfType(type);
+ const size_t copySize = typeSize * dim;
+ const size_t partSize = sizeof(struct part);
+ char* temp_c = 0;
+ hsize_t shape[2];
+ char buffer[150];
+
+ /* Create data space */
+ h_space = H5Screate(H5S_SIMPLE);
+ if(h_space < 0)
+ {
+ error( "Error while creating data space for field '%s'." , name );
+ }
+
+ if(dim > 1)
+ {
+ rank = 2;
+ shape[0] = N_total; shape[1] = dim;
+ }
+ else
+ {
+ rank = 1;
+ shape[0] = N_total; shape[1] = 0;
+ }
+
+ /* Change shape of data space */
+ h_err = H5Sset_extent_simple(h_space, rank, shape, NULL);
+ if(h_err < 0)
+ {
+ error( "Error while changing data space shape for field '%s'." , name );
+ }
+
+ /* Create dataset */
+ h_data = H5Dcreate1(grp, name, hdf5Type(type), h_space, H5P_DEFAULT);
+ if(h_data < 0)
+ {
+ error( "Error while creating dataspace '%s'." , name );
+ }
+
+ /* Write XMF description for this data set */
+ writeXMFline(xmfFile, fileName, name, N_total, dim, type);
+
+ /* Write unit conversion factors for this data set */
+ conversionString( buffer, us, convFactor );
+ writeAttribute_d( h_data, "CGS conversion factor", conversionFactor( us, convFactor ) );
+ writeAttribute_f( h_data, "h-scale exponant", hFactor( us, convFactor ) );
+ writeAttribute_f( h_data, "a-scale exponant", aFactor( us, convFactor ) );
+ writeAttribute_s( h_data, "Conversion factor", buffer );
+
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+}
+
/**
* @brief Writes a data array in given HDF5 group.
*
@@ -335,22 +393,19 @@ void read_ic_serial ( char* fileName, double dim[3], struct part **parts, int*
* @param us The UnitSystem currently in use
* @param convFactor The UnitConversionFactor for this array
*
- * @todo A better version using HDF5 hyperslabs to write the file directly from the part array
- * will be written once the strucutres have been stabilized.
*
* Calls #error() if an error occurs.
*/
-void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, enum DATA_TYPE type, int N, int dim, char* part_c, struct UnitSystem* us, enum UnitConversionFactor convFactor)
+void writeArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N, int dim, long long N_total, long long offset, char* part_c)
{
- hid_t h_data=0, h_err=0, h_space=0;
+ hid_t h_data=0, h_err=0, h_memspace=0, h_filespace=0;
+ hsize_t shape[2], shape_total[2], offsets[2];
void* temp = 0;
int i=0, rank=0;
const size_t typeSize = sizeOfType(type);
const size_t copySize = typeSize * dim;
const size_t partSize = sizeof(struct part);
char* temp_c = 0;
- hsize_t shape[2];
- char buffer[150];
/* message("Writing '%s' array...", name); */
@@ -364,59 +419,52 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, enu
for(i=0; i<N; ++i)
memcpy(&temp_c[i*copySize], part_c+i*partSize, copySize);
- /* Create data space */
- h_space = H5Screate(H5S_SIMPLE);
- if(h_space < 0)
- {
- error( "Error while creating data space for field '%s'." , name );
- }
-
+ /* Construct information for the hyperslab */
if(dim > 1)
{
rank = 2;
shape[0] = N; shape[1] = dim;
+ shape_total[0] = N_total; shape_total[1] = dim;
+ offsets[0] = offset; offsets[1] = 0;
}
else
{
rank = 1;
shape[0] = N; shape[1] = 0;
+ shape_total[0] = N_total; shape_total[1] = 0;
+ offsets[0] = offset; offsets[1] = 0;
}
+
- /* Change shape of data space */
- h_err = H5Sset_extent_simple(h_space, rank, shape, NULL);
+ /* Create data space in memory */
+ h_memspace = H5Screate(H5S_SIMPLE);
+ if(h_memspace < 0)
+ error( "Error while creating data space (memory) for field '%s'." , name );
+
+ /* Change shape of memory data space */
+ h_err = H5Sset_extent_simple(h_memspace, rank, shape, NULL);
if(h_err < 0)
- {
- error( "Error while changing data space shape for field '%s'." , name );
- }
+ error( "Error while changing data space (memory) shape for field '%s'." , name );
- /* Create dataset */
- h_data = H5Dcreate1(grp, name, hdf5Type(type), h_space, H5P_DEFAULT);
+ /* Open pre-existing data set */
+ h_data = H5Dopen(grp, name, H5P_DEFAULT);
if(h_data < 0)
- {
- error( "Error while creating dataspace '%s'." , name );
- }
-
- /* Write temporary buffer to HDF5 dataspace */
- h_err = H5Dwrite(h_data, hdf5Type(type), h_space, H5S_ALL, H5P_DEFAULT, temp);
- if(h_err < 0)
- {
- error( "Error while writing data array '%s'." , name );
- }
+ error( "Error while opening dataset '%s'." , name );
- /* Write XMF description for this data set */
- writeXMFline(xmfFile, fileName, name, N, dim, type);
+ /* Select data space in that data set */
+ h_filespace = H5Dget_space(h_data);
+ H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape, NULL);
- /* Write unit conversion factors for this data set */
- conversionString( buffer, us, convFactor );
- writeAttribute_d( h_data, "CGS conversion factor", conversionFactor( us, convFactor ) );
- writeAttribute_f( h_data, "h-scale exponant", hFactor( us, convFactor ) );
- writeAttribute_f( h_data, "a-scale exponant", aFactor( us, convFactor ) );
- writeAttribute_s( h_data, "Conversion factor", buffer );
-
+ /* Write temporary buffer to HDF5 dataspace */
+ h_err = H5Dwrite(h_data, hdf5Type(type), h_memspace, h_filespace, H5P_DEFAULT, temp);
+ if(h_err < 0)
+ error( "Error while writing data array '%s'." , name );
+
/* Free and close everything */
free(temp);
H5Dclose(h_data);
- H5Sclose(h_space);
+ H5Sclose(h_memspace);
+ H5Sclose(h_filespace);
}
/**
@@ -435,7 +483,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, enu
* @param convFactor The UnitConversionFactor for this array
*
*/
-#define writeArray(grp, fileName, xmfFile, name, type, N, dim, part, field, us, convFactor) writeArrayBackEnd(grp, fileName, xmfFile, name, type, N, dim, (char*)(&(part[0]).field), us, convFactor)
+#define writeArray(grp, name, type, N, dim, N_total, offset, part, field) writeArrayBackEnd(grp, name, type, N, dim, N_total, offset, (char*)(&(part[0]).field))
+
/**
* @brief Writes an HDF5 output file (GADGET-3 type) with its XMF descriptor
@@ -452,14 +501,17 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, enu
*
*/
void write_output_serial ( struct engine* e, struct UnitSystem* us, int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info )
-{
-
+{
hid_t h_file=0, h_grp=0;
int N = e->s->nr_parts;
int periodic = e->s->periodic;
int numParticles[6]={N,0};
int numParticlesHighWord[6]={0};
+ unsigned int flagEntropy[6]={0};
+ long long N_total = 0, offset = 0;
+ double offset_d = 0., N_d = 0., N_total_d = 0.;
int numFiles = 1;
+ int rank = 0;
struct part* parts = e->s->parts;
FILE* xmfFile = 0;
static int outputCount = 0;
@@ -468,93 +520,151 @@ void write_output_serial ( struct engine* e, struct UnitSystem* us, int mpi_rank
char fileName[200];
sprintf(fileName, "output_%03i.hdf5", outputCount);
- /* First time, we need to create the XMF file */
- if(outputCount == 0)
- createXMFfile();
-
- /* Prepare the XMF file for the new entry */
- xmfFile = prepareXMFfile();
+ /* Compute offset in the file and total number of particles */
+ /* Done using double to allow for up to 2^50=10^15 particles */
+ N_d = (double)N;
+ MPI_Exscan(&N_d, &offset_d, 1, MPI_DOUBLE, MPI_SUM, comm);
+ N_total_d = offset_d + N_d;
+ MPI_Bcast(&N_total_d, 1, MPI_DOUBLE, mpi_size-1, comm);
+ if(N_total_d > 1.e15)
+ error("Error while computing the offest for parallel output: Simulation has more than 10^15 particles.\n");
+ N_total = (long long) N_total_d;
+ offset = (long long) offset_d;
- /* Write the part corresponding to this specific output */
- writeXMFheader(xmfFile, N, fileName, e->time);
+ /* Do common stuff first */
+ if ( mpi_rank == 0 ) {
- /* Open file */
- /* message("Opening file '%s'.", fileName); */
- h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT,H5P_DEFAULT);
- if(h_file < 0)
- {
- error( "Error while opening file '%s'." , fileName );
- }
+ /* First time, we need to create the XMF file */
+ if(outputCount == 0)
+ createXMFfile();
+
+ /* Prepare the XMF file for the new entry */
+ xmfFile = prepareXMFfile();
+
+ /* Write the part corresponding to this specific output */
+ writeXMFheader(xmfFile, N_total, fileName, e->time);
- /* Open header to write simulation properties */
- /* message("Writing runtime parameters..."); */
- h_grp = H5Gcreate1(h_file, "/RuntimePars", 0);
- if(h_grp < 0)
- error("Error while creating runtime parameters group\n");
+ /* Open file */
+ /* message("Opening file '%s'.", fileName); */
+ h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT,H5P_DEFAULT);
+ if(h_file < 0)
+ {
+ error( "Error while opening file '%s'." , fileName );
+ }
- /* Write the relevant information */
- writeAttribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
+ /* Open header to write simulation properties */
+ /* message("Writing runtime parameters..."); */
+ h_grp = H5Gcreate1(h_file, "/RuntimePars", 0);
+ if(h_grp < 0)
+ error("Error while creating runtime parameters group\n");
- /* Close runtime parameters */
- H5Gclose(h_grp);
+ /* Write the relevant information */
+ writeAttribute(h_grp, "PeriodicBoundariesOn", INT, &periodic, 1);
+
+ /* Close runtime parameters */
+ H5Gclose(h_grp);
- /* Open header to write simulation properties */
- /* message("Writing file header..."); */
- h_grp = H5Gcreate1(h_file, "/Header", 0);
- if(h_grp < 0)
- error("Error while creating file header\n");
+ /* Open header to write simulation properties */
+ /* message("Writing file header..."); */
+ h_grp = H5Gcreate1(h_file, "/Header", 0);
+ if(h_grp < 0)
+ error("Error while creating file header\n");
- /* Print the relevant information and print status */
- writeAttribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
- writeAttribute(h_grp, "NumPart_ThisFile", UINT, numParticles, 6);
- writeAttribute(h_grp, "Time", DOUBLE, &e->time, 1);
-
- /* GADGET-2 legacy values */
- writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, 6);
- writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord, 6);
- double MassTable[6] = {0., 0., 0., 0., 0., 0.};
- writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, 6);
- writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, numParticlesHighWord, 6);
- writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
-
- /* Close header */
- H5Gclose(h_grp);
-
- /* Print the SPH parameters */
- writeSPHflavour(h_file);
-
- /* Print the system of Units */
- writeUnitSystem(h_file, us);
+ /* Print the relevant information and print status */
+ writeAttribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3);
+ writeAttribute(h_grp, "Time", DOUBLE, &e->time, 1);
+
+ /* GADGET-2 legacy values */
+ numParticles[0] = (unsigned int) N_total ;
+ writeAttribute(h_grp, "NumPart_ThisFile", UINT, numParticles, 6);
+ writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, 6);
+ numParticlesHighWord[0] = (unsigned int) (N_total >> 32);
+ writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord, 6);
+ double MassTable[6] = {0., 0., 0., 0., 0., 0.};
+ writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, 6);
+ writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy, 6);
+ writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
+
+ /* Close header */
+ H5Gclose(h_grp);
+
+ /* Print the SPH parameters */
+ writeSPHflavour(h_file);
+
+ /* Print the system of Units */
+ writeUnitSystem(h_file, us);
- /* Create SPH particles group */
- /* message("Writing particle arrays..."); */
- h_grp = H5Gcreate1(h_file, "/PartType0", 0);
- if(h_grp < 0)
- error( "Error while creating particle group.\n");
-
- /* Write arrays */
- writeArray(h_grp, fileName, xmfFile, "Coordinates", DOUBLE, N, 3, parts, x, us, UNIT_CONV_LENGTH);
- writeArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, N, 3, parts, v, us, UNIT_CONV_SPEED);
- writeArray(h_grp, fileName, xmfFile, "Masses", FLOAT, N, 1, parts, mass, us, UNIT_CONV_MASS);
- writeArray(h_grp, fileName, xmfFile, "SmoothingLength", FLOAT, N, 1, parts, h, us, UNIT_CONV_LENGTH);
- writeArray(h_grp, fileName, xmfFile, "InternalEnergy", FLOAT, N, 1, parts, u, us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
- writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, parts, id, us, UNIT_CONV_NO_UNITS);
- writeArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N, 1, parts, dt, us, UNIT_CONV_TIME);
- writeArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N, 3, parts, a, us, UNIT_CONV_ACCELERATION);
- writeArray(h_grp, fileName, xmfFile, "Density", FLOAT, N, 1, parts, rho, us, UNIT_CONV_DENSITY);
-
- /* Close particle group */
- H5Gclose(h_grp);
-
- /* Write LXMF file descriptor */
- writeXMFfooter(xmfFile);
+ /* Create SPH particles group */
+ /* message("Writing particle arrays..."); */
+ h_grp = H5Gcreate1(h_file, "/PartType0", 0);
+ if(h_grp < 0)
+ error( "Error while creating particle group.\n");
+
+ /* Prepare the arrays in the file */
+ prepareArray(h_grp, fileName, xmfFile, "Coordinates", DOUBLE, N_total, 3, us, UNIT_CONV_LENGTH);
+ prepareArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, N_total, 3, us, UNIT_CONV_SPEED);
+ prepareArray(h_grp, fileName, xmfFile, "Masses", FLOAT, N_total, 1, us, UNIT_CONV_MASS);
+ prepareArray(h_grp, fileName, xmfFile, "SmoothingLength", FLOAT, N_total, 1, us, UNIT_CONV_LENGTH);
+ prepareArray(h_grp, fileName, xmfFile, "InternalEnergy", FLOAT, N_total, 1, us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ prepareArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N_total, 1, us, UNIT_CONV_NO_UNITS);
+ prepareArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N_total, 1, us, UNIT_CONV_TIME);
+ prepareArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N_total, 3, us, UNIT_CONV_ACCELERATION);
+ prepareArray(h_grp, fileName, xmfFile, "Density", FLOAT, N_total, 1, us, UNIT_CONV_DENSITY);
+
+ /* Close particle group */
+ H5Gclose(h_grp);
- /* message("Done writing particles..."); */
+ /* Close file */
+ H5Fclose(h_file);
+
+ /* Write footer of LXMF file descriptor */
+ writeXMFfooter(xmfFile);
+ }
- /* Close file */
- H5Fclose(h_file);
+
+ /* Now loop over ranks and write the data */
+ for ( rank = 0; rank < mpi_size ; ++ rank ) {
+
+ /* Is it this rank's turn to write ? */
+ if ( rank == mpi_rank ) {
+
+ h_file = H5Fopen(fileName, H5F_ACC_RDWR, H5P_DEFAULT);
+ if(h_file < 0)
+ error( "Error while opening file '%s' on rank %d." , fileName, mpi_rank );
+
+ /* Open SPH particles group */
+ /* message("Reading particle arrays..."); */
+ h_grp = H5Gopen1(h_file, "/PartType0");
+ if(h_grp < 0)
+ error( "Error while opening particle group on rank %d.\n", mpi_rank);
+
+ /* Write arrays */
+ writeArray(h_grp, "Coordinates", DOUBLE, N, 3, N_total, offset, parts, x);
+ writeArray(h_grp, "Velocities", FLOAT, N, 3, N_total, offset, parts, v);
+ writeArray(h_grp, "Masses", FLOAT, N, 1, N_total, offset, parts, mass);
+ writeArray(h_grp, "SmoothingLength", FLOAT, N, 1, N_total, offset, parts, h);
+ writeArray(h_grp, "InternalEnergy", FLOAT, N, 1, N_total, offset, parts, u);
+ writeArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, N_total, offset, parts, id);
+ writeArray(h_grp, "TimeStep", FLOAT, N, 1, N_total, offset, parts, dt);
+ writeArray(h_grp, "Acceleration", FLOAT, N, 3, N_total, offset, parts, a);
+ writeArray(h_grp, "Density", FLOAT, N, 1, N_total, offset, parts, rho);
+
+ /* Close particle group */
+ H5Gclose(h_grp);
+
+ /* Close file */
+ H5Fclose(h_file);
+
+ }
+
+ /* Wait for the read of the reading to complete */
+ MPI_Barrier(comm);
+
+ }
+
+ /* message("Done writing particles..."); */
++outputCount;
}