Merge branch 'io_fixes' into 'master'

Io fixes

That fixes some issues #52 that appeared on the BlueGene when reading more than 1>>32 particles in parallel or serial mode.

I have also modified the uniform box script to be much faster and lightweight. 

See merge request !48

examples/UniformBox/makeIC.py

@@ -19,6 +19,7 @@
  ##############################################################################
 
 import h5py
+import sys
 from numpy import *
 
 # Generates a swift IC file containing a cartesian distribution of particles
@@ -27,7 +28,7 @@ from numpy import *
 # Parameters
 periodic= 1           # 1 For periodic box
 boxSize = 1.
-L = 50           # Number of particles along one axis
+L = int(sys.argv[1])  # Number of particles along one axis
 rho = 2.              # Density
 P = 1.                # Pressure
 gamma = 5./3.         # Gas adiabatic index
@@ -39,34 +40,6 @@ numPart = L**3
 mass = boxSize**3 * rho / numPart
 internalEnergy = P / ((gamma - 1.)*rho)
 
-#Generate particles
-coords = zeros((numPart, 3))
-v      = zeros((numPart, 3))
-m      = zeros((numPart, 1))
-h      = zeros((numPart, 1))
-u      = zeros((numPart, 1))
-ids    = zeros((numPart, 1), dtype='L')
-
-for i in range(L):
-    for j in range(L):
-        for k in range(L):
-            index = i*L*L + j*L + k
-            x = i * boxSize / L + boxSize / (2*L)
-            y = j * boxSize / L + boxSize / (2*L)
-            z = k * boxSize / L + boxSize / (2*L)
-            coords[index,0] = x
-            coords[index,1] = y
-            coords[index,2] = z
-            v[index,0] = 0.
-            v[index,1] = 0.
-            v[index,2] = 0.
-            m[index] = mass
-            h[index] = 2.251 * boxSize / L
-            u[index] = internalEnergy
-            ids[index] = index
-            
-            
-
 #--------------------------------------------------
 
 #File
@@ -90,17 +63,39 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
 
 #Particle group
 grp = file.create_group("/PartType0")
-ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
-ds[()] = coords
+
+v  = zeros((numPart, 3))
 ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
 ds[()] = v
+v = zeros(1)
+
+m = full((numPart, 1), mass)
 ds = grp.create_dataset('Masses', (numPart,1), 'f')
 ds[()] = m
+m = zeros(1)
+
+h = full((numPart, 1), 2.251 * boxSize / L)
 ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
 ds[()] = h
+h = zeros(1)
+
+u = full((numPart, 1), internalEnergy)
 ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
 ds[()] = u
+u = zeros(1)
+
+
+ids = linspace(0, numPart, numPart, endpoint=False, dtype='L').reshape((numPart,1))
 ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
 ds[()] = ids
+x      = ids % L;
+y      = ((ids - x) / L) % L;
+z      = (ids - x - L * y) / L**2;
+coords = zeros((numPart, 3))
+coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
+coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
+coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = coords
 
 file.close()

src/parallel_io.c

@@ -110,9 +110,9 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
     offsets[0] = offset;
     offsets[1] = 0;
   } else {
-    rank = 1;
+    rank = 2;
     shape[0] = N;
-    shape[1] = 0;
+    shape[1] = 1;
     offsets[0] = offset;
     offsets[1] = 0;
   }
@@ -238,13 +238,13 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
   dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
   dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];
 
+  /* message("Found %d particles in a %speriodic box of size [%f %f %f].",  */
+  /* 	 N_total, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
+
   /* Divide the particles among the tasks. */
   offset = mpi_rank * N_total / mpi_size;
   *N = (mpi_rank + 1) * N_total / mpi_size - offset;
 
-  /* message("Found %d particles in a %speriodic box of size [%f %f %f].",  */
-  /* 	 *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
-
   /* Close header */
   H5Gclose(h_grp);
 

src/serial_io.c

@@ -87,13 +87,15 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
     if (importance == COMPULSORY) {
       error("Compulsory data set '%s' not present in the file.", name);
     } else {
-      for (i = 0; i < N; ++i) memset(part_c + i * partSize, 0, copySize);
+      for (i = 0; i < N; ++i)
+        memset(part_c + i * partSize, 0, copySize);
       return;
     }
   }
 
-  /* message( "Reading %s '%s' array...", importance == COMPULSORY ?
-   * "compulsory": "optional  ", name); */
+  /* message( "Reading %s '%s' array...", importance == COMPULSORY ? */
+  /* 	   "compulsory": "optional  ", name); */
+  /* fflush(stdout); */
 
   /* Open data space */
   h_data = H5Dopen1(grp, name);
@@ -117,9 +119,9 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
     offsets[0] = offset;
     offsets[1] = 0;
   } else {
-    rank = 1;
+    rank = 2;
     shape[0] = N;
-    shape[1] = 0;
+    shape[1] = 1;
     offsets[0] = offset;
     offsets[1] = 0;
   }
@@ -131,6 +133,20 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
   h_filespace = H5Dget_space(h_data);
   H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape, NULL);
 
+  /* int rank_memspace = H5Sget_simple_extent_ndims(h_memspace); */
+  /* int rank_filespace = H5Sget_simple_extent_ndims(h_filespace); */
+
+  /* message("Memspace rank: %d", rank_memspace); */
+  /* message("Filespace rank: %d", rank_filespace); */
+  /* fflush(stdout); */
+
+  /* hsize_t dims_memspace[2], max_dims_memspace[2]; */
+  /* hsize_t dims_filespace[2], max_dims_filespace[2]; */
+
+  /* H5Sget_simple_extent_dims(h_memspace, dims_memspace, max_dims_memspace); */
+  /* H5Sget_simple_extent_dims(h_filespace, dims_filespace, max_dims_filespace);
+   */
+
   /* Read HDF5 dataspace in temporary buffer */
   /* Dirty version that happens to work for vectors but should be improved */
   /* Using HDF5 dataspaces would be better */
@@ -196,11 +212,11 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, int* N,
                     int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
                     MPI_Info info) {
   hid_t h_file = 0, h_grp = 0;
-  double boxSize[3] = {
-      0.0, -1.0, -1.0}; /* GADGET has only cubic boxes (in cosmological mode) */
-  int numParticles[6] = {
-      0}; /* GADGET has 6 particle types. We only keep the type 0*/
-  int numParticles_highWord[6] = {0};
+  double boxSize[3] = { 0.0, -1.0, -1.0 };
+      /* GADGET has only cubic boxes (in cosmological mode) */
+  int numParticles[6] = { 0 };
+      /* GADGET has 6 particle types. We only keep the type 0*/
+  int numParticles_highWord[6] = { 0 };
   long long offset = 0;
   long long N_total = 0;
   int rank;
@@ -235,14 +251,17 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, int* N,
     readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
     readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord);
 
-    N_total = ((long long)numParticles[0]) +
-              ((long long)numParticles_highWord[0] << 32);
+    N_total = ((long long) numParticles[0]) +
+              ((long long) numParticles_highWord[0] << 32);
     dim[0] = boxSize[0];
     dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
     dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];
 
-    /* message("Found %d particles in a %speriodic box of size [%f %f %f].",  */
-    /* 	 *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
+    /* message("Found %lld particles in a %speriodic box of size [%f %f %f].",
+     */
+    /* 	    N_total, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
+
+    fflush(stdout);
 
     /* Close header */
     H5Gclose(h_grp);
@@ -264,8 +283,8 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, int* N,
   if (posix_memalign((void*)parts, part_align, (*N) * sizeof(struct part)) != 0)
     error("Error while allocating memory for particles");
   bzero(*parts, *N * sizeof(struct part));
-  /* message("Allocated %8.2f MB for particles.", *N * sizeof(struct part) /
-   * (1024.*1024.)); */
+  /* message("Allocated %8.2f MB for particles.", *N * sizeof(struct part) / */
+  /* 	  (1024.*1024.)); */
 
   /* Now loop over ranks and read the data */
   for (rank = 0; rank < mpi_size; ++rank) {
@@ -499,9 +518,9 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
   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};
+  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;
@@ -516,7 +535,7 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
 
   /* 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;
+  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);
@@ -524,8 +543,8 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
     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;
+  N_total = (long long) N_total_d;
+  offset = (long long) offset_d;
 
   /* Do common stuff first */
   if (mpi_rank == 0) {
@@ -568,13 +587,13 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
     writeAttribute(h_grp, "Time", DOUBLE, &dblTime, 1);
 
     /* GADGET-2 legacy values */
-    numParticles[0] = (unsigned int)N_total;
+    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.};
+    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);

src/single_io.c

@@ -83,7 +83,8 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
       /* message("Optional data set '%s' not present. Zeroing this particle
        * field...", name);	   */
 
-      for (i = 0; i < N; ++i) memset(part_c + i * partSize, 0, copySize);
+      for (i = 0; i < N; ++i)
+        memset(part_c + i * partSize, 0, copySize);
 
       return;
     }
@@ -166,10 +167,10 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
 void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
                     int* periodic) {
   hid_t h_file = 0, h_grp = 0;
-  double boxSize[3] = {
-      0.0, -1.0, -1.0}; /* GADGET has only cubic boxes (in cosmological mode) */
-  int numParticles[6] = {
-      0}; /* GADGET has 6 particle types. We only keep the type 0*/
+  double boxSize[3] = { 0.0, -1.0, -1.0 };
+      /* GADGET has only cubic boxes (in cosmological mode) */
+  int numParticles[6] = { 0 };
+      /* GADGET has 6 particle types. We only keep the type 0*/
 
   /* Open file */
   /* message("Opening file '%s' as IC.", fileName); */
@@ -384,8 +385,8 @@ void write_output_single(struct engine* e, struct UnitSystem* us) {
   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};
+  int numParticles[6] = { N, 0 };
+  int numParticlesHighWord[6] = { 0 };
   int numFiles = 1;
   struct part* parts = e->s->parts;
   FILE* xmfFile = 0;
@@ -437,7 +438,7 @@ void write_output_single(struct engine* e, struct UnitSystem* us) {
   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.};
+  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);