Add HDF5 output support

One HDF5 file per variable, and one extensible dataset per radius per mode.  HDF5 is required only optionally, so this thorn can be compiled without it.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Multipole/trunk@74 4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843

6 files changed, 165 insertions, 8 deletions

diff --git a/configuration.ccl b/configuration.ccl
new file mode 100644
index 0000000..b3f112e
--- /dev/null
+++ b/configuration.ccl
@@ -0,0 +1,3 @@
+OPTIONAL HDF5
+{
+}
diff --git a/interface.ccl b/interface.ccl
index 463a8ce..a8eda59 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,6 +2,11 @@
 implements: multipole
 inherits:   Grid 
 
+CCTK_INT FUNCTION IO_TruncateOutputFiles \
+  (CCTK_POINTER_TO_CONST IN GH)
+
+REQUIRES FUNCTION IO_TruncateOutputFiles
+
 public:
 
 CCTK_REAL harmonics type=GF timelevels=1
diff --git a/param.ccl b/param.ccl
index 2f089dc..9af3650 100644
--- a/param.ccl
+++ b/param.ccl
@@ -137,3 +137,16 @@ CCTK_INT m_mode "all modes: Up to which m mode to calculate/ specific mode: whic
 {
   -100:* :: "Deprecated"
 } -100
+
+CCTK_BOOLEAN output_ascii "Output a simple ASCII file for each mode at each radius"
+{
+} "yes"
+
+CCTK_BOOLEAN output_hdf5 "Output an HDF5 file for each variable containing one dataset per mode at each radius"
+{
+} "no"
+
+CCTK_INT hdf5_chunk_size "How many iterations to preallocate in extensible HDF5 datasets"
+{
+  1:* :: "Any integer"
+} 200
diff --git a/src/multipole.cc b/src/multipole.cc
index 975d589..8a6a6ea 100644
--- a/src/multipole.cc
+++ b/src/multipole.cc
@@ -97,15 +97,28 @@ static void parse_variables_string(const string &var_string, variable_desc v[max
 static void output_mode(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad, 
                      CCTK_REAL real_lm, CCTK_REAL imag_lm, int l, int m)
 {
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
 
-  if (CCTK_MyProc(cctkGH) == 0)
+  if (output_ascii)
+  {
+    if (CCTK_MyProc(cctkGH) == 0)
+    {
+      ostringstream name;
+      name << "mp_" << v->name << "_l" << l << "_m" << m <<
+        "_r" << setiosflags(ios::fixed) << setprecision(2) <<  rad << ".asc";
+      Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name.str(), real_lm, imag_lm);
+    }
+  }
+  if (output_hdf5)
   {
-    ostringstream name;
-    name << "mp_" << v->name << "_l" << l << "_m" << m <<
-      "_r" << setiosflags(ios::fixed) << setprecision(2) <<  rad << ".asc";
-    Multipole_OutputComplexToFile(CCTK_PASS_CTOC, name.str(), real_lm, imag_lm);
+    if (CCTK_MyProc(cctkGH) == 0)
+    {
+      ostringstream h5datasetname;
+      h5datasetname << "l" << l << "_m" << m << "_r" << setiosflags(ios::fixed) << setprecision(2) << rad;
+      Multipole_OutputComplexToH5File(CCTK_PASS_CTOC, "mp_" + v->name + ".h5", h5datasetname.str(),
+                                      real_lm, imag_lm);
+    }
   }
 }
 
@@ -116,7 +129,7 @@ static void output_1D(CCTK_ARGUMENTS, const variable_desc *v, CCTK_REAL rad,
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
 
-  if (CCTK_MyProc(cctkGH) == 0)
+  if (CCTK_MyProc(cctkGH) == 0 && output_ascii)
   {
     if (out_1d_every != 0 && cctk_iteration % out_1d_every == 0)
     {
diff --git a/src/utils.cc b/src/utils.cc
index 168acc2..7140718 100644
--- a/src/utils.cc
+++ b/src/utils.cc
@@ -2,12 +2,25 @@
 #include <string.h>
 #include <math.h>
 #include <string>
+#include <sys/stat.h>
+#include <errno.h>
 #include <iostream>
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
+#ifdef HAVE_CAPABILITY_HDF5
+// We currently support the HDF5 1.6 API (and when using 1.8 the
+// compatibility mode introduced by H5_USE_16_API).  Several machines
+// in SimFactory use HDF5 1.6, so we cannot drop support for it.  It
+// seems it is hard to support both the 1.6 and 1.8 API
+// simultaneously; for example H5Fopen takes a different number of
+// arguments in the two versions.
+#define H5_USE_16_API
+#include <hdf5.h>
+#endif
+
 #include "utils.hh"
 #include "integrate.hh"
 
@@ -133,6 +146,113 @@ void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL
   }
 }
 
+////////////////////////////////////////////////////////////////////////
+// HDF5 complex output
+////////////////////////////////////////////////////////////////////////
+
+#ifdef HAVE_CAPABILITY_HDF5
+
+bool dataset_exists(hid_t file, const string &dataset_name)
+{
+  // To test whether a dataset exists, the recommended way in API 1.6
+  // is to use H5Gget_objinfo, but this prints an error to stderr if
+  // the dataset does not exist.  We explicitly avoid this by wrapping
+  // the call in H5E_BEGIN_TRY/H5E_END_TRY statements.  In 1.8,
+  // H5Gget_objinfo is deprecated, and H5Lexists does the job.  See
+  // http://www.mail-archive.com/hdf-forum@hdfgroup.org/msg00125.html
+
+  #if 1
+  bool exists;
+  H5E_BEGIN_TRY
+  {
+    exists = H5Gget_objinfo(file, dataset_name.c_str(), 1, NULL) >= 0;
+  }
+  H5E_END_TRY;
+  return exists;
+  #else
+  return H5Lexists(file, dataset_name.c_str(), H5P_DEFAULT);
+  #endif
+}
+
+void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &basename, const string &datasetname,
+                                     CCTK_REAL redata, CCTK_REAL imdata)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  string output_name = out_dir + string("/") + basename;
+  int status = 0;
+  static int first_time = true;
+
+  hid_t file =
+    (first_time && IO_TruncateOutputFiles(cctkGH)) ?
+    H5Fcreate(output_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) :
+    H5Fopen(output_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
+
+  first_time = false;
+
+  hid_t dataset = -1;
+
+  if (dataset_exists(file, datasetname))
+  {
+    dataset = H5Dopen(file, datasetname.c_str());
+  }
+  else
+  {
+    hsize_t dims[2] = {0,3};
+    hsize_t maxdims[2] = {H5S_UNLIMITED, 3};
+    hid_t dataspace = H5Screate_simple(2, dims, maxdims);
+
+    hid_t cparms = -1;
+    hsize_t chunk_dims[2] = {hdf5_chunk_size,3};
+    cparms = H5Pcreate (H5P_DATASET_CREATE);
+    status = H5Pset_chunk(cparms, 2, chunk_dims);
+
+    dataset = H5Dcreate(file, datasetname.c_str(), H5T_NATIVE_DOUBLE, dataspace, cparms);
+    H5Pclose(cparms);
+  }
+
+  hid_t filespace = H5Dget_space (dataset);
+  hsize_t filedims[2];
+  hsize_t maxdims[2];
+  status = H5Sget_simple_extent_dims(filespace, filedims, maxdims);
+
+  filedims[0] += 1;
+  hsize_t size[2] = {filedims[0], filedims[1]};
+  status = H5Dextend (dataset, size);
+  H5Sclose(filespace);
+
+  /* Select a hyperslab  */
+  hsize_t offset[2] = {filedims[0]-1, 0};
+  hsize_t dims2[2] = {1,3};
+  filespace = H5Dget_space (dataset);
+  status = H5Sselect_hyperslab (filespace, H5S_SELECT_SET, offset, NULL,
+                                dims2, NULL);
+
+  CCTK_REAL data[] = {cctk_time, redata, imdata};
+
+  hid_t memdataspace = H5Screate_simple(2, dims2, NULL);
+
+  /* Write the data to the hyperslab  */
+  status = H5Dwrite (dataset, H5T_NATIVE_DOUBLE, memdataspace, filespace,
+                     H5P_DEFAULT, data);
+
+  H5Dclose(dataset);
+  H5Sclose(filespace);
+  H5Sclose(memdataspace);
+
+  status = H5Fclose(file);
+}
+
+#else
+
+void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &basename, const string &datasetname,
+                                     CCTK_REAL redata, CCTK_REAL imdata)
+{
+  CCTK_WARN(0,"HDF5 output has been requested but Cactus has been compiled without HDF5 support");
+}
+
+#endif
 
 ////////////////////////////////////////////////////////////////////////
 // Coordinates
diff --git a/src/utils.hh b/src/utils.hh
index 51c76a7..7d10e8b 100644
--- a/src/utils.hh
+++ b/src/utils.hh
@@ -20,6 +20,9 @@ void Multipole_Output1D(CCTK_ARGUMENTS, const string &name, int array_size,
 
 void Multipole_OutputComplexToFile(CCTK_ARGUMENTS, const string &name, CCTK_REAL redata, CCTK_REAL imdata);
 
+void Multipole_OutputComplexToH5File(CCTK_ARGUMENTS, const string &filename, const string &dataset,
+                                     CCTK_REAL redata, CCTK_REAL imdata);
+
 void Multipole_CoordSetup(int ntheta, int nphi, 
                           CCTK_REAL xhat[], CCTK_REAL yhat[], 
                           CCTK_REAL zhat[], CCTK_REAL th[],