Commit 320bf527 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

The code can now write output in serial with MPI. Each rank writes its...

The code can now write output in serial with MPI. Each rank writes its particle after the other in the same file.


Former-commit-id: 835922aa8ad8a8f961a238002c820d67f368a3bb
parent 10d3f223
......@@ -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;
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment