Merge branch 'InitFromHDF5' into slabUpdate

Author: Martin D. Weinberg

src/Component.H

@@ -184,6 +184,43 @@ struct loadb_datum
   specifying alternative parameters, using an old-style PSP file, or
   starting a new job with a previous output from another simulation)
   </ol>
+
+  <br> <b>Input data format</b> Component now supports reading from
+  HDF5 files on input.  The scheme for HDF5 input files is:
+  @verbatim
+  /
+  ├── Attributes
+  │ ├── num_particles (int)
+  │ ├── num_aux_ints (int)
+  │ └── num_aux_floats (int)
+  └── particles/ (Group)
+  ├── m (Dataset)
+  ├── x, y, z (Datasets)
+  ├── u, v, w (Datasets)
+  ├── index (optional Dataset)
+  ├── aux_int_0, aux_int_1, ... (Datasets)
+  └── aux_float_0, aux_float_1, ... (Datasets)
+  @endverbatim
+
+  @note
+  <ol>
+
+  <li> The <code>index</code> dataset is optional.  If it is not
+  present, the particles are assigned indices in the order they are read.
+
+  <li> The <code>aux_int_*</code> and <code>aux_float_*</code>
+  datasets are optional.  If they are not present, the corresponding
+  vectors are empty.  The number of auxiliary integer and float
+  datasets is determined by the <code>num_aux_ints</code> and
+  <code>num_aux_floats</code> attributes, respectively.
+
+  <li> The datatype for the <code>m</code>, <code>x</code>,
+  <code>y</code>, <code>z</code>, <code>u</code>, <code>v</code>,
+  <code>w</code> and <code>aux_float_*</code> datasets maybe float or
+  double (i.e. float32 or float64) but they all must are assumed to
+  have the same datatype as <code>m</code>.
+
+  </ol>
 */
 class Component
 {

src/Component.cc

@@ -1572,6 +1572,8 @@ void Component::read_bodies_and_distribute_hdf5(void)
   // Phase-space data
   //
   FloatData m, x, y, z, u, v, w;
+  std::vector<unsigned long> index;
+  bool has_index = particles_group.exist("index");
   std::vector<std::vector<int>> aux_ints;
   std::vector<FloatData> aux_floats;
 
@@ -1589,11 +1591,9 @@ void Component::read_bodies_and_distribute_hdf5(void)
       throw std::runtime_error("Component ERROR: dataset <" + name + "> shorter than expected");
     }
 
-    // Sanity check to ensure we don't read beyond the dataset
     size_t current = std::min(batch, total - offset);
-    // Wrap returned vector in FloatData
-    if (precision == 4) {
-      return FloatData(dataset.select({offset}, {current}).read<std::vector<float>>());
+    if (current != batch) {
+      throw std::runtime_error("Component ERROR: dataset <" + name + "> shorter than expected");
     }
     else if (precision == 8) {
       return FloatData(dataset.select({offset}, {current}).read<std::vector<double>>());
@@ -1627,6 +1627,10 @@ void Component::read_bodies_and_distribute_hdf5(void)
     size_t offset = 0;
 
     // Read Node 0's PSP fields
+    if (has_index) {
+      index = particles_group.getDataSet("index")
+        .select({offset}, {nbodies_table[0]}).read<std::vector<unsigned long>>();
+    }
     m = read_dataset("m", offset, nbodies_table[0]);
 
     x = read_dataset("x", offset, nbodies_table[0]);
@@ -1675,8 +1679,8 @@ void Component::read_bodies_and_distribute_hdf5(void)
 	part->dattrib[j] = aux_floats[j][i];
       }
 
-      // Make a sequential index for the particle
-      part->indx = i + 1;
+      // Set particle index from optional dataset, otherwise sequential
+      part->indx = has_index ? index[i] : i + 1;
 
       // Get the radius
       r2 = 0.0;
@@ -1699,6 +1703,10 @@ void Component::read_bodies_and_distribute_hdf5(void)
     for (int n=1; n<numprocs; n++) {
 
       // Read core PSP fields
+      if (has_index) {
+        index = particles_group.getDataSet("index")
+          .select({offset}, {nbodies_table[n]}).read<std::vector<unsigned long>>();
+      }
       m = read_dataset("m", offset, nbodies_table[n]);
 
       x = read_dataset("x", offset, nbodies_table[n]);
@@ -1732,8 +1740,8 @@ void Component::read_bodies_and_distribute_hdf5(void)
 
 	PartPtr part = std::make_shared<Particle>(niattrib, ndattrib);
 
-	// Make a sequential index for the particle
-	part->indx = nbodies_index[n-1] + 1 + icount;
+	// Set particle index from optional dataset, otherwise sequential
+	part->indx = has_index ? index[icount] : nbodies_index[n-1] + 1 + icount;
 
 	part->mass   = m[icount];